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