MOAB: Mesh Oriented datABase  (version 5.2.1)
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 <math.h>
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                 { new_initial_histogram[i] += ( *initial )->histogram[j]; }
00895             }
00896         }
00897     }
00898 
00899     // put the optimal quality values into the correct new bins
00900     std::vector< int > new_optimal_histogram;
00901     new_optimal_histogram.resize( ( *optimal )->histogram.size(), 0 );
00902     new_optimal_histogram[0] = ( *optimal )->histogram[0];
00903     new_optimal_histogram[new_optimal_histogram.size() - 1] =
00904         ( *optimal )->histogram[( *optimal )->histogram.size() - 1];
00905 
00906     // populate optimal histogram
00907     if( optimal_range == combined_range )
00908     {
00909         // just copy histogram
00910         new_optimal_histogram = ( *optimal )->histogram;
00911     }
00912     else
00913     {
00914         for( size_t i = 1; i < new_optimal_histogram.size() - 1; i++ )
00915         {
00916             double combined_bin_value = combined_min + ( combined_step * ( i - 1 ) );
00917             for( size_t j = 1; j < new_optimal_histogram.size() - 1; j++ )
00918             {
00919                 double optimal_bin_value = optimal_min + ( optimal_step * ( j - 1 ) );
00920                 if( optimal_bin_value >= combined_bin_value && optimal_bin_value < combined_bin_value + combined_step )
00921                 { new_optimal_histogram[i] += ( *optimal )->histogram[j]; }
00922             }
00923         }
00924     }
00925 
00926     // determine largest number of values in a 'bin' for both histograms
00927     unsigned i;
00928     int max_interval_num = 1;
00929     for( i = 0; i < new_initial_histogram.size(); ++i )
00930     {
00931         if( new_initial_histogram[i] > max_interval_num ) max_interval_num = new_initial_histogram[i];
00932     }
00933     for( i = 0; i < new_optimal_histogram.size(); ++i )
00934     {
00935         if( new_optimal_histogram[i] > max_interval_num ) max_interval_num = new_optimal_histogram[i];
00936     }
00937 
00938     // calculate how many bar graph characters will represent the
00939     // largest 'bin' value.
00940     // create the 'before' histogram
00941     int termwidth         = get_terminal_width();
00942     const char indent[]   = "   ";
00943     const char GRAPH_CHAR = '=';                              // Character used to create bar graphs
00944     const int TOTAL_WIDTH = termwidth > 30 ? termwidth : 70;  // Width of histogram
00945     int GRAPHW            = TOTAL_WIDTH - sizeof( indent );
00946 
00947     if( 0 == max_interval_num ) return;  // no data
00948 
00949     // Calculate width of field containing counts for
00950     // histogram intervals (log10(max_interval)).
00951     int num_width = 1;
00952     for( int temp = max_interval_num; temp > 0; temp /= 10 )
00953         ++num_width;
00954     GRAPHW -= num_width;
00955 
00956     // Create an array of bar graph characters for use in output
00957     std::vector< char > graph_chars( GRAPHW + 1, GRAPH_CHAR );
00958 
00959     // Check if bar-graph should be linear or log10 plot
00960     // Do log plot if standard deviation is less that 1.5
00961     // histogram intervals.
00962 
00963     bool log_plot = false;
00964     double stddev = ( *initial )->get_stddev();
00965     if( stddev > 0 && stddev < 2.0 * combined_step )
00966     {
00967         int new_interval = (int)( log10( (double)( 1 + max_interval_num ) ) );
00968         if( new_interval > 0 )
00969         {
00970             log_plot         = true;
00971             max_interval_num = new_interval;
00972         }
00973     }
00974     // output the 'before' histogram
00975 
00976     // Write title
00977     outputStream << std::endl << "************** Common-scale Histograms **************" << std::endl << std::endl;
00978     outputStream << indent << ( *initial )->get_label() << " histogram (initial mesh):";
00979     if( log_plot ) outputStream << " (log10 plot)";
00980     outputStream << std::endl;
00981 
00982     // Calculate width of a single quality interval value
00983     double interval_value  = 0.0;
00984     int max_interval_width = 0;
00985     std::stringstream str_stream;
00986     std::string interval_string;
00987     for( i = 0; i < new_initial_histogram.size(); ++i )
00988     {
00989         interval_value = combined_min + (i)*combined_step;
00990         if( combined_step >= .001 ) interval_value = round_to_3_significant_digits( interval_value );
00991         str_stream.clear();
00992         str_stream.str( "" );
00993         interval_string = "";
00994         str_stream << interval_value;
00995         interval_string = str_stream.str();
00996         if( interval_string.length() > (size_t)max_interval_width ) max_interval_width = interval_string.length();
00997     }
00998 
00999     // adjust graph width for actual size of interval values
01000     GRAPHW = GRAPHW - ( max_interval_width * 2 ) - 5;
01001 
01002     // For each interval of histogram
01003     for( i = 0; i < new_initial_histogram.size(); ++i )
01004     {
01005         // First value is the count of the number of values that
01006         // were below the minimum value of the histogram.
01007         if( 0 == i )
01008         {
01009             if( 0 == new_initial_histogram[i] ) continue;
01010             outputStream << indent << std::setw( max_interval_width ) << "under min";
01011         }
01012         // Last value is the count of the number of values that
01013         // were above the maximum value of the histogram.
01014         else if( i + 1 == new_initial_histogram.size() )
01015         {
01016             if( 0 == new_initial_histogram[i] ) continue;
01017             outputStream << indent << std::setw( max_interval_width ) << "over max";
01018         }
01019         // Anything else is a valid interval of the histogram.
01020         // Print the range for each interval.
01021         else
01022         {
01023             double start_value = combined_min + ( i - 1 ) * combined_step;
01024             double end_value   = combined_min + (i)*combined_step;
01025 
01026             if( combined_step >= 0.001 )
01027             {
01028                 start_value = round_to_3_significant_digits( start_value );
01029                 end_value   = round_to_3_significant_digits( end_value );
01030             }
01031 
01032             outputStream << indent << "(" << std::setw( max_interval_width ) << std::right << start_value << "-"
01033                          << std::setw( max_interval_width ) << std::left << end_value << ") |";
01034         }
01035 
01036         // Print bar graph
01037 
01038         // First calculate the number of characters to output
01039         int num_graph;
01040         if( log_plot )
01041             num_graph = GRAPHW * (int)log10( (double)( 1 + new_initial_histogram[i] ) ) / max_interval_num;
01042         else
01043             num_graph = GRAPHW * new_initial_histogram[i] / max_interval_num;
01044 
01045         // print num_graph characters using array of fill characters.
01046         graph_chars[num_graph] = '\0';
01047         outputStream << arrptr( graph_chars );
01048         graph_chars[num_graph] = GRAPH_CHAR;
01049 
01050         // Print interval count.
01051         outputStream << new_initial_histogram[i] << std::endl;
01052     }
01053 
01054     outputStream << "  metric was evaluated " << ( *initial )->count << " times." << std::endl << std::endl;
01055 
01056     // output the 'after' histogram
01057     outputStream << std::endl << indent << ( *optimal )->get_label() << " histogram (optimal mesh):";
01058     if( log_plot ) outputStream << " (log10 plot)";
01059     outputStream << std::endl;
01060 
01061     // For each interval of histogram
01062     for( i = 0; i < new_optimal_histogram.size(); ++i )
01063     {
01064         // First value is the count of the number of values that
01065         // were below the minimum value of the histogram.
01066         if( 0 == i )
01067         {
01068             if( 0 == new_optimal_histogram[i] ) continue;
01069             outputStream << indent << std::setw( max_interval_width ) << "under min";
01070         }
01071         // Last value is the count of the number of values that
01072         // were above the maximum value of the histogram.
01073         else if( i + 1 == new_optimal_histogram.size() )
01074         {
01075             if( 0 == new_optimal_histogram[i] ) continue;
01076             outputStream << indent << std::setw( max_interval_width ) << "over max";
01077         }
01078         // Anything else is a valid interval of the histogram.
01079         // Print the range for each interval.
01080         else
01081         {
01082             double start_value = combined_min + ( i - 1 ) * combined_step;
01083             double end_value   = combined_min + (i)*combined_step;
01084 
01085             if( combined_step >= 0.001 )
01086             {
01087                 start_value = round_to_3_significant_digits( start_value );
01088                 end_value   = round_to_3_significant_digits( end_value );
01089             }
01090 
01091             outputStream << indent << "(" << std::setw( max_interval_width ) << std::right << start_value << "-"
01092                          << std::setw( max_interval_width ) << std::left << end_value << ") |";
01093         }
01094 
01095         // Print bar graph
01096 
01097         // First calculate the number of characters to output
01098         int num_graph;
01099         if( log_plot )
01100             num_graph = GRAPHW * (int)log10( (double)( 1 + new_optimal_histogram[i] ) ) / max_interval_num;
01101         else
01102             num_graph = GRAPHW * new_optimal_histogram[i] / max_interval_num;
01103 
01104         // print num_graph characters using array of fill characters.
01105         graph_chars[num_graph] = '\0';
01106         outputStream << arrptr( graph_chars );
01107         graph_chars[num_graph] = GRAPH_CHAR;
01108 
01109         // Print interval count.
01110         outputStream << new_optimal_histogram[i] << std::endl;
01111     }
01112 
01113     outputStream << "  metric was evaluated " << ( *optimal )->count << " times." << std::endl << std::endl;
01114 
01115     return;
01116 }
01117 
01118 QualityAssessor::Assessor::Assessor( QualityMetric* metric, const char* label )
01119     : qualMetric( metric ), mLabel( label ? std::string( label ) : metric->get_name() ), pMean( 0.0 ),
01120       haveHistRange( false ), histMin( 1.0 ), histMax( 0.0 ), tagHandle( 0 ), stoppingFunction( false ),
01121       referenceCount( 0 ), assessScheme( NO_SCHEME )
01122 {
01123     reset_data();
01124 }
01125 
01126 const QualityAssessor::Assessor* QualityAssessor::get_results( QualityMetric* metric ) const
01127 {
01128     list_type::const_iterator iter;
01129     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01130         if( ( *iter )->get_metric() == metric ) return *iter;
01131     return 0;
01132 }
01133 
01134 const QualityAssessor::Assessor* QualityAssessor::get_results( const char* name ) const
01135 {
01136     list_type::const_iterator iter;
01137     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01138         if( ( *iter )->get_label() == name ) return *iter;
01139     return 0;
01140 }
01141 
01142 void QualityAssessor::Assessor::get_histogram( double& lower_bound_out, double& upper_bound_out,
01143                                                std::vector< int >& counts_out, MsqError& err ) const
01144 {
01145     if( !have_histogram() )
01146     {
01147         MSQ_SETERR( err )( "No histogram calculated.", MsqError::INVALID_STATE );
01148         return;
01149     }
01150 
01151     lower_bound_out = histMin;
01152     upper_bound_out = histMax;
01153     counts_out      = histogram;
01154 }
01155 
01156 void QualityAssessor::Assessor::reset_data()
01157 {
01158     count        = 0;
01159     sum          = 0;
01160     maximum      = -HUGE_VAL;
01161     minimum      = HUGE_VAL;
01162     sqrSum       = 0;
01163     pSum         = 0;
01164     numInvalid   = 0;
01165     assessScheme = NO_SCHEME;
01166     // zero histogram data
01167     size_t hist_size = histogram.size();
01168     histogram.clear();
01169     histogram.resize( hist_size, 0 );
01170 }
01171 
01172 void QualityAssessor::Assessor::add_value( double metric_value )
01173 {
01174     sum += metric_value;
01175     sqrSum += metric_value * metric_value;
01176     if( metric_value > maximum ) maximum = metric_value;
01177     if( metric_value < minimum ) minimum = metric_value;
01178     // Only add value to histogram data from this function if
01179     // the user has specified the range.  If user has not
01180     // specified the range, QualityAssessor will call add_hist_value()
01181     // directly once the range has been calculated.
01182     if( have_histogram() && haveHistRange ) add_hist_value( metric_value );
01183 
01184     if( have_power_mean() ) pSum += pow( metric_value, pMean );
01185 
01186     ++count;
01187 }
01188 
01189 void QualityAssessor::Assessor::add_invalid_value()
01190 {
01191     ++numInvalid;
01192 }
01193 
01194 void QualityAssessor::Assessor::add_hist_value( double metric_value )
01195 {
01196     // First and last values in array are counts of values
01197     // outside the user-specified range of the histogram
01198     // (below and above, respectively.)
01199     if( metric_value < histMin )
01200         ++histogram[0];
01201     else if( metric_value > histMax )
01202         ++histogram[histogram.size() - 1];
01203     else
01204     {
01205         // Calculate which interval the value is in.  Add one
01206         // because first entry is for values below user-specifed
01207         // minimum value for histogram.
01208         double range = histMax - histMin;
01209         double fract;
01210         if( range > DBL_EPSILON )
01211             fract = ( metric_value - histMin ) / range;
01212         else
01213             fract = 0.0;
01214         unsigned cell;
01215         if( fabs( fract - 1.0 ) < histMax * DBL_EPSILON )
01216             cell = histogram.size() - 1;
01217         else
01218             cell = 1 + (int)( ( fract * ( histogram.size() - 2 ) ) );
01219 
01220         // Add value to interval.
01221         ++histogram[cell];
01222     }
01223 }
01224 
01225 void QualityAssessor::Assessor::calculate_histogram_range()
01226 {
01227     double lower_bound = minimum;
01228     int num_intervals  = histogram.size();
01229     double step        = ( maximum - lower_bound ) / num_intervals;
01230     if( step == 0 ) step = 1.0;
01231     double size = pow( 10.0, floor( log10( step / ( num_intervals - 1 ) ) ) );
01232     if( size < 1e-6 ) size = 1.0;
01233     histMin = lower_bound;
01234     histMax = lower_bound + num_intervals * size * ceil( step / size );
01235 }
01236 
01237 void QualityAssessor::print_summary( std::ostream& stream ) const
01238 {
01239     const int NAMEW = 19;  // Width of name column in table output
01240     const int NUMW  = 12;  // Width of value columns in table output
01241 
01242     // Print title
01243     stream << std::endl
01244            << "************** " << qualityAssessorName << " Summary **************" << std::endl
01245            << std::endl;
01246 
01247     if( myData->freeElementCount != myData->elementCount )
01248         stream << "  Evaluating quality for " << myData->freeElementCount << " free elements of "
01249                << myData->elementCount << " total elements." << std::endl;
01250     else
01251         stream << "  Evaluating quality for " << myData->elementCount << " elements." << std::endl;
01252 
01253     // Print element type totals
01254     std::string str_value = "";
01255     for( int i = POLYGON; i <= MIXED; i++ )
01256     {
01257         if( elementTypeCount[i - POLYGON] )
01258         {
01259             str_value = element_name_as_string( i );
01260             if( str_value != "" )
01261                 stream << "  This mesh had " << elementTypeCount[i - POLYGON] << " " << str_value << " elements."
01262                        << std::endl;
01263         }
01264     }
01265 
01266     if( myData->invertedElementCount )
01267     {
01268         stream << "  THERE ARE " << myData->invertedElementCount << " INVERTED ELEMENTS. " << std::endl
01269                << "  (Elements invalid at " << myData->invertedSampleCount << " of " << myData->sampleCount
01270                << " sample locations.)" << std::endl
01271                << std::endl;
01272     }
01273     else
01274     {
01275         stream << "  There were no inverted elements detected. " << std::endl;
01276     }
01277 
01278     // Check if there are invalid values for any metric
01279     list_type::const_iterator iter;
01280     bool some_invalid = false;
01281     for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01282     {
01283         if( ( *iter )->get_invalid_element_count() )
01284         {
01285             some_invalid = true;
01286             stream << "  " << ( *iter )->get_invalid_element_count() << " OF " << ( *iter )->get_count()
01287                    << " ENTITIES EVALUATED TO AN UNDEFINED VALUE FOR " << ( *iter )->get_label() << std::endl
01288                    << std::endl;
01289         }
01290     }
01291     if( !some_invalid )
01292     { stream << "  No entities had undefined values for any computed metric." << std::endl << std::endl; }
01293 
01294     if( !invalid_values )
01295     {
01296         // Check if a user-define power-mean was calculated for any of the metrics
01297         std::set< double > pmeans;
01298         for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01299             if( ( *iter )->have_power_mean() ) pmeans.insert( ( *iter )->get_power() );
01300 
01301         // If power-mean of 1 or 2 was requested, change titles rather
01302         // than printing redundant values
01303         std::string average_str( "average" ), rms_str( "rms" );
01304         if( pmeans.find( 1.0 ) != pmeans.end() )
01305         {
01306             pmeans.erase( pmeans.find( 1.0 ) );
01307             average_str = "1-mean";
01308         }
01309         if( pmeans.find( 2.0 ) != pmeans.end() )
01310         {
01311             pmeans.erase( pmeans.find( 2.0 ) );
01312             rms_str = "2-mean";
01313         }
01314 
01315         // Number of values in table
01316         unsigned num_values = pmeans.size() + 5;
01317 
01318         // Decide how wide of a table field should be used for the metric name
01319         unsigned twidth    = get_terminal_width();
01320         unsigned maxnwidth = NAMEW;
01321         if( twidth )
01322         {
01323             unsigned valwidth = NUMW * num_values;
01324             maxnwidth         = valwidth < twidth ? twidth - valwidth : 0;
01325         }
01326         unsigned namewidth = 0;
01327         for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01328             if( ( *iter )->get_label().size() > namewidth ) namewidth = ( *iter )->get_label().size();
01329         if( namewidth > maxnwidth ) namewidth = maxnwidth;
01330         if( namewidth < 7 )  // at least enough width for the column header
01331             namewidth = 7;
01332 
01333         int number_of_assessments = 0;
01334 
01335         // print metric values
01336         for( iter = assessList.begin(); iter != assessList.end(); ++iter )
01337         {
01338             if( number_of_assessments > 0 ) stream << "    -------------------------------------------" << std::endl;
01339 
01340             // print assessment method used to calculate the statistics
01341             if( ( *iter )->assessScheme == TMP_QUALITY_METRIC )
01342                 stream << "     Sample Point Quality Statistics" << std::endl << std::endl;
01343             else
01344                 stream << "     Element Quality Statistics" << std::endl << std::endl;
01345 
01346             // print comlumn label line
01347             std::set< double >::const_iterator piter;
01348             stream << std::setw( NUMW ) << "minimum";
01349             for( piter = pmeans.begin(); piter != pmeans.end() && *piter < 1.0; ++piter )
01350                 stream << std::setw( NUMW - 6 ) << *piter << "-mean ";
01351             stream << std::setw( NUMW ) << average_str;
01352             for( ; piter != pmeans.end() && *piter < 2.0; ++piter )
01353                 stream << std::setw( NUMW - 6 ) << *piter << "-mean ";
01354             stream << std::setw( NUMW ) << rms_str;
01355             for( ; piter != pmeans.end(); ++piter )
01356                 stream << std::setw( NUMW - 6 ) << *piter << "-mean ";
01357             stream << std::setw( NUMW ) << "maximum";
01358             stream << std::setw( NUMW ) << "std.dev.";
01359             stream << std::endl;
01360 
01361             stream << std::setw( NUMW ) << ( *iter )->get_minimum();
01362             // print power-means with P less than 1.0
01363             for( piter = pmeans.begin(); piter != pmeans.end() && *piter < 1.0; ++piter )
01364             {
01365                 if( ( *iter )->have_power_mean() && ( *iter )->get_power() == *piter )
01366                     stream << std::setw( NUMW ) << ( *iter )->get_power_mean();
01367                 else
01368                     stream << std::setw( NUMW ) << " ";
01369             }
01370             // print average
01371             stream << std::setw( NUMW ) << ( *iter )->get_average();
01372             // print power-means with P less than 2.0
01373             for( ; piter != pmeans.end() && *piter < 2.0; ++piter )
01374             {
01375                 if( ( *iter )->have_power_mean() && ( *iter )->get_power() == *piter )
01376                     stream << std::setw( NUMW ) << ( *iter )->get_power_mean();
01377                 else
01378                     stream << std::setw( NUMW ) << " ";
01379             }
01380             // print RMS
01381             stream << std::setw( NUMW ) << ( *iter )->get_rms();
01382             // print power-means with P greater than 2.0
01383             for( ; piter != pmeans.end(); ++piter )
01384             {
01385                 if( ( *iter )->have_power_mean() && ( *iter )->get_power() == *piter )
01386                     stream << std::setw( NUMW ) << ( *iter )->get_power_mean();
01387                 else
01388                     stream << std::setw( NUMW ) << " ";
01389             }
01390             // print maximum and standard deviation
01391             stream << std::setw( NUMW ) << ( *iter )->get_maximum();
01392             stream << std::setw( NUMW ) << ( *iter )->get_stddev();
01393             stream << std::endl << std::endl;
01394 
01395             stream << "     Number of statistics = " << ( *iter )->get_count() << std::endl;
01396 
01397             // print name
01398             stream << "     Metric = " << ( *iter )->get_label() << std::endl;
01399 
01400             // Output the method used to calcualte the quality values
01401             switch( ( *iter )->assessScheme )
01402             {
01403                 case ELEMENT_AVG_QM:
01404                     stream << "     Element Quality = average over metric values at the elements' "
01405                               "sample points"
01406                            << std::endl
01407                            << std::endl;
01408                     break;
01409                 case ELEMENT_MAX_QM:
01410                     stream << "     Element Quality = maximum over metric values at the elements' "
01411                               "sample points"
01412                            << std::endl
01413                            << std::endl;
01414                     break;
01415                 case TMP_QUALITY_METRIC:
01416                     stream << std::endl << std::endl;
01417                     break;
01418                 case QUALITY_METRIC:
01419                     stream << "     Element Quality not based on sample points." << std::endl << std::endl;
01420                     break;
01421                 default:
01422                     stream << "     Scheme used for deriving qualitiy values unknown" << std::endl << std::endl;
01423             }
01424             if( ( *iter )->have_histogram() ) ( *iter )->print_histogram( stream, get_terminal_width() );
01425             number_of_assessments++;
01426         }
01427     }
01428     else
01429     {
01430         stream << "  Element Quality Statistics are Undefined for this Metric because" << std::endl
01431                << "  there are Inverted Elements" << std::endl
01432                << std::endl;
01433     }
01434 }
01435 
01436 void QualityAssessor::Assessor::print_histogram( std::ostream& stream, int termwidth ) const
01437 {
01438     // Portability notes:
01439     //  Use log10 rather than log10f because the float variations require
01440     //  including platform-dependent headers on some platforms.
01441     //  Explicitly cast log10 argument to double because some platforms
01442     //  have overloaded float and double variations in C++ making an
01443     //  implicit cast from an integer ambiguous.
01444 
01445     const char indent[]   = "   ";
01446     const char GRAPH_CHAR = '=';                              // Character used to create bar graphs
01447     const int TOTAL_WIDTH = termwidth > 30 ? termwidth : 70;  // Width of histogram
01448     int GRAPHW            = TOTAL_WIDTH - sizeof( indent );
01449 
01450     // range is either user-specified (histMin & histMax) or
01451     // calculated (minimum & maximum)
01452     double min, max;
01453     min = histMin;
01454     max = histMax;
01455     // Witdh of one interval of histogram
01456     double step = ( max - min ) / ( histogram.size() - 2 );
01457     // round step to 3 significant digits
01458 
01459     if( step >= 0.001 ) step = round_to_3_significant_digits( step );
01460 
01461     // Find maximum value for an interval bin of the histogram
01462     unsigned i;
01463     int max_bin_value = 1;
01464     for( i = 0; i < histogram.size(); ++i )
01465         if( histogram[i] > max_bin_value ) max_bin_value = histogram[i];
01466 
01467     if( 0 == max_bin_value ) return;  // no data
01468 
01469     // Calculate width of field containing counts for
01470     // histogram intervals (log10(max_bin_value)).
01471     int num_width = 1;
01472     for( int temp = max_bin_value; temp > 0; temp /= 10 )
01473         ++num_width;
01474     GRAPHW -= num_width;
01475 
01476     // Create an array of bar graph characters for use in output
01477     std::vector< char > graph_chars( GRAPHW + 1, GRAPH_CHAR );
01478 
01479     // Check if bar-graph should be linear or log10 plot
01480     // Do log plot if standard deviation is less that 1.5
01481     // histogram intervals.
01482     bool log_plot = false;
01483     double stddev = get_stddev();
01484     if( stddev > 0 && stddev < 2.0 * step )
01485     {
01486         int new_interval = (int)( log10( (double)( 1 + max_bin_value ) ) );
01487         if( new_interval > 0 )
01488         {
01489             log_plot      = true;
01490             max_bin_value = new_interval;
01491         }
01492     }
01493 
01494     // Write title
01495     stream << indent << get_label() << " histogram:";
01496     if( log_plot ) stream << " (log10 plot)";
01497     stream << std::endl;
01498 
01499     // Calculate width of a single quality interval value
01500 
01501     double interval_value  = 0.0;
01502     int max_interval_width = 0;
01503     std::stringstream str_stream;
01504     std::string interval_string;
01505     for( i = 0; i < histogram.size(); ++i )
01506     {
01507         interval_value = min + (i)*step;
01508         if( step >= 0.001 ) interval_value = round_to_3_significant_digits( interval_value );
01509         str_stream.clear();
01510         str_stream.str( "" );
01511         interval_string = "";
01512         str_stream << interval_value;
01513         interval_string = str_stream.str();
01514         if( interval_string.length() > (size_t)max_interval_width ) max_interval_width = interval_string.length();
01515     }
01516 
01517     // adjust graph width for actual size of interval values
01518     GRAPHW = GRAPHW - ( max_interval_width * 2 ) - 5;
01519 
01520     // For each interval of histogram
01521     for( i = 0; i < histogram.size(); ++i )
01522     {
01523         // First value is the count of the number of values that
01524         // were below the minimum value of the histogram.
01525         if( 0 == i )
01526         {
01527             if( 0 == histogram[i] ) continue;
01528             stream << indent << std::setw( max_interval_width ) << "under min";
01529         }
01530         // Last value is the count of the number of values that
01531         // were above the maximum value of the histogram.
01532         else if( i + 1 == histogram.size() )
01533         {
01534             if( 0 == histogram[i] ) continue;
01535             stream << indent << std::setw( max_interval_width ) << "over max";
01536         }
01537         // Anything else is a valid interval of the histogram.
01538         // Print the range for each interval.
01539         else
01540         {
01541             double start_value = min + ( i - 1 ) * step;
01542             double end_value   = min + (i)*step;
01543             if( step >= 0.001 )
01544             {
01545                 start_value = round_to_3_significant_digits( start_value );
01546                 end_value   = round_to_3_significant_digits( end_value );
01547             }
01548 
01549             stream << indent << "(" << std::setw( max_interval_width ) << std::right << start_value << "-"
01550                    << std::setw( max_interval_width ) << std::left << end_value << ") |";
01551 
01552             // reset stream alignment to right (the default)
01553             stream << std::right;
01554         }
01555 
01556         // Print bar graph
01557 
01558         // First calculate the number of characters to output
01559         int num_graph;
01560         if( log_plot )
01561             num_graph = GRAPHW * (int)log10( (double)( 1 + histogram[i] ) ) / max_bin_value;
01562         else
01563             num_graph = GRAPHW * histogram[i] / max_bin_value;
01564 
01565         // print num_graph characters using array of fill characters.
01566         graph_chars[num_graph] = '\0';
01567         stream << arrptr( graph_chars );
01568         graph_chars[num_graph] = GRAPH_CHAR;
01569 
01570         // Print interval count.
01571         stream << histogram[i] << std::endl;
01572     }
01573     stream << "  metric was evaluated " << count << " times." << std::endl << std::endl;
01574 }
01575 
01576 #ifdef _WIN32
01577 #define fileno( A ) _fileno( ( A ) )
01578 #endif
01579 int QualityAssessor::get_terminal_width() const
01580 {
01581 #ifdef TIOCGWINSZ
01582     int fd = -1;
01583     if( &outputStream == &std::cout )
01584         fd = fileno( stdout );
01585     else if( &outputStream == &std::cerr )
01586         fd = fileno( stderr );
01587 #ifdef FSTREAM_HAS_FD
01588     else if( std::ofstream* f = dynamic_cast< std::ofstream* >( &outputStream ) )
01589         fd = f->rdbuf()->fd();
01590 #endif
01591 
01592     if( fd >= 0 )
01593     {
01594         struct winsize ws;
01595         if( ioctl( fd, TIOCGWINSZ, &ws ) >= 0 ) return ws.ws_col;
01596     }
01597 #endif
01598 
01599     return 0;
01600 }
01601 
01602 std::string QualityAssessor::element_name_as_string( int enum_name )
01603 {
01604     std::string str_value = "";
01605 
01606     switch( enum_name )
01607     {
01608         case POLYGON:
01609             str_value.assign( "polygon" );
01610             break;
01611         case TRIANGLE:
01612             str_value.assign( "triangle" );
01613             break;
01614         case QUADRILATERAL:
01615             str_value.assign( "quadrilateral" );
01616             break;
01617         case POLYHEDRON:
01618             str_value.assign( "polyhedron" );
01619             break;
01620         case TETRAHEDRON:
01621             str_value.assign( "tetrahedron" );
01622             break;
01623         case HEXAHEDRON:
01624             str_value.assign( "hexahedron" );
01625             break;
01626         case PRISM:
01627             str_value.assign( "prism" );
01628             break;
01629         case SEPTAHEDRON:
01630             str_value.assign( "septahedron" );
01631             break;
01632         case MIXED:
01633             str_value.assign( "mixed" );
01634             break;
01635         case PYRAMID:
01636             str_value.assign( "pyramid" );
01637             break;
01638     }
01639 
01640     return str_value;
01641 }
01642 
01643 double QualityAssessor::round_to_3_significant_digits( double number )
01644 {
01645     double result = number;
01646     double p      = 1.0;
01647 
01648     if( number > 0 && number < 1000 )
01649     {
01650         if( number < 1000 )
01651         {
01652             while( number < 100 )
01653             {
01654                 number *= 10;
01655                 p *= 10;
01656             }
01657         }
01658         else
01659         {
01660             while( number >= 1000 )
01661             {
01662                 number /= 10;
01663                 p /= 10;
01664             }
01665         }
01666         int z  = int( number + 0.5 );
01667         result = z / p;
01668     }
01669 
01670     return result;
01671 }
01672 
01673 double QualityAssessor::Assessor::stopping_function_value() const
01674 {
01675     return have_power_mean() ? get_power_mean() : get_average();
01676 }
01677 
01678 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines