MOAB: Mesh Oriented datABase  (version 5.4.0)
QualityAssessor.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025 
00026   ***************************************************************** */
00027 /*!
00028   \file   QualityAssessor.cpp
00029   \brief  Member function of the MBMesquite::QualityAssessor class
00030 
00031   \author Thomas Leurent
00032   \date   2002-05-23
00033 */
00034 
00035 #include "QualityAssessor.hpp"
00036 #include "QualityMetric.hpp"
00037 #include "TMPQualityMetric.hpp"
00038 #include "ElementMaxQM.hpp"
00039 #include "ElementAvgQM.hpp"
00040 #include "PatchData.hpp"
00041 #include "MsqMeshEntity.hpp"
00042 #include "MsqVertex.hpp"
00043 #include "MsqDebug.hpp"
00044 #include "MeshInterface.hpp"
00045 #include "VertexPatches.hpp"
00046 #include "ElementPatches.hpp"
00047 #include "ParallelMeshInterface.hpp"
00048 #include "ParallelHelperInterface.hpp"
00049 
00050 #include <list>
00051 #include <vector>
00052 #include <iostream>
00053 #include <fstream>
00054 #include <iomanip>
00055 #include <set>
00056 #include <cmath>
00057 
00058 #ifdef HAVE_SYS_IOCTL_H
00059 #include <sys/ioctl.h>
00060 #endif
00061 #ifdef HAVE_UNISTD_H
00062 #include <unistd.h>
00063 #endif
00064 #ifdef HAVE_TERMIOS_H
00065 #include <termios.h>
00066 #endif
00067 
00068 namespace MBMesquite
00069 {
00070 
00071 const char* default_name( bool free_only )
00072 {
00073     static const char all_name[]  = "QualityAssessor";
00074     static const char free_name[] = "QualityAssessor(free only)";
00075     return free_only ? free_name : all_name;
00076 }
00077 
00078 QualityAssessor::QualityAssessor( bool p_print_summary,
00079                                   bool free_only,
00080                                   const char* inverted_tag_name,
00081                                   std::string p_name )
00082     : myData( new Data ), qualityAssessorName( p_name ), outputStream( std::cout ), printSummary( p_print_summary ),
00083       skipFixedSamples( free_only )
00084 {
00085     if( inverted_tag_name ) tag_inverted_elements( inverted_tag_name );
00086 
00087     if( qualityAssessorName.empty() ) qualityAssessorName = default_name( free_only );
00088 
00089     for( int i = POLYGON; i <= MIXED; i++ )
00090         elementTypeCount[i - POLYGON] = 0;
00091 }
00092 
00093 QualityAssessor::QualityAssessor( std::ostream& stream,
00094                                   bool free_only,
00095                                   const char* inverted_tag_name,
00096                                   std::string name )
00097     : myData( new Data ), qualityAssessorName( name ), outputStream( stream ), printSummary( true ),
00098       skipFixedSamples( free_only )
00099 {
00100     if( inverted_tag_name ) tag_inverted_elements( inverted_tag_name );
00101 
00102     if( qualityAssessorName.empty() ) qualityAssessorName = default_name( free_only );
00103 
00104     for( int i = POLYGON; i <= MIXED; i++ )
00105         elementTypeCount[i - POLYGON] = 0;
00106 }
00107 
00108 QualityAssessor::QualityAssessor( std::ostream& output_stream,
00109                                   QualityMetric* metric,
00110                                   int histogram_intervals,
00111                                   double power_mean,
00112                                   bool free_only,
00113                                   const char* metric_value_tag_name,
00114                                   const char* inverted_tag_name,
00115                                   std::string name )
00116     : myData( new Data ), qualityAssessorName( name ), outputStream( output_stream ), printSummary( true ),
00117       skipFixedSamples( free_only )
00118 {
00119     if( inverted_tag_name ) tag_inverted_elements( inverted_tag_name );
00120 
00121     if( qualityAssessorName.empty() ) qualityAssessorName = default_name( free_only );
00122 
00123     set_stopping_assessment( metric, histogram_intervals, power_mean, metric_value_tag_name );
00124 
00125     for( int i = POLYGON; i <= MIXED; i++ )
00126         elementTypeCount[i - POLYGON] = 0;
00127 }
00128 
00129 QualityAssessor::QualityAssessor( QualityMetric* metric,
00130                                   int histogram_intervals,
00131                                   double power_mean,
00132                                   bool free_only,
00133                                   const char* metric_value_tag_name,
00134                                   bool p_print_summary,
00135                                   const char* inverted_tag_name,
00136                                   std::string name )
00137     : myData( new Data ), qualityAssessorName( name ), outputStream( std::cout ), printSummary( p_print_summary ),
00138       skipFixedSamples( free_only )
00139 {
00140     if( inverted_tag_name ) tag_inverted_elements( inverted_tag_name );
00141 
00142     if( qualityAssessorName.empty() ) qualityAssessorName = default_name( free_only );
00143 
00144     set_stopping_assessment( metric, histogram_intervals, power_mean, metric_value_tag_name );
00145 
00146     for( int i = POLYGON; i <= MIXED; i++ )
00147         elementTypeCount[i - POLYGON] = 0;
00148 }
00149 
00150 QualityAssessor::QualityAssessor( const QualityAssessor& copy )
00151     : myData( copy.myData ), qualityAssessorName( copy.qualityAssessorName ), assessList( copy.assessList ),
00152       outputStream( copy.outputStream ), printSummary( copy.printSummary ), invertedTagName( copy.invertedTagName ),
00153       fixedTagName( copy.fixedTagName ), skipFixedSamples( copy.skipFixedSamples )
00154 
00155 {
00156     list_type::iterator iter;
00157     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00158         ( *iter )->referenceCount++;
00159     myData->referenceCount++;
00160 
00161     for( int i = POLYGON; i <= MIXED; i++ )
00162         elementTypeCount[i - POLYGON] = copy.elementTypeCount[i - POLYGON];
00163 }
00164 
00165 QualityAssessor& QualityAssessor::operator=( const QualityAssessor& copy )
00166 {
00167     list_type::iterator iter;
00168     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00169     {
00170         Assessor* assessor = *iter;
00171         if( 0 == --assessor->referenceCount ) delete assessor;
00172     }
00173 
00174     myData              = copy.myData;
00175     qualityAssessorName = copy.qualityAssessorName;
00176     assessList          = copy.assessList;
00177     printSummary        = copy.printSummary;
00178     invertedTagName     = copy.invertedTagName;
00179     fixedTagName        = copy.fixedTagName;
00180     skipFixedSamples    = copy.skipFixedSamples;
00181 
00182     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00183         ( *iter )->referenceCount++;
00184     myData->referenceCount++;
00185 
00186     for( int i = POLYGON; i <= MIXED; i++ )
00187         elementTypeCount[i - POLYGON] = copy.elementTypeCount[i - POLYGON];
00188 
00189     return *this;
00190 }
00191 
00192 QualityAssessor::~QualityAssessor()
00193 {
00194     list_type::iterator iter;
00195     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00196     {
00197         Assessor* assessor = *iter;
00198         if( 0 == --assessor->referenceCount ) delete assessor;
00199     }
00200     if( 0 == --myData->referenceCount ) delete myData;
00201 }
00202 
00203 double QualityAssessor::Assessor::get_average() const
00204 {
00205     return count ? sum / count : 0;
00206 }
00207 
00208 double QualityAssessor::Assessor::get_rms() const
00209 {
00210     return count ? sqrt( sqrSum / count ) : 0;
00211 }
00212 
00213 double QualityAssessor::Assessor::get_stddev() const
00214 {
00215     double sqr = count ? sqrSum / count - sum * sum / ( (double)count * count ) : 0;
00216     return sqr < 0 ? 0 : sqrt( sqr );
00217 }
00218 
00219 double QualityAssessor::Assessor::get_power_mean() const
00220 {
00221     return ( count && pMean ) ? pow( pSum / count, 1 / pMean ) : 0;
00222 }
00223 
00224 bool QualityAssessor::get_inverted_element_count( int& inverted_elems, int& inverted_samples, MsqError& err )
00225 {
00226     if( myData->invertedElementCount == -1 )
00227     {
00228         MSQ_SETERR( err )
00229         ( "Number of inverted elements has not yet been calculated.", MsqError::INVALID_STATE );
00230         return false;
00231     }
00232     inverted_elems   = myData->invertedElementCount;
00233     inverted_samples = myData->invertedSampleCount;
00234     return true;
00235 }
00236 
00237 void QualityAssessor::add_quality_assessment( QualityMetric* metric,
00238                                               int histogram_intervals,
00239                                               double power_mean,
00240                                               const char* tag_name,
00241                                               const char* label )
00242 {
00243     list_type::iterator i = find_or_add( metric, label );
00244     ( *i )->pMean         = power_mean;
00245     if( histogram_intervals > 0 )
00246         ( *i )->histogram.resize( histogram_intervals + 2 );
00247     else
00248         ( *i )->histogram.clear();
00249     ( *i )->haveHistRange = false;
00250     if( !tag_name )
00251         ( *i )->tagName.clear();
00252     else
00253         ( *i )->tagName = tag_name;
00254 }
00255 
00256 QualityAssessor::list_type::iterator QualityAssessor::find_or_add( QualityMetric* qm, const char* label )
00257 {
00258     list_type::iterator iter;
00259 
00260     // If metric is already in list, find it
00261     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00262         if( ( *iter )->qualMetric == qm ) break;
00263 
00264     // If metric not found in list, add it
00265     if( iter == assessList.end() )
00266     {
00267         Assessor* new_assessor       = new Assessor( qm, label );
00268         new_assessor->referenceCount = 1;
00269         if( qm->get_metric_type() == QualityMetric::VERTEX_BASED )
00270         {
00271             assessList.push_back( new_assessor );
00272             iter = --assessList.end();
00273         }
00274         else
00275         {
00276             assessList.push_front( new_assessor );
00277             iter = assessList.begin();
00278         }
00279     }
00280 
00281     return iter;
00282 }
00283 
00284 QualityAssessor::list_type::iterator QualityAssessor::find_stopping_assessment()
00285 {
00286     list_type::iterator iter;
00287     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00288         if( ( *iter )->stopping_function() ) break;
00289     return iter;
00290 }
00291 
00292 /*!Sets which QualityMetric and QAFunction
00293 combination is used to determine the value return from assess_mesh_quality().
00294 It first ensures that the inputed QAFunction was not HISTOGRAM.  It then
00295 calls add_quality_assessment with the given QualityMetric and QAFunction,
00296 to ensure that this combination will be computed.  Finally, it sets
00297 the stoppingMetric pointer and the stoppingFunction data members.
00298 \param qm Pointer to QualityMetric.
00299 \param func (QAFUNCTION) Wrapper function for qm (e.g. MINIMUM, MAXIMUM,...).
00300     */
00301 void QualityAssessor::set_stopping_assessment( QualityMetric* metric,
00302                                                int histogram_intervals,
00303                                                double power_mean,
00304                                                const char* tag_name,
00305                                                const char* label )
00306 {
00307     list_type::iterator i = find_stopping_assessment();
00308     if( i != assessList.end() ) ( *i )->set_stopping_function( false );
00309 
00310     i             = find_or_add( metric, label );
00311     ( *i )->pMean = power_mean;
00312     if( histogram_intervals > 0 )
00313         ( *i )->histogram.resize( histogram_intervals + 2 );
00314     else
00315         ( *i )->histogram.clear();
00316     ( *i )->haveHistRange = false;
00317     if( !tag_name )
00318         ( *i )->tagName.clear();
00319     else
00320         ( *i )->tagName = tag_name;
00321     ( *i )->set_stopping_function( true );
00322 }
00323 
00324 /*!
00325 Checks first to see if the QualityMetric, qm, has been added to this
00326 QualityAssessor, and if it has not, adds it.  It then adds HISTOGRAM as a
00327 QAFunciton for that metric.  It then sets the minimum and maximum values
00328 for the histogram.
00329 \param qm Pointer to the QualityMetric to be used in histogram.
00330 \param min_val (double) Minimum range of histogram.
00331 \param max_val (double) Maximum range of histogram.
00332 \param intervals Number of histogram intervals
00333     */
00334 void QualityAssessor::add_histogram_assessment( QualityMetric* metric,
00335                                                 double min_val,
00336                                                 double max_val,
00337                                                 int intervals,
00338                                                 double power_mean,
00339                                                 const char* tag_name,
00340                                                 const char* label )
00341 {
00342     if( intervals < 1 ) intervals = 1;
00343     list_type::iterator i = find_or_add( metric, label );
00344     ( *i )->pMean         = power_mean;
00345     ( *i )->histMin       = min_val;
00346     ( *i )->histMax       = max_val;
00347     ( *i )->haveHistRange = min_val < max_val;
00348     ( *i )->histogram.resize( intervals + 2 );
00349     if( !tag_name )
00350         ( *i )->tagName.clear();
00351     else
00352         ( *i )->tagName = tag_name;
00353 }
00354 
00355 TagHandle QualityAssessor::get_tag( Mesh* mesh, std::string name, Mesh::TagType type, unsigned size, MsqError& err )
00356 {
00357     TagHandle tag = mesh->tag_get( name, err );
00358     if( !err )
00359     {
00360         Mesh::TagType exist_type;
00361         std::string junk;
00362         unsigned exist_size;
00363         mesh->tag_properties( tag, junk, exist_type, exist_size, err );
00364         MSQ_ERRZERO( err );
00365         if( type != exist_type || size != exist_size )
00366         {
00367             MSQ_SETERR( err )
00368             ( MsqError::TAG_ALREADY_EXISTS, "Tag \"%s\" exists with incorrect type or length.", name.c_str() );
00369         }
00370     }
00371     else if( err.error_code() == MsqError::TAG_NOT_FOUND )
00372     {
00373         err.clear();
00374         tag = mesh->tag_create( name, type, size, 0, err );
00375         MSQ_ERRZERO( err );
00376     }
00377     else
00378     {
00379         MSQ_ERRZERO( err );
00380     }
00381 
00382     return tag;
00383 }
00384 
00385 void QualityAssessor::initialize_queue( MeshDomainAssoc* mesh_and_domain, const Settings* settings, MsqError& err )
00386 {
00387     for( list_type::iterator i = assessList.begin(); i != assessList.end(); ++i )
00388     {
00389         ( *i )->get_metric()->initialize_queue( mesh_and_domain, settings, err );MSQ_ERRRTN( err );
00390     }
00391 }
00392 
00393 /*!
00394   Computes the quality data for a given
00395   MeshSet, ms. What quality information is calculated, depends
00396   on what has been requested through the use of the QualityAssessor
00397   constructor, add_quality_assessment(), and set_stopping_assessment().
00398   The resulting data is printed in a table unless disable_printing_results()
00399   has been called.  The double returned depends on the QualityMetric
00400   and QAFunction "return" combination, which can be set using
00401   set_stopping_assessemnt().
00402   \param ms (const MeshSet &) MeshSet used for quality assessment.
00403  */
00404 double QualityAssessor::loop_over_mesh( MeshDomainAssoc* mesh_and_domain, const Settings* settings, MsqError& err )
00405 {
00406     return loop_over_mesh_internal( mesh_and_domain, settings, NULL, err );
00407 }
00408 
00409 /*!
00410   Computes the quality data for a given
00411   MeshSet, ms. What quality information is calculated, depends
00412   on what has been requested through the use of the QualityAssessor
00413   constructor, add_quality_assessment(), and set_stopping_assessment().
00414   The resulting data is printed in a table unless disable_printing_results()
00415   has been called.  The double returned depends on the QualityMetric
00416   and QAFunction "return" combination, which can be set using
00417   set_stopping_assessemnt().
00418   \param ms (const MeshSet &) MeshSet used for quality assessment.
00419  */
00420 double QualityAssessor::loop_over_mesh( ParallelMesh* mesh,
00421                                         MeshDomain* domain,
00422                                         const Settings* settings,
00423                                         MsqError& err )
00424 {
00425     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( (Mesh*)mesh, domain, false, true );
00426     return loop_over_mesh_internal( &mesh_and_domain, settings, mesh->get_parallel_helper(), err );
00427 }
00428 
00429 double QualityAssessor::loop_over_mesh_internal( MeshDomainAssoc* mesh_and_domain,
00430                                                  const Settings* settings,
00431                                                  ParallelHelper* helper,
00432                                                  MsqError& err )
00433 {
00434     // Clear out any previous data
00435     reset_data();
00436 
00437     // Clear culling flag, set hard fixed flag, etc on all vertices
00438     initialize_vertex_byte( mesh_and_domain, settings, err );
00439     MSQ_ERRZERO( err );
00440 
00441     Mesh* mesh         = mesh_and_domain->get_mesh();
00442     MeshDomain* domain = mesh_and_domain->get_domain();
00443     invalid_values     = false;
00444 
00445     PatchData patch;
00446     patch.set_mesh( mesh );
00447     patch.set_domain( domain );
00448     if( settings ) patch.attach_settings( settings );
00449 
00450     ElementPatches elem_patches;
00451     elem_patches.set_mesh( mesh );
00452     VertexPatches vert_patches( 1, false );
00453     vert_patches.set_mesh( mesh );
00454 
00455     std::vector< PatchSet::PatchHandle > patches;
00456     std::vector< PatchSet::PatchHandle >::iterator p;
00457     std::vector< Mesh::VertexHandle > patch_verts;
00458     std::vector< Mesh::ElementHandle > patch_elems;
00459     std::vector< size_t > metric_handles;
00460 
00461     // Check if we really need the helper
00462     if( helper && helper->get_nprocs() == 1 ) helper = 0;
00463 
00464     // Get any necessary tag handles.
00465     TagHandle invertedTag = 0;
00466     if( tagging_inverted_elements() )
00467     {
00468         invertedTag = get_tag( mesh, invertedTagName, Mesh::INT, 1, err );
00469         MSQ_ERRZERO( err );
00470     }
00471     list_type::iterator iter;
00472     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00473     {
00474         if( ( *iter )->write_to_tag() )
00475         {
00476             ( *iter )->tagHandle = get_tag( mesh, ( *iter )->tagName, Mesh::DOUBLE, 1, err );
00477             MSQ_ERRZERO( err );
00478         }
00479     }
00480 
00481     TagHandle fixedTag = 0;
00482     if( tagging_fixed_elements() )
00483     {
00484         fixedTag = get_tag( mesh, fixedTagName, Mesh::INT, 1, err );
00485         MSQ_ERRZERO( err );
00486     }
00487 
00488     // Record the type of metric for each assessment so that it can be
00489     // included in the QualitySummary report.
00490     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00491     {
00492         ElementAvgQM* avg_ptr    = dynamic_cast< ElementAvgQM* >( ( *iter )->get_metric() );
00493         ElementMaxQM* max_ptr    = dynamic_cast< ElementMaxQM* >( ( *iter )->get_metric() );
00494         TMPQualityMetric* tq_ptr = dynamic_cast< TMPQualityMetric* >( ( *iter )->get_metric() );
00495         if( avg_ptr )
00496             ( *iter )->assessScheme = ELEMENT_AVG_QM;
00497         else if( max_ptr )
00498             ( *iter )->assessScheme = ELEMENT_MAX_QM;
00499         else if( tq_ptr )
00500             ( *iter )->assessScheme = TMP_QUALITY_METRIC;
00501         else
00502             ( *iter )->assessScheme = QUALITY_METRIC;
00503     }
00504 
00505     // Check for any metrics for which a histogram is to be
00506     // calculated and for which the user has not specified
00507     // minimum and maximum values.
00508     // Element-based metrics are first in list, followed
00509     // by vertex-based metrics.  Find first vertex-based
00510     // metric also such that element metrics go from
00511     // assessList.begin() to elem_end and vertex metrics
00512     // go from elem_end to assessList.end()
00513     list_type::iterator elem_end       = assessList.end();
00514     bool need_second_pass_for_elements = false;
00515     bool need_second_pass_for_vertices = false;
00516     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00517     {
00518         if( ( *iter )->get_metric()->get_metric_type() == QualityMetric::VERTEX_BASED ) break;
00519 
00520         if( ( *iter )->have_histogram() && !( *iter )->haveHistRange ) need_second_pass_for_elements = true;
00521     }
00522     elem_end = iter;
00523     for( ; iter != assessList.end(); ++iter )
00524     {
00525         if( ( *iter )->have_histogram() && !( *iter )->haveHistRange ) need_second_pass_for_vertices = true;
00526     }
00527 
00528     // Do element-based metrics
00529     elem_patches.get_patch_handles( patches, err );
00530     MSQ_ERRZERO( err );
00531 
00532     myData->invertedElementCount = myData->invertedSampleCount = myData->sampleCount = 0;
00533     myData->elementCount                                                             = patches.size();
00534     myData->freeElementCount                                                         = 0;
00535 
00536     int inverted, samples;
00537     bool first_pass = false;
00538     do
00539     {  // might need to loop twice to calculate histograms
00540         first_pass = !first_pass;
00541 
00542         // until there are no more patches
00543         // there is another get_next_patch at
00544         // the end of this loop
00545         for( p = patches.begin(); p != patches.end(); ++p )
00546         {
00547             elem_patches.get_patch( *p, patch_elems, patch_verts, err );
00548             MSQ_ERRZERO( err );
00549 
00550             if( helper )
00551             {
00552                 bool ours = helper->is_our_element( patch_elems[0], err );
00553                 MSQ_ERRZERO( err );
00554                 if( !ours )
00555                 {
00556                     --myData->elementCount;
00557                     continue;
00558                 }
00559             }
00560             patch.set_mesh_entities( patch_elems, patch_verts, err );
00561             MSQ_ERRZERO( err );
00562 
00563             if( first_pass )
00564             {
00565                 // record element type for use in print_summary
00566                 for( size_t i = 0; i < patch.num_elements(); i++ )
00567                 {
00568                     const MsqMeshEntity* elem = &patch.element_by_index( i );
00569                     EntityTopology topo       = elem->get_element_type();
00570                     elementTypeCount[topo - POLYGON]++;
00571                 }
00572             }
00573 
00574             if( 0 == patch.num_free_vertices() )
00575             {
00576                 if( tagging_fixed_elements() )
00577                 {
00578                     Mesh::ElementHandle h = patch.get_element_handles_array()[0];
00579                     int val               = 1;
00580                     mesh->tag_set_element_data( fixedTag, 1, &h, &val, err );
00581                     MSQ_ERRZERO( err );
00582                 }
00583                 if( skipFixedSamples ) continue;
00584             }
00585 
00586             // first check for inverted elements
00587             if( first_pass )
00588             {
00589                 ++myData->freeElementCount;
00590                 inverted = samples = 0;
00591                 patch.element_by_index( 0 ).check_element_orientation( patch, inverted, samples, err );
00592                 MSQ_ERRZERO( err );
00593                 myData->invertedElementCount += ( inverted != 0 );
00594                 myData->invertedSampleCount += inverted;
00595                 myData->sampleCount += samples;
00596 
00597                 if( tagging_inverted_elements() )
00598                 {
00599                     Mesh::ElementHandle h = patch.get_element_handles_array()[0];
00600                     int val               = ( inverted != 0 );
00601                     mesh->tag_set_element_data( invertedTag, 1, &h, &val, err );
00602                     MSQ_ERRZERO( err );
00603                 }
00604             }
00605 
00606             // now process all element-based metrics
00607             for( iter = assessList.begin(); iter != elem_end; ++iter )
00608             {
00609                 // If first pass, get values for all metrics
00610                 if( first_pass )
00611                 {
00612                     double value = 0.0;
00613                     metric_handles.clear();
00614                     QualityMetric* qm = ( *iter )->get_metric();
00615                     qm->get_single_pass( patch, metric_handles, skipFixedSamples, err );
00616                     MSQ_ERRZERO( err );
00617                     for( std::vector< size_t >::iterator j = metric_handles.begin(); j != metric_handles.end(); ++j )
00618                     {
00619                         bool valid = ( *iter )->get_metric()->evaluate( patch, *j, value, err );  // MSQ_ERRZERO(err);
00620                         if( err.error_code() == err.BARRIER_VIOLATED )
00621                         {
00622                             err.clear();
00623                             invalid_values = true;
00624                         }
00625                         ( *iter )->add_value( value );
00626                         if( !valid ) ( *iter )->add_invalid_value();
00627                     }
00628                     // we don't do tag stuff unless metric is truely element-based
00629                     // (only one value per element)
00630                     if( ( *iter )->write_to_tag() && metric_handles.size() == 1 )
00631                     {
00632                         Mesh::ElementHandle h = patch.get_element_handles_array()[0];
00633                         mesh->tag_set_element_data( ( *iter )->tagHandle, 1, &h, &value, err );
00634                         MSQ_ERRZERO( err );
00635                     }
00636                 }
00637                 // If second pass, only do metrics for which the
00638                 // histogram hasn't been calculated yet.
00639                 else if( ( *iter )->have_histogram() && !( *iter )->haveHistRange )
00640                 {
00641                     metric_handles.clear();
00642                     QualityMetric* qm = ( *iter )->get_metric();
00643                     qm->get_evaluations( patch, metric_handles, skipFixedSamples, err );
00644                     MSQ_ERRZERO( err );
00645                     for( std::vector< size_t >::iterator j = metric_handles.begin(); j != metric_handles.end(); ++j )
00646                     {
00647                         double value;
00648                         ( *iter )->get_metric()->evaluate( patch, *j, value, err );
00649                         MSQ_ERRZERO( err );
00650                         ( *iter )->add_hist_value( value );
00651                     }
00652                 }
00653             }
00654         }
00655         if( MSQ_CHKERR( err ) ) return 0.0;
00656 
00657         // Fix up any histogram ranges which were calculated
00658         for( iter = assessList.begin(); iter != elem_end; ++iter )
00659             if( ( *iter )->have_histogram() && !( *iter )->haveHistRange )
00660                 if( first_pass )
00661                 {
00662                     if( helper )
00663                     {
00664                         helper->communicate_min_max_to_all( &( ( *iter )->minimum ), &( ( *iter )->maximum ), err );
00665                         MSQ_ERRZERO( err );
00666                     }
00667 
00668                     ( *iter )->calculate_histogram_range();
00669                     // Uncomment the following to have the QA keep the first
00670                     // calculated histogram range for all subsequent iterations.
00671                     //          else
00672                     //            (*iter)->haveHistRange = true;
00673                 }
00674     } while( first_pass && need_second_pass_for_elements && !invalid_values );
00675 
00676     // Do vertex-based metrics
00677     if( assessList.end() != elem_end )
00678     {
00679         vert_patches.get_patch_handles( patches, err );
00680         MSQ_ERRZERO( err );
00681 
00682         first_pass = false;
00683         do
00684         {  // might need to loop twice to calculate histograms
00685             first_pass = !first_pass;
00686 
00687             // until there are no more patches
00688             // there is another get_next_patch at
00689             // the end of this loop
00690             for( p = patches.begin(); p != patches.end(); ++p )
00691             {
00692                 vert_patches.get_patch( *p, patch_elems, patch_verts, err );
00693                 MSQ_ERRZERO( err );
00694                 patch.set_mesh_entities( patch_elems, patch_verts, err );
00695                 MSQ_ERRZERO( err );
00696                 if( skipFixedSamples && 0 == patch.num_free_vertices() ) continue;
00697 
00698                 Mesh::VertexHandle vert_handle = reinterpret_cast< Mesh::VertexHandle >( *p );
00699 
00700                 for( iter = elem_end; iter != assessList.end(); ++iter )
00701                 {
00702                     // If first pass, get values for all metrics
00703                     if( first_pass )
00704                     {
00705                         double value = 0.0;
00706                         metric_handles.clear();
00707                         QualityMetric* qm = ( *iter )->get_metric();
00708                         qm->get_single_pass( patch, metric_handles, skipFixedSamples, err );
00709                         MSQ_ERRZERO( err );
00710                         for( std::vector< size_t >::iterator j = metric_handles.begin(); j != metric_handles.end();
00711                              ++j )
00712                         {
00713                             bool valid = ( *iter )->get_metric()->evaluate( patch, *j, value, err );
00714                             MSQ_ERRZERO( err );
00715                             ( *iter )->add_value( value );
00716                             if( !valid ) ( *iter )->add_invalid_value();
00717                         }
00718                         // we don't do tag stuff unless metric is truely vertex-based
00719                         // (only one value per vertex)
00720                         if( ( *iter )->write_to_tag() && metric_handles.size() == 1 )
00721                         {
00722                             mesh->tag_set_vertex_data( ( *iter )->tagHandle, 1, &vert_handle, &value, err );
00723                             MSQ_ERRZERO( err );
00724                         }
00725                     }
00726                     // If second pass, only do metrics for which the
00727                     // histogram hasn't been calculated yet.
00728                     else if( ( *iter )->have_histogram() && !( *iter )->haveHistRange )
00729                     {
00730                         metric_handles.clear();
00731                         QualityMetric* qm = ( *iter )->get_metric();
00732                         qm->get_evaluations( patch, metric_handles, skipFixedSamples, err );
00733                         MSQ_ERRZERO( err );
00734                         for( std::vector< size_t >::iterator j = metric_handles.begin(); j != metric_handles.end();
00735                              ++j )
00736                         {
00737                             double value;
00738                             ( *iter )->get_metric()->evaluate( patch, *j, value, err );
00739                             MSQ_ERRZERO( err );
00740                             ( *iter )->add_hist_value( value );
00741                         }
00742                     }
00743                 }
00744             }
00745             if( MSQ_CHKERR( err ) ) return 0.0;
00746 
00747             // Fix up any histogram ranges which were calculated
00748             for( iter = elem_end; iter != assessList.end(); ++iter )
00749                 if( ( *iter )->have_histogram() && !( *iter )->haveHistRange )
00750                     if( first_pass )
00751                     {
00752                         if( helper )
00753                         {
00754                             helper->communicate_min_max_to_all( &( ( *iter )->minimum ), &( ( *iter )->maximum ), err );
00755                             MSQ_ERRZERO( err );
00756                         }
00757                         ( *iter )->calculate_histogram_range();
00758                         // Uncomment the following to have the QA keep the first
00759                         // calculated histogram range for all subsequent iterations.
00760                         //          else
00761                         //            (*iter)->haveHistRange = true;
00762                     }
00763         } while( first_pass && need_second_pass_for_vertices && !invalid_values );
00764     }
00765 
00766     if( helper )
00767     {
00768         for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00769         {
00770 
00771             helper->communicate_min_max_to_zero( &( ( *iter )->minimum ), &( ( *iter )->maximum ), err );
00772             MSQ_ERRZERO( err );
00773 
00774             helper->communicate_sums_to_zero( &myData->freeElementCount, &myData->invertedElementCount,
00775                                               &myData->elementCount, &myData->invertedSampleCount, &myData->sampleCount,
00776                                               &( ( *iter )->count ), &( ( *iter )->numInvalid ), &( ( *iter )->sum ),
00777                                               &( ( *iter )->sqrSum ), err );
00778             MSQ_ERRZERO( err );
00779 
00780             if( ( *iter )->have_power_mean() )
00781             {
00782                 helper->communicate_power_sum_to_zero( &( ( *iter )->pMean ), err );
00783                 MSQ_ERRZERO( err );
00784             }
00785 
00786             if( ( *iter )->have_histogram() )
00787             {
00788                 helper->communicate_histogram_to_zero( ( *iter )->histogram, err );
00789                 MSQ_ERRZERO( err );
00790             }
00791         }
00792     }
00793 
00794     // Print results, if requested
00795     if( printSummary && ( !helper || helper->get_rank() == 0 ) ) print_summary( this->outputStream );
00796 
00797     list_type::iterator i = find_stopping_assessment();
00798     return i == assessList.end() ? 0 : ( *i )->stopping_function_value();
00799 }
00800 
00801 bool QualityAssessor::invalid_elements() const
00802 {
00803     bool result = false;
00804     list_type::const_iterator iter;
00805     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00806         if( ( *iter )->get_invalid_element_count() ) result = true;
00807     return result;
00808 }
00809 
00810 void QualityAssessor::reset_data()
00811 {
00812     list_type::iterator iter;
00813     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00814         ( *iter )->reset_data();
00815     myData->invertedElementCount = -1;
00816     myData->invertedSampleCount  = -1;
00817 
00818     for( int i = POLYGON; i <= MIXED; i++ )
00819         elementTypeCount[i - POLYGON] = 0;
00820 }
00821 
00822 void QualityAssessor::scale_histograms( QualityAssessor* optimized )
00823 {
00824     // find the histograms to scale
00825     list_type::iterator iter;
00826     list_type::iterator initial;
00827     list_type::iterator optimal;
00828     bool initial_found = false, optimal_found = false;
00829     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
00830     {
00831         if( ( *iter )->histogram.size() > 0 )
00832         {
00833             if( initial_found == false )
00834             {
00835                 initial       = iter;
00836                 initial_found = true;
00837                 break;
00838             }
00839         }
00840     }
00841     for( iter = optimized->assessList.begin(); iter != optimized->assessList.end(); ++iter )
00842     {
00843         if( ( *iter )->histogram.size() > 0 )
00844         {
00845             optimal       = iter;
00846             optimal_found = true;
00847             break;
00848         }
00849     }
00850     if( !initial_found || !optimal_found )
00851     {
00852         // issue warning: orig histograms not found
00853         if( !initial_found )
00854             outputStream << "WARNING: 'before' histogram not found" << std::endl;
00855         else
00856             outputStream << "WARNING: 'after' histogram not found" << std::endl;
00857         return;
00858     }
00859 
00860     // check number of intervals (bins) for each histogram
00861     int num_intervals = ( *initial )->histogram.size() - 2;
00862     if( num_intervals != int( ( *optimal )->histogram.size() - 2 ) )
00863     {
00864         // issue warning: number of intervals not the same
00865         outputStream << "WARNING: histogram intervals are not the same" << std::endl;
00866         return;
00867     }
00868 
00869     // calculate new max and min values
00870     double combined_min, combined_max;
00871     if( ( *initial )->histMin < ( *optimal )->histMin )
00872         combined_min = ( *initial )->histMin;
00873     else
00874         combined_min = ( *optimal )->histMin;
00875 
00876     if( ( *initial )->histMax > ( *optimal )->histMax )
00877         combined_max = ( *initial )->histMax;
00878     else
00879         combined_max = ( *optimal )->histMax;
00880 
00881     // put the before quality values into the correct new bins
00882 
00883     // First and last values in array are counts of valuesnum_intervals+1
00884     // outside the user-specified range of the histogram
00885     // (below and above, respectively.)
00886     std::vector< int > new_initial_histogram;
00887     new_initial_histogram.resize( ( *initial )->histogram.size(), 0 );
00888     new_initial_histogram[0] = ( *initial )->histogram[0];
00889     new_initial_histogram[new_initial_histogram.size() - 1] =
00890         ( *initial )->histogram[( *initial )->histogram.size() - 1];
00891 
00892     // Re-calculate which interval the value is in.  Add one
00893     // because first entry is for values below user-specifed
00894     // minimum value for histogram.
00895     double combined_range = combined_max - combined_min;
00896     double initial_min    = ( *initial )->histMin;
00897     double optimal_min    = ( *optimal )->histMin;
00898     double initial_range  = ( *initial )->histMax - ( *initial )->histMin;
00899     double optimal_range  = ( *optimal )->histMax - ( *optimal )->histMin;
00900     double combined_step  = combined_range / num_intervals;
00901     double initial_step   = initial_range / num_intervals;
00902     double optimal_step   = optimal_range / num_intervals;
00903     // round steps to 3 significant digits
00904     if( combined_step >= 0.001 ) combined_step = round_to_3_significant_digits( combined_step );
00905     if( initial_step >= 0.001 ) initial_step = round_to_3_significant_digits( initial_step );
00906     if( optimal_step >= 0.001 ) optimal_step = round_to_3_significant_digits( optimal_step );
00907 
00908     // populate initial histogram
00909     if( initial_range == combined_range )
00910     {
00911         // just copy histogram
00912         new_initial_histogram = ( *initial )->histogram;
00913     }
00914     else
00915     {
00916         for( size_t i = 1; i < new_initial_histogram.size() - 1; i++ )
00917         {
00918             double combined_bin_value = combined_min + ( combined_step * ( i - 1 ) );
00919             for( size_t j = 1; j < new_initial_histogram.size() - 1; j++ )
00920             {
00921                 double initial_bin_value = initial_min + ( initial_step * ( j - 1 ) );
00922                 if( initial_bin_value >= combined_bin_value && initial_bin_value < combined_bin_value + combined_step )
00923                 {
00924                     new_initial_histogram[i] += ( *initial )->histogram[j];
00925                 }
00926             }
00927         }
00928     }
00929 
00930     // put the optimal quality values into the correct new bins
00931     std::vector< int > new_optimal_histogram;
00932     new_optimal_histogram.resize( ( *optimal )->histogram.size(), 0 );
00933     new_optimal_histogram[0] = ( *optimal )->histogram[0];
00934     new_optimal_histogram[new_optimal_histogram.size() - 1] =
00935         ( *optimal )->histogram[( *optimal )->histogram.size() - 1];
00936 
00937     // populate optimal histogram
00938     if( optimal_range == combined_range )
00939     {
00940         // just copy histogram
00941         new_optimal_histogram = ( *optimal )->histogram;
00942     }
00943     else
00944     {
00945         for( size_t i = 1; i < new_optimal_histogram.size() - 1; i++ )
00946         {
00947             double combined_bin_value = combined_min + ( combined_step * ( i - 1 ) );
00948             for( size_t j = 1; j < new_optimal_histogram.size() - 1; j++ )
00949             {
00950                 double optimal_bin_value = optimal_min + ( optimal_step * ( j - 1 ) );
00951                 if( optimal_bin_value >= combined_bin_value && optimal_bin_value < combined_bin_value + combined_step )
00952                 {
00953                     new_optimal_histogram[i] += ( *optimal )->histogram[j];
00954                 }
00955             }
00956         }
00957     }
00958 
00959     // determine largest number of values in a 'bin' for both histograms
00960     unsigned i;
00961     int max_interval_num = 1;
00962     for( i = 0; i < new_initial_histogram.size(); ++i )
00963     {
00964         if( new_initial_histogram[i] > max_interval_num ) max_interval_num = new_initial_histogram[i];
00965     }
00966     for( i = 0; i < new_optimal_histogram.size(); ++i )
00967     {
00968         if( new_optimal_histogram[i] > max_interval_num ) max_interval_num = new_optimal_histogram[i];
00969     }
00970 
00971     // calculate how many bar graph characters will represent the
00972     // largest 'bin' value.
00973     // create the 'before' histogram
00974     int termwidth         = get_terminal_width();
00975     const char indent[]   = "   ";
00976     const char GRAPH_CHAR = '=';                              // Character used to create bar graphs
00977     const int TOTAL_WIDTH = termwidth > 30 ? termwidth : 70;  // Width of histogram
00978     int GRAPHW            = TOTAL_WIDTH - sizeof( indent );
00979 
00980     if( 0 == max_interval_num ) return;  // no data
00981 
00982     // Calculate width of field containing counts for
00983     // histogram intervals (log10(max_interval)).
00984     int num_width = 1;
00985     for( int temp = max_interval_num; temp > 0; temp /= 10 )
00986         ++num_width;
00987     GRAPHW -= num_width;
00988 
00989     // Create an array of bar graph characters for use in output
00990     std::vector< char > graph_chars( GRAPHW + 1, GRAPH_CHAR );
00991 
00992     // Check if bar-graph should be linear or log10 plot
00993     // Do log plot if standard deviation is less that 1.5
00994     // histogram intervals.
00995 
00996     bool log_plot = false;
00997     double stddev = ( *initial )->get_stddev();
00998     if( stddev > 0 && stddev < 2.0 * combined_step )
00999     {
01000         int new_interval = (int)( log10( (double)( 1 + max_interval_num ) ) );
01001         if( new_interval > 0 )
01002         {
01003             log_plot         = true;
01004             max_interval_num = new_interval;
01005         }
01006     }
01007     // output the 'before' histogram
01008 
01009     // Write title
01010     outputStream << std::endl << "************** Common-scale Histograms **************" << std::endl << std::endl;
01011     outputStream << indent << ( *initial )->get_label() << " histogram (initial mesh):";
01012     if( log_plot ) outputStream << " (log10 plot)";
01013     outputStream << std::endl;
01014 
01015     // Calculate width of a single quality interval value
01016     double interval_value  = 0.0;
01017     int max_interval_width = 0;
01018     std::stringstream str_stream;
01019     std::string interval_string;
01020     for( i = 0; i < new_initial_histogram.size(); ++i )
01021     {
01022         interval_value = combined_min + (i)*combined_step;
01023         if( combined_step >= .001 ) interval_value = round_to_3_significant_digits( interval_value );
01024         str_stream.clear();
01025         str_stream.str( "" );
01026         interval_string = "";
01027         str_stream << interval_value;
01028         interval_string = str_stream.str();
01029         if( interval_string.length() > (size_t)max_interval_width ) max_interval_width = interval_string.length();
01030     }
01031 
01032     // adjust graph width for actual size of interval values
01033     GRAPHW = GRAPHW - ( max_interval_width * 2 ) - 5;
01034 
01035     // For each interval of histogram
01036     for( i = 0; i < new_initial_histogram.size(); ++i )
01037     {
01038         // First value is the count of the number of values that
01039         // were below the minimum value of the histogram.
01040         if( 0 == i )
01041         {
01042             if( 0 == new_initial_histogram[i] ) continue;
01043             outputStream << indent << std::setw( max_interval_width ) << "under min";
01044         }
01045         // Last value is the count of the number of values that
01046         // were above the maximum value of the histogram.
01047         else if( i + 1 == new_initial_histogram.size() )
01048         {
01049             if( 0 == new_initial_histogram[i] ) continue;
01050             outputStream << indent << std::setw( max_interval_width ) << "over max";
01051         }
01052         // Anything else is a valid interval of the histogram.
01053         // Print the range for each interval.
01054         else
01055         {
01056             double start_value = combined_min + ( i - 1 ) * combined_step;
01057             double end_value   = combined_min + (i)*combined_step;
01058 
01059             if( combined_step >= 0.001 )
01060             {
01061                 start_value = round_to_3_significant_digits( start_value );
01062                 end_value   = round_to_3_significant_digits( end_value );
01063             }
01064 
01065             outputStream << indent << "(" << std::setw( max_interval_width ) << std::right << start_value << "-"
01066                          << std::setw( max_interval_width ) << std::left << end_value << ") |";
01067         }
01068 
01069         // Print bar graph
01070 
01071         // First calculate the number of characters to output
01072         int num_graph;
01073         if( log_plot )
01074             num_graph = GRAPHW * (int)log10( (double)( 1 + new_initial_histogram[i] ) ) / max_interval_num;
01075         else
01076             num_graph = GRAPHW * new_initial_histogram[i] / max_interval_num;
01077 
01078         // print num_graph characters using array of fill characters.
01079         graph_chars[num_graph] = '\0';
01080         outputStream << arrptr( graph_chars );
01081         graph_chars[num_graph] = GRAPH_CHAR;
01082 
01083         // Print interval count.
01084         outputStream << new_initial_histogram[i] << std::endl;
01085     }
01086 
01087     outputStream << "  metric was evaluated " << ( *initial )->count << " times." << std::endl << std::endl;
01088 
01089     // output the 'after' histogram
01090     outputStream << std::endl << indent << ( *optimal )->get_label() << " histogram (optimal mesh):";
01091     if( log_plot ) outputStream << " (log10 plot)";
01092     outputStream << std::endl;
01093 
01094     // For each interval of histogram
01095     for( i = 0; i < new_optimal_histogram.size(); ++i )
01096     {
01097         // First value is the count of the number of values that
01098         // were below the minimum value of the histogram.
01099         if( 0 == i )
01100         {
01101             if( 0 == new_optimal_histogram[i] ) continue;
01102             outputStream << indent << std::setw( max_interval_width ) << "under min";
01103         }
01104         // Last value is the count of the number of values that
01105         // were above the maximum value of the histogram.
01106         else if( i + 1 == new_optimal_histogram.size() )
01107         {
01108             if( 0 == new_optimal_histogram[i] ) continue;
01109             outputStream << indent << std::setw( max_interval_width ) << "over max";
01110         }
01111         // Anything else is a valid interval of the histogram.
01112         // Print the range for each interval.
01113         else
01114         {
01115             double start_value = combined_min + ( i - 1 ) * combined_step;
01116             double end_value   = combined_min + (i)*combined_step;
01117 
01118             if( combined_step >= 0.001 )
01119             {
01120                 start_value = round_to_3_significant_digits( start_value );
01121                 end_value   = round_to_3_significant_digits( end_value );
01122             }
01123 
01124             outputStream << indent << "(" << std::setw( max_interval_width ) << std::right << start_value << "-"
01125                          << std::setw( max_interval_width ) << std::left << end_value << ") |";
01126         }
01127 
01128         // Print bar graph
01129 
01130         // First calculate the number of characters to output
01131         int num_graph;
01132         if( log_plot )
01133             num_graph = GRAPHW * (int)log10( (double)( 1 + new_optimal_histogram[i] ) ) / max_interval_num;
01134         else
01135             num_graph = GRAPHW * new_optimal_histogram[i] / max_interval_num;
01136 
01137         // print num_graph characters using array of fill characters.
01138         graph_chars[num_graph] = '\0';
01139         outputStream << arrptr( graph_chars );
01140         graph_chars[num_graph] = GRAPH_CHAR;
01141 
01142         // Print interval count.
01143         outputStream << new_optimal_histogram[i] << std::endl;
01144     }
01145 
01146     outputStream << "  metric was evaluated " << ( *optimal )->count << " times." << std::endl << std::endl;
01147 
01148     return;
01149 }
01150 
01151 QualityAssessor::Assessor::Assessor( QualityMetric* metric, const char* label )
01152     : qualMetric( metric ), mLabel( label ? std::string( label ) : metric->get_name() ), pMean( 0.0 ),
01153       haveHistRange( false ), histMin( 1.0 ), histMax( 0.0 ), tagHandle( 0 ), stoppingFunction( false ),
01154       referenceCount( 0 ), assessScheme( NO_SCHEME )
01155 {
01156     reset_data();
01157 }
01158 
01159 const QualityAssessor::Assessor* QualityAssessor::get_results( QualityMetric* metric ) const
01160 {
01161     list_type::const_iterator iter;
01162     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01163         if( ( *iter )->get_metric() == metric ) return *iter;
01164     return 0;
01165 }
01166 
01167 const QualityAssessor::Assessor* QualityAssessor::get_results( const char* name ) const
01168 {
01169     list_type::const_iterator iter;
01170     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01171         if( ( *iter )->get_label() == name ) return *iter;
01172     return 0;
01173 }
01174 
01175 void QualityAssessor::Assessor::get_histogram( double& lower_bound_out,
01176                                                double& upper_bound_out,
01177                                                std::vector< int >& counts_out,
01178                                                MsqError& err ) const
01179 {
01180     if( !have_histogram() )
01181     {
01182         MSQ_SETERR( err )( "No histogram calculated.", MsqError::INVALID_STATE );
01183         return;
01184     }
01185 
01186     lower_bound_out = histMin;
01187     upper_bound_out = histMax;
01188     counts_out      = histogram;
01189 }
01190 
01191 void QualityAssessor::Assessor::reset_data()
01192 {
01193     count        = 0;
01194     sum          = 0;
01195     maximum      = -HUGE_VAL;
01196     minimum      = HUGE_VAL;
01197     sqrSum       = 0;
01198     pSum         = 0;
01199     numInvalid   = 0;
01200     assessScheme = NO_SCHEME;
01201     // zero histogram data
01202     size_t hist_size = histogram.size();
01203     histogram.clear();
01204     histogram.resize( hist_size, 0 );
01205 }
01206 
01207 void QualityAssessor::Assessor::add_value( double metric_value )
01208 {
01209     sum += metric_value;
01210     sqrSum += metric_value * metric_value;
01211     if( metric_value > maximum ) maximum = metric_value;
01212     if( metric_value < minimum ) minimum = metric_value;
01213     // Only add value to histogram data from this function if
01214     // the user has specified the range.  If user has not
01215     // specified the range, QualityAssessor will call add_hist_value()
01216     // directly once the range has been calculated.
01217     if( have_histogram() && haveHistRange ) add_hist_value( metric_value );
01218 
01219     if( have_power_mean() ) pSum += pow( metric_value, pMean );
01220 
01221     ++count;
01222 }
01223 
01224 void QualityAssessor::Assessor::add_invalid_value()
01225 {
01226     ++numInvalid;
01227 }
01228 
01229 void QualityAssessor::Assessor::add_hist_value( double metric_value )
01230 {
01231     // First and last values in array are counts of values
01232     // outside the user-specified range of the histogram
01233     // (below and above, respectively.)
01234     if( metric_value < histMin )
01235         ++histogram[0];
01236     else if( metric_value > histMax )
01237         ++histogram[histogram.size() - 1];
01238     else
01239     {
01240         // Calculate which interval the value is in.  Add one
01241         // because first entry is for values below user-specifed
01242         // minimum value for histogram.
01243         double range = histMax - histMin;
01244         double fract;
01245         if( range > DBL_EPSILON )
01246             fract = ( metric_value - histMin ) / range;
01247         else
01248             fract = 0.0;
01249         unsigned cell;
01250         if( fabs( fract - 1.0 ) < histMax * DBL_EPSILON )
01251             cell = histogram.size() - 1;
01252         else
01253             cell = 1 + (int)( ( fract * ( histogram.size() - 2 ) ) );
01254 
01255         // Add value to interval.
01256         ++histogram[cell];
01257     }
01258 }
01259 
01260 void QualityAssessor::Assessor::calculate_histogram_range()
01261 {
01262     double lower_bound = minimum;
01263     int num_intervals  = histogram.size();
01264     double step        = ( maximum - lower_bound ) / num_intervals;
01265     if( step == 0 ) step = 1.0;
01266     double size = pow( 10.0, floor( log10( step / ( num_intervals - 1 ) ) ) );
01267     if( size < 1e-6 ) size = 1.0;
01268     histMin = lower_bound;
01269     histMax = lower_bound + num_intervals * size * ceil( step / size );
01270 }
01271 
01272 void QualityAssessor::print_summary( std::ostream& stream ) const
01273 {
01274     const int NAMEW = 19;  // Width of name column in table output
01275     const int NUMW  = 12;  // Width of value columns in table output
01276 
01277     // Print title
01278     stream << std::endl
01279            << "************** " << qualityAssessorName << " Summary **************" << std::endl
01280            << std::endl;
01281 
01282     if( myData->freeElementCount != myData->elementCount )
01283         stream << "  Evaluating quality for " << myData->freeElementCount << " free elements of "
01284                << myData->elementCount << " total elements." << std::endl;
01285     else
01286         stream << "  Evaluating quality for " << myData->elementCount << " elements." << std::endl;
01287 
01288     // Print element type totals
01289     std::string str_value = "";
01290     for( int i = POLYGON; i <= MIXED; i++ )
01291     {
01292         if( elementTypeCount[i - POLYGON] )
01293         {
01294             str_value = element_name_as_string( i );
01295             if( str_value != "" )
01296                 stream << "  This mesh had " << elementTypeCount[i - POLYGON] << " " << str_value << " elements."
01297                        << std::endl;
01298         }
01299     }
01300 
01301     if( myData->invertedElementCount )
01302     {
01303         stream << "  THERE ARE " << myData->invertedElementCount << " INVERTED ELEMENTS. " << std::endl
01304                << "  (Elements invalid at " << myData->invertedSampleCount << " of " << myData->sampleCount
01305                << " sample locations.)" << std::endl
01306                << std::endl;
01307     }
01308     else
01309     {
01310         stream << "  There were no inverted elements detected. " << std::endl;
01311     }
01312 
01313     // Check if there are invalid values for any metric
01314     list_type::const_iterator iter;
01315     bool some_invalid = false;
01316     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01317     {
01318         if( ( *iter )->get_invalid_element_count() )
01319         {
01320             some_invalid = true;
01321             stream << "  " << ( *iter )->get_invalid_element_count() << " OF " << ( *iter )->get_count()
01322                    << " ENTITIES EVALUATED TO AN UNDEFINED VALUE FOR " << ( *iter )->get_label() << std::endl
01323                    << std::endl;
01324         }
01325     }
01326     if( !some_invalid )
01327     {
01328         stream << "  No entities had undefined values for any computed metric." << std::endl << std::endl;
01329     }
01330 
01331     if( !invalid_values )
01332     {
01333         // Check if a user-define power-mean was calculated for any of the metrics
01334         std::set< double > pmeans;
01335         for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01336             if( ( *iter )->have_power_mean() ) pmeans.insert( ( *iter )->get_power() );
01337 
01338         // If power-mean of 1 or 2 was requested, change titles rather
01339         // than printing redundant values
01340         std::string average_str( "average" ), rms_str( "rms" );
01341         if( pmeans.find( 1.0 ) != pmeans.end() )
01342         {
01343             pmeans.erase( pmeans.find( 1.0 ) );
01344             average_str = "1-mean";
01345         }
01346         if( pmeans.find( 2.0 ) != pmeans.end() )
01347         {
01348             pmeans.erase( pmeans.find( 2.0 ) );
01349             rms_str = "2-mean";
01350         }
01351 
01352         // Number of values in table
01353         unsigned num_values = pmeans.size() + 5;
01354 
01355         // Decide how wide of a table field should be used for the metric name
01356         unsigned twidth    = get_terminal_width();
01357         unsigned maxnwidth = NAMEW;
01358         if( twidth )
01359         {
01360             unsigned valwidth = NUMW * num_values;
01361             maxnwidth         = valwidth < twidth ? twidth - valwidth : 0;
01362         }
01363         unsigned namewidth = 0;
01364         for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01365             if( ( *iter )->get_label().size() > namewidth ) namewidth = ( *iter )->get_label().size();
01366         if( namewidth > maxnwidth ) namewidth = maxnwidth;
01367         if( namewidth < 7 )  // at least enough width for the column header
01368             namewidth = 7;
01369 
01370         int number_of_assessments = 0;
01371 
01372         // print metric values
01373         for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01374         {
01375             if( number_of_assessments > 0 ) stream << "    -------------------------------------------" << std::endl;
01376 
01377             // print assessment method used to calculate the statistics
01378             if( ( *iter )->assessScheme == TMP_QUALITY_METRIC )
01379                 stream << "     Sample Point Quality Statistics" << std::endl << std::endl;
01380             else
01381                 stream << "     Element Quality Statistics" << std::endl << std::endl;
01382 
01383             // print comlumn label line
01384             std::set< double >::const_iterator piter;
01385             stream << std::setw( NUMW ) << "minimum";
01386             for( piter = pmeans.begin(); piter != pmeans.end() && *piter < 1.0; ++piter )
01387                 stream << std::setw( NUMW - 6 ) << *piter << "-mean ";
01388             stream << std::setw( NUMW ) << average_str;
01389             for( ; piter != pmeans.end() && *piter < 2.0; ++piter )
01390                 stream << std::setw( NUMW - 6 ) << *piter << "-mean ";
01391             stream << std::setw( NUMW ) << rms_str;
01392             for( ; piter != pmeans.end(); ++piter )
01393                 stream << std::setw( NUMW - 6 ) << *piter << "-mean ";
01394             stream << std::setw( NUMW ) << "maximum";
01395             stream << std::setw( NUMW ) << "std.dev.";
01396             stream << std::endl;
01397 
01398             stream << std::setw( NUMW ) << ( *iter )->get_minimum();
01399             // print power-means with P less than 1.0
01400             for( piter = pmeans.begin(); piter != pmeans.end() && *piter < 1.0; ++piter )
01401             {
01402                 if( ( *iter )->have_power_mean() && ( *iter )->get_power() == *piter )
01403                     stream << std::setw( NUMW ) << ( *iter )->get_power_mean();
01404                 else
01405                     stream << std::setw( NUMW ) << " ";
01406             }
01407             // print average
01408             stream << std::setw( NUMW ) << ( *iter )->get_average();
01409             // print power-means with P less than 2.0
01410             for( ; piter != pmeans.end() && *piter < 2.0; ++piter )
01411             {
01412                 if( ( *iter )->have_power_mean() && ( *iter )->get_power() == *piter )
01413                     stream << std::setw( NUMW ) << ( *iter )->get_power_mean();
01414                 else
01415                     stream << std::setw( NUMW ) << " ";
01416             }
01417             // print RMS
01418             stream << std::setw( NUMW ) << ( *iter )->get_rms();
01419             // print power-means with P greater than 2.0
01420             for( ; piter != pmeans.end(); ++piter )
01421             {
01422                 if( ( *iter )->have_power_mean() && ( *iter )->get_power() == *piter )
01423                     stream << std::setw( NUMW ) << ( *iter )->get_power_mean();
01424                 else
01425                     stream << std::setw( NUMW ) << " ";
01426             }
01427             // print maximum and standard deviation
01428             stream << std::setw( NUMW ) << ( *iter )->get_maximum();
01429             stream << std::setw( NUMW ) << ( *iter )->get_stddev();
01430             stream << std::endl << std::endl;
01431 
01432             stream << "     Number of statistics = " << ( *iter )->get_count() << std::endl;
01433 
01434             // print name
01435             stream << "     Metric = " << ( *iter )->get_label() << std::endl;
01436 
01437             // Output the method used to calcualte the quality values
01438             switch( ( *iter )->assessScheme )
01439             {
01440                 case ELEMENT_AVG_QM:
01441                     stream << "     Element Quality = average over metric values at the elements' "
01442                               "sample points"
01443                            << std::endl
01444                            << std::endl;
01445                     break;
01446                 case ELEMENT_MAX_QM:
01447                     stream << "     Element Quality = maximum over metric values at the elements' "
01448                               "sample points"
01449                            << std::endl
01450                            << std::endl;
01451                     break;
01452                 case TMP_QUALITY_METRIC:
01453                     stream << std::endl << std::endl;
01454                     break;
01455                 case QUALITY_METRIC:
01456                     stream << "     Element Quality not based on sample points." << std::endl << std::endl;
01457                     break;
01458                 default:
01459                     stream << "     Scheme used for deriving qualitiy values unknown" << std::endl << std::endl;
01460             }
01461             if( ( *iter )->have_histogram() ) ( *iter )->print_histogram( stream, get_terminal_width() );
01462             number_of_assessments++;
01463         }
01464     }
01465     else
01466     {
01467         stream << "  Element Quality Statistics are Undefined for this Metric because" << std::endl
01468                << "  there are Inverted Elements" << std::endl
01469                << std::endl;
01470     }
01471 }
01472 
01473 void QualityAssessor::Assessor::print_histogram( std::ostream& stream, int termwidth ) const
01474 {
01475     // Portability notes:
01476     //  Use log10 rather than log10f because the float variations require
01477     //  including platform-dependent headers on some platforms.
01478     //  Explicitly cast log10 argument to double because some platforms
01479     //  have overloaded float and double variations in C++ making an
01480     //  implicit cast from an integer ambiguous.
01481 
01482     const char indent[]   = "   ";
01483     const char GRAPH_CHAR = '=';                              // Character used to create bar graphs
01484     const int TOTAL_WIDTH = termwidth > 30 ? termwidth : 70;  // Width of histogram
01485     int GRAPHW            = TOTAL_WIDTH - sizeof( indent );
01486 
01487     // range is either user-specified (histMin & histMax) or
01488     // calculated (minimum & maximum)
01489     double min, max;
01490     min = histMin;
01491     max = histMax;
01492     // Witdh of one interval of histogram
01493     double step = ( max - min ) / ( histogram.size() - 2 );
01494     // round step to 3 significant digits
01495 
01496     if( step >= 0.001 ) step = round_to_3_significant_digits( step );
01497 
01498     // Find maximum value for an interval bin of the histogram
01499     unsigned i;
01500     int max_bin_value = 1;
01501     for( i = 0; i < histogram.size(); ++i )
01502         if( histogram[i] > max_bin_value ) max_bin_value = histogram[i];
01503 
01504     if( 0 == max_bin_value ) return;  // no data
01505 
01506     // Calculate width of field containing counts for
01507     // histogram intervals (log10(max_bin_value)).
01508     int num_width = 1;
01509     for( int temp = max_bin_value; temp > 0; temp /= 10 )
01510         ++num_width;
01511     GRAPHW -= num_width;
01512 
01513     // Create an array of bar graph characters for use in output
01514     std::vector< char > graph_chars( GRAPHW + 1, GRAPH_CHAR );
01515 
01516     // Check if bar-graph should be linear or log10 plot
01517     // Do log plot if standard deviation is less that 1.5
01518     // histogram intervals.
01519     bool log_plot = false;
01520     double stddev = get_stddev();
01521     if( stddev > 0 && stddev < 2.0 * step )
01522     {
01523         int new_interval = (int)( log10( (double)( 1 + max_bin_value ) ) );
01524         if( new_interval > 0 )
01525         {
01526             log_plot      = true;
01527             max_bin_value = new_interval;
01528         }
01529     }
01530 
01531     // Write title
01532     stream << indent << get_label() << " histogram:";
01533     if( log_plot ) stream << " (log10 plot)";
01534     stream << std::endl;
01535 
01536     // Calculate width of a single quality interval value
01537 
01538     double interval_value  = 0.0;
01539     int max_interval_width = 0;
01540     std::stringstream str_stream;
01541     std::string interval_string;
01542     for( i = 0; i < histogram.size(); ++i )
01543     {
01544         interval_value = min + (i)*step;
01545         if( step >= 0.001 ) interval_value = round_to_3_significant_digits( interval_value );
01546         str_stream.clear();
01547         str_stream.str( "" );
01548         interval_string = "";
01549         str_stream << interval_value;
01550         interval_string = str_stream.str();
01551         if( interval_string.length() > (size_t)max_interval_width ) max_interval_width = interval_string.length();
01552     }
01553 
01554     // adjust graph width for actual size of interval values
01555     GRAPHW = GRAPHW - ( max_interval_width * 2 ) - 5;
01556 
01557     // For each interval of histogram
01558     for( i = 0; i < histogram.size(); ++i )
01559     {
01560         // First value is the count of the number of values that
01561         // were below the minimum value of the histogram.
01562         if( 0 == i )
01563         {
01564             if( 0 == histogram[i] ) continue;
01565             stream << indent << std::setw( max_interval_width ) << "under min";
01566         }
01567         // Last value is the count of the number of values that
01568         // were above the maximum value of the histogram.
01569         else if( i + 1 == histogram.size() )
01570         {
01571             if( 0 == histogram[i] ) continue;
01572             stream << indent << std::setw( max_interval_width ) << "over max";
01573         }
01574         // Anything else is a valid interval of the histogram.
01575         // Print the range for each interval.
01576         else
01577         {
01578             double start_value = min + ( i - 1 ) * step;
01579             double end_value   = min + (i)*step;
01580             if( step >= 0.001 )
01581             {
01582                 start_value = round_to_3_significant_digits( start_value );
01583                 end_value   = round_to_3_significant_digits( end_value );
01584             }
01585 
01586             stream << indent << "(" << std::setw( max_interval_width ) << std::right << start_value << "-"
01587                    << std::setw( max_interval_width ) << std::left << end_value << ") |";
01588 
01589             // reset stream alignment to right (the default)
01590             stream << std::right;
01591         }
01592 
01593         // Print bar graph
01594 
01595         // First calculate the number of characters to output
01596         int num_graph;
01597         if( log_plot )
01598             num_graph = GRAPHW * (int)log10( (double)( 1 + histogram[i] ) ) / max_bin_value;
01599         else
01600             num_graph = GRAPHW * histogram[i] / max_bin_value;
01601 
01602         // print num_graph characters using array of fill characters.
01603         graph_chars[num_graph] = '\0';
01604         stream << arrptr( graph_chars );
01605         graph_chars[num_graph] = GRAPH_CHAR;
01606 
01607         // Print interval count.
01608         stream << histogram[i] << std::endl;
01609     }
01610     stream << "  metric was evaluated " << count << " times." << std::endl << std::endl;
01611 }
01612 
01613 #ifdef _WIN32
01614 #define fileno( A ) _fileno( ( A ) )
01615 #endif
01616 int QualityAssessor::get_terminal_width() const
01617 {
01618 #ifdef TIOCGWINSZ
01619     int fd = -1;
01620     if( &outputStream == &std::cout )
01621         fd = fileno( stdout );
01622     else if( &outputStream == &std::cerr )
01623         fd = fileno( stderr );
01624 #ifdef FSTREAM_HAS_FD
01625     else if( std::ofstream* f = dynamic_cast< std::ofstream* >( &outputStream ) )
01626         fd = f->rdbuf()->fd();
01627 #endif
01628 
01629     if( fd >= 0 )
01630     {
01631         struct winsize ws;
01632         if( ioctl( fd, TIOCGWINSZ, &ws ) >= 0 ) return ws.ws_col;
01633     }
01634 #endif
01635 
01636     return 0;
01637 }
01638 
01639 std::string QualityAssessor::element_name_as_string( int enum_name )
01640 {
01641     std::string str_value = "";
01642 
01643     switch( enum_name )
01644     {
01645         case POLYGON:
01646             str_value.assign( "polygon" );
01647             break;
01648         case TRIANGLE:
01649             str_value.assign( "triangle" );
01650             break;
01651         case QUADRILATERAL:
01652             str_value.assign( "quadrilateral" );
01653             break;
01654         case POLYHEDRON:
01655             str_value.assign( "polyhedron" );
01656             break;
01657         case TETRAHEDRON:
01658             str_value.assign( "tetrahedron" );
01659             break;
01660         case HEXAHEDRON:
01661             str_value.assign( "hexahedron" );
01662             break;
01663         case PRISM:
01664             str_value.assign( "prism" );
01665             break;
01666         case SEPTAHEDRON:
01667             str_value.assign( "septahedron" );
01668             break;
01669         case MIXED:
01670             str_value.assign( "mixed" );
01671             break;
01672         case PYRAMID:
01673             str_value.assign( "pyramid" );
01674             break;
01675     }
01676 
01677     return str_value;
01678 }
01679 
01680 double QualityAssessor::round_to_3_significant_digits( double number )
01681 {
01682     double result = number;
01683     double p      = 1.0;
01684 
01685     if( number > 0 && number < 1000 )
01686     {
01687         if( number < 1000 )
01688         {
01689             while( number < 100 )
01690             {
01691                 number *= 10;
01692                 p *= 10;
01693             }
01694         }
01695         else
01696         {
01697             while( number >= 1000 )
01698             {
01699                 number /= 10;
01700                 p /= 10;
01701             }
01702         }
01703         int z  = int( number + 0.5 );
01704         result = z / p;
01705     }
01706 
01707     return result;
01708 }
01709 
01710 double QualityAssessor::Assessor::stopping_function_value() const
01711 {
01712     return have_power_mean() ? get_power_mean() : get_average();
01713 }
01714 
01715 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines