MOAB: Mesh Oriented datABase  (version 5.4.0)
QualityAssessorTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2007 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to 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     (2007) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file QualityAssessorTest.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "QualityAssessor.hpp"
00034 #include "AspectRatioGammaQualityMetric.hpp"
00035 #include "LocalSizeQualityMetric.hpp"
00036 #include "ConditionNumberQualityMetric.hpp"
00037 #include "TShapeNB1.hpp"
00038 #include "IdealShapeTarget.hpp"
00039 #include "TQualityMetric.hpp"
00040 #include "PlanarDomain.hpp"
00041 #include "Settings.hpp"
00042 #include "MeshImpl.hpp"
00043 #include "PatchData.hpp"
00044 #include "MsqMeshEntity.hpp"
00045 #include "ArrayMesh.hpp"
00046 #include "InstructionQueue.hpp"
00047 
00048 #include "UnitUtil.hpp"
00049 
00050 #include <algorithm>
00051 #include <numeric>
00052 #include <iomanip>
00053 
00054 using namespace MBMesquite;
00055 using namespace std;
00056 
00057 class QualityAssessorTest : public CppUnit::TestFixture
00058 {
00059   private:
00060     CPPUNIT_TEST_SUITE( QualityAssessorTest );
00061     CPPUNIT_TEST( test_basic_stats_element );
00062     CPPUNIT_TEST( test_basic_stats_vertex );
00063     CPPUNIT_TEST( test_basic_stats_sample );
00064     CPPUNIT_TEST( test_histogram_known_range );
00065     CPPUNIT_TEST( test_histogram_unknown_range );
00066     CPPUNIT_TEST( test_power_mean );
00067     CPPUNIT_TEST( test_invalid_count );
00068     CPPUNIT_TEST( test_inverted_count );
00069     CPPUNIT_TEST( test_output_control );
00070     CPPUNIT_TEST( test_tag_element );
00071     CPPUNIT_TEST( test_tag_vertex );
00072     CPPUNIT_TEST( test_tag_inverted );
00073     CPPUNIT_TEST( test_print_inverted );
00074     CPPUNIT_TEST( test_print_stats );
00075     CPPUNIT_TEST( test_print_name );
00076     CPPUNIT_TEST( test_modify_metric );
00077     CPPUNIT_TEST( test_free_only );
00078     CPPUNIT_TEST_SUITE_END();
00079 
00080     vector< double > vertCoords, invertCoords;
00081     vector< int > fixedFlags;
00082     vector< unsigned long > triConn, invertConn;
00083 
00084     MeshImpl myMesh, invertedMesh;
00085     PlanarDomain myDomain;
00086     Settings mySettings;
00087 
00088   public:
00089     void setUp();
00090     void tearDown();
00091 
00092   public:
00093     QualityAssessorTest() : myDomain( PlanarDomain::XY ) {}
00094 
00095     void test_basic_stats_element();
00096     void test_basic_stats_vertex();
00097     void test_basic_stats_sample();
00098     void test_histogram_known_range();
00099     void test_histogram_unknown_range();
00100     void test_power_mean();
00101     void test_invalid_count();
00102     void test_inverted_count();
00103     void test_output_control();
00104     void test_tag_element();
00105     void test_tag_vertex();
00106     void test_tag_inverted();
00107     void test_print_inverted();
00108     void test_print_stats();
00109     void test_print_name();
00110     void test_modify_metric();
00111     void test_free_only();
00112 };
00113 
00114 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( QualityAssessorTest, "Unit" );
00115 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( QualityAssessorTest, "QualityAssessorTest" );
00116 
00117 void QualityAssessorTest::setUp()
00118 {
00119     MsqError err;
00120 
00121     // create the following mesh:
00122     /*2.0
00123      *  (3)------------------(2)
00124      *   |\                __/|
00125      *   | \       3    __/ / |
00126      *   |  \        __/   /  |
00127      *   |   \    __/  4  /   |
00128      *   |    \ _/       /    |
00129      *1.0|  0 (4)------(5)  2 |
00130      *   |    /. \__    .\    |
00131      *   |   / .    \__5. \   |
00132      *   |  /  .       \__ \  |
00133      *   | /   .   1    . \_\ |
00134      *   |/    .        .    \|
00135      *  (0)------------------(1)
00136      *0.0     1.0      2.0   3.0
00137      */
00138     const char vtk_file[] = "# vtk DataFile Version 3.0\n"
00139                             "QualityAssessorTest input 1\n"
00140                             "ASCII\n"
00141                             "DATASET UNSTRUCTURED_GRID\n"
00142                             "POINTS 6 float\n"
00143                             "0 0 0\n"
00144                             "3 0 0\n"
00145                             "3 2 0\n"
00146                             "0 2 0\n"
00147                             "1 1 0\n"
00148                             "2 1 0\n"
00149                             "CELLS 6 24\n"
00150                             "3 4 3 0\n"
00151                             "3 0 1 4\n"
00152                             "3 1 2 5\n"
00153                             "3 2 3 4\n"
00154                             "3 4 5 2\n"
00155                             "3 5 4 1\n"
00156                             "CELL_TYPES 6\n"
00157                             "5 5 5 5 5 5\n"
00158                             "POINT_DATA 6\n"
00159                             "SCALARS fixed int\n"
00160                             "LOOKUP_TABLE default\n"
00161                             "1 1 1 1 0 0\n";
00162 
00163     // define a mesh with two triangles, one of which is inverted wrt the other
00164     const char invert_file[] = "# vtk DataFile Version 3.0\n"
00165                                "QualityAssessorTest input 2\n"
00166                                "ASCII\n"
00167                                "DATASET UNSTRUCTURED_GRID\n"
00168                                "POINTS 3 float\n"
00169                                "0 0 0\n"
00170                                "1 0 0\n"
00171                                "0 1 0\n"
00172                                "CELLS 2 8\n"
00173                                "3 0 1 2\n"
00174                                "3 0 2 1\n"
00175                                "CELL_TYPES 2\n"
00176                                "5 5\n"
00177                                "POINT_DATA 3\n"
00178                                "SCALARS fixed int\n"
00179                                "LOOKUP_TABLE default\n"
00180                                "1 1 1\n";
00181 
00182     myMesh.clear();
00183     invertedMesh.clear();
00184 
00185     const char* tmpname = "QualityAssessorTest.vtk";
00186     FILE* filp          = fopen( tmpname, "w" );
00187     CPPUNIT_ASSERT( NULL != filp );
00188     size_t w = fwrite( vtk_file, sizeof( vtk_file ), 1, filp );
00189     fclose( filp );
00190     myMesh.clear();
00191     myMesh.read_vtk( tmpname, err );
00192     remove( tmpname );
00193     CPPUNIT_ASSERT_EQUAL( (size_t)1, w );
00194     ASSERT_NO_ERROR( err );
00195 
00196     filp = fopen( tmpname, "w" );
00197     CPPUNIT_ASSERT( NULL != filp );
00198     w = fwrite( invert_file, sizeof( invert_file ), 1, filp );
00199     fclose( filp );
00200     invertedMesh.clear();
00201     invertedMesh.read_vtk( tmpname, err );
00202     remove( tmpname );
00203     CPPUNIT_ASSERT_EQUAL( (size_t)1, w );
00204     ASSERT_NO_ERROR( err );
00205 }
00206 
00207 void QualityAssessorTest::tearDown() {}
00208 
00209 // Define a decorator (wrapper) for a QualityMetric
00210 // instance that records each metric value.
00211 class MetricLogger : public QualityMetric
00212 {
00213   public:
00214     MetricLogger( QualityMetric* metric ) : invalidCount( 0 ), mMetric( metric ) {}
00215 
00216     vector< double > mValues;
00217     vector< Mesh::EntityHandle > mHandles;
00218     int invalidCount;
00219 
00220     double min() const
00221     {
00222         return *min_element( mValues.begin(), mValues.end() );
00223     }
00224     double max() const
00225     {
00226         return *max_element( mValues.begin(), mValues.end() );
00227     }
00228     double avg() const
00229     {
00230         double x = 0;
00231         return accumulate( mValues.begin(), mValues.end(), x ) / mValues.size();
00232     }
00233     double sqrsum() const
00234     {
00235         double result = 0.0;
00236         for( unsigned i = 0; i < mValues.size(); ++i )
00237             result += mValues[i] * mValues[i];
00238         return result;
00239     }
00240     double rms() const
00241     {
00242         return sqrt( sqrsum() / mValues.size() );
00243     }
00244     double dev() const
00245     {
00246         return sqrt( sqrsum() / mValues.size() - avg() * avg() );
00247     }
00248     double pmean( double p ) const
00249     {
00250         double result = 0.0;
00251         for( unsigned i = 0; i < mValues.size(); ++i )
00252             result += pow( mValues[i], p );
00253         return pow( result / mValues.size(), 1. / p );
00254     }
00255 
00256     bool no_duplicate_evals() const
00257     {
00258         vector< Mesh::EntityHandle > handles( mHandles );
00259         sort( handles.begin(), handles.end() );
00260         return handles.end() == unique( handles.begin(), handles.end() );
00261     }
00262 
00263     MetricType get_metric_type() const
00264     {
00265         return mMetric->get_metric_type();
00266     }
00267 
00268     std::string get_name() const
00269     {
00270         return mMetric->get_name();
00271     }
00272 
00273     int get_negate_flag() const
00274     {
00275         return mMetric->get_negate_flag();
00276     }
00277 
00278     void get_evaluations( PatchData& pd, vector< size_t >& handles, bool free, MsqError& err )
00279     {
00280         return mMetric->get_evaluations( pd, handles, free, err );
00281     }
00282 
00283     bool evaluate( PatchData& pd, size_t handle, double& value, MsqError& err )
00284     {
00285         bool rval = mMetric->evaluate( pd, handle, value, err );
00286         mValues.push_back( value );
00287         if( get_metric_type() == QualityMetric::VERTEX_BASED )
00288             mHandles.push_back( pd.get_vertex_handles_array()[handle] );
00289         else if( handle < pd.num_elements() )
00290             mHandles.push_back( pd.get_element_handles_array()[handle] );
00291         if( !rval ) ++invalidCount;
00292         return rval;
00293     }
00294 
00295     bool evaluate_with_indices( PatchData&, size_t, double&, vector< size_t >&, MsqError& err )
00296     {  // shouldn't be called by QualityAssessor
00297         CPPUNIT_ASSERT( false );
00298         MSQ_SETERR( err )( MsqError::NOT_IMPLEMENTED );
00299         return false;
00300     }
00301     bool evaluate_with_gradient( PatchData&, size_t, double&, vector< size_t >&, vector< Vector3D >&, MsqError& err )
00302     {  // shouldn't be called by QualityAssessor
00303         CPPUNIT_ASSERT( false );
00304         MSQ_SETERR( err )( MsqError::NOT_IMPLEMENTED );
00305         return false;
00306     }
00307     bool evaluate_with_Hessian( PatchData&,
00308                                 size_t,
00309                                 double&,
00310                                 vector< size_t >&,
00311                                 vector< Vector3D >&,
00312                                 vector< Matrix3D >&,
00313                                 MsqError& err )
00314     {  // shouldn't be called by QualityAssessor
00315         CPPUNIT_ASSERT( false );
00316         MSQ_SETERR( err )( MsqError::NOT_IMPLEMENTED );
00317         return false;
00318     }
00319 
00320   private:
00321     QualityMetric* mMetric;
00322 };
00323 
00324 void QualityAssessorTest::test_basic_stats_element()
00325 {
00326     AspectRatioGammaQualityMetric metric;
00327     MetricLogger logger( &metric );
00328     QualityAssessor qa( &logger );
00329     qa.disable_printing_results();
00330 
00331     MsqError err;
00332     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00333     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00334     ASSERT_NO_ERROR( err );
00335 
00336     // check didn't evaluate any element more than once
00337     CPPUNIT_ASSERT( logger.no_duplicate_evals() );
00338 
00339     // check one eval for each element
00340     vector< Mesh::ElementHandle > elems;
00341     myMesh.get_all_elements( elems, err );
00342     ASSERT_NO_ERROR( err );
00343     CPPUNIT_ASSERT_EQUAL( elems.size(), logger.mValues.size() );
00344 
00345     // check values
00346     const QualityAssessor::Assessor* results = qa.get_results( &logger );
00347     CPPUNIT_ASSERT( results != NULL );
00348     CPPUNIT_ASSERT_EQUAL( (int)elems.size(), results->get_count() );
00349     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.min(), results->get_minimum(), 1e-6 );
00350     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.avg(), results->get_average(), 1e-6 );
00351     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.rms(), results->get_rms(), 1e-6 );
00352     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.max(), results->get_maximum(), 1e-6 );
00353     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.dev(), results->get_stddev(), 1e-6 );
00354 }
00355 
00356 void QualityAssessorTest::test_basic_stats_vertex()
00357 {
00358     LocalSizeQualityMetric metric;
00359     MetricLogger logger( &metric );
00360     QualityAssessor qa( &logger );
00361     qa.measure_free_samples_only( false );
00362     qa.disable_printing_results();
00363 
00364     MsqError err;
00365     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00366     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00367     ASSERT_NO_ERROR( err );
00368 
00369     // check didn't evaluate any vertex more than once
00370     CPPUNIT_ASSERT( logger.no_duplicate_evals() );
00371 
00372     // check one eval for each element
00373     vector< Mesh::VertexHandle > verts;
00374     myMesh.get_all_vertices( verts, err );
00375     ASSERT_NO_ERROR( err );
00376     CPPUNIT_ASSERT_EQUAL( verts.size(), logger.mValues.size() );
00377 
00378     // check values
00379     const QualityAssessor::Assessor* results = qa.get_results( &logger );
00380     CPPUNIT_ASSERT( results != NULL );
00381     CPPUNIT_ASSERT_EQUAL( (int)verts.size(), results->get_count() );
00382     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.min(), results->get_minimum(), 1e-6 );
00383     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.avg(), results->get_average(), 1e-6 );
00384     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.rms(), results->get_rms(), 1e-6 );
00385     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.max(), results->get_maximum(), 1e-6 );
00386     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.dev(), results->get_stddev(), 1e-6 );
00387 }
00388 
00389 void QualityAssessorTest::test_basic_stats_sample()
00390 {
00391     TShapeNB1 tm;
00392     IdealShapeTarget tc;
00393     TQualityMetric metric( &tc, &tm );
00394     MetricLogger logger( &metric );
00395     QualityAssessor qa( &logger );
00396     qa.disable_printing_results();
00397 
00398     MsqError err;
00399     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00400     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00401     ASSERT_NO_ERROR( err );
00402 
00403     // check didn't evaluate any sample more than once
00404     CPPUNIT_ASSERT( logger.no_duplicate_evals() );
00405 
00406     // count total number of sample points in mesh
00407     vector< Mesh::ElementHandle > elems;
00408     myMesh.get_all_elements( elems, err );
00409     ASSERT_NO_ERROR( err );
00410 
00411     vector< Mesh::VertexHandle > free_verts;
00412     PatchData global_patch;
00413     global_patch.set_mesh( &myMesh );
00414     global_patch.set_domain( &myDomain );
00415     global_patch.attach_settings( &mySettings );
00416     global_patch.set_mesh_entities( elems, free_verts, err );
00417     ASSERT_NO_ERROR( err );
00418 
00419     size_t num_samples = 0;
00420     for( size_t i = 0; i < global_patch.num_elements(); ++i )
00421         num_samples += global_patch.get_samples( i ).num_nodes();
00422 
00423     // check correct number of metric evaluations
00424     CPPUNIT_ASSERT_EQUAL( num_samples, logger.mValues.size() );
00425 
00426     // check values
00427     const QualityAssessor::Assessor* results = qa.get_results( &logger );
00428     CPPUNIT_ASSERT( results != NULL );
00429     CPPUNIT_ASSERT_EQUAL( num_samples, (size_t)results->get_count() );
00430     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.min(), results->get_minimum(), 1e-6 );
00431     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.avg(), results->get_average(), 1e-6 );
00432     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.rms(), results->get_rms(), 1e-6 );
00433     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.max(), results->get_maximum(), 1e-6 );
00434     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.dev(), results->get_stddev(), 1e-6 );
00435 }
00436 
00437 void QualityAssessorTest::test_histogram_known_range()
00438 {
00439     const double lower = 1.0, upper = 3.0;
00440     const int intervals = 5;
00441 
00442     AspectRatioGammaQualityMetric metric;
00443     MetricLogger logger( &metric );
00444     QualityAssessor qa;
00445     qa.add_histogram_assessment( &logger, lower, upper, intervals );
00446     qa.disable_printing_results();
00447 
00448     MsqError err;
00449     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00450     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00451     ASSERT_NO_ERROR( err );
00452 
00453     // calculate expected histogram
00454     std::vector< int > counts( intervals + 2, 0 );
00455     for( vector< double >::iterator i = logger.mValues.begin(); i != logger.mValues.end(); ++i )
00456     {
00457         int bucket = (int)( ( *i - lower ) * intervals / ( upper - lower ) );
00458         if( bucket < 0 )
00459             ++counts.front();
00460         else if( bucket >= intervals )
00461             ++counts.back();
00462         else
00463             ++counts[bucket + 1];
00464     }
00465 
00466     // check values
00467     const QualityAssessor::Assessor* results = qa.get_results( &logger );
00468     CPPUNIT_ASSERT( results != NULL );
00469     CPPUNIT_ASSERT( results->have_histogram() );
00470 
00471     double r_lower, r_upper;
00472     vector< int > r_counts;
00473     results->get_histogram( r_lower, r_upper, r_counts, err );
00474     ASSERT_NO_ERROR( err );
00475     CPPUNIT_ASSERT_DOUBLES_EQUAL( lower, r_lower, DBL_EPSILON );
00476     CPPUNIT_ASSERT_DOUBLES_EQUAL( upper, r_upper, DBL_EPSILON );
00477     CPPUNIT_ASSERT_EQUAL( intervals + 2, (int)r_counts.size() );
00478 
00479     CPPUNIT_ASSERT_EQUAL( counts.front(), r_counts.front() );
00480     CPPUNIT_ASSERT_EQUAL( counts.back(), r_counts.back() );
00481     switch( intervals )
00482     {
00483         default:
00484             for( unsigned i = 11; i < intervals + 1; ++i )
00485                 CPPUNIT_ASSERT_EQUAL( counts[i], r_counts[i] );
00486         case 10:
00487             CPPUNIT_ASSERT_EQUAL( counts[10], r_counts[10] );
00488         case 9:
00489             CPPUNIT_ASSERT_EQUAL( counts[9], r_counts[9] );
00490         case 8:
00491             CPPUNIT_ASSERT_EQUAL( counts[8], r_counts[8] );
00492         case 7:
00493             CPPUNIT_ASSERT_EQUAL( counts[7], r_counts[7] );
00494         case 6:
00495             CPPUNIT_ASSERT_EQUAL( counts[6], r_counts[6] );
00496         case 5:
00497             CPPUNIT_ASSERT_EQUAL( counts[5], r_counts[5] );
00498         case 4:
00499             CPPUNIT_ASSERT_EQUAL( counts[4], r_counts[4] );
00500         case 3:
00501             CPPUNIT_ASSERT_EQUAL( counts[3], r_counts[3] );
00502         case 2:
00503             CPPUNIT_ASSERT_EQUAL( counts[2], r_counts[2] );
00504         case 1:
00505             CPPUNIT_ASSERT_EQUAL( counts[1], r_counts[1] );
00506         case 0:;
00507     }
00508 }
00509 
00510 void QualityAssessorTest::test_histogram_unknown_range()
00511 {
00512     const int intervals = 10;
00513 
00514     AspectRatioGammaQualityMetric metric;
00515     MetricLogger logger( &metric );
00516     QualityAssessor qa( &logger, intervals );
00517     qa.disable_printing_results();
00518 
00519     MsqError err;
00520     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00521     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00522     ASSERT_NO_ERROR( err );
00523 
00524     // check values
00525     const QualityAssessor::Assessor* results = qa.get_results( &logger );
00526     CPPUNIT_ASSERT( results != NULL );
00527     CPPUNIT_ASSERT( results->have_histogram() );
00528 
00529     double r_lower, r_upper;
00530     vector< int > r_counts;
00531     results->get_histogram( r_lower, r_upper, r_counts, err );
00532     ASSERT_NO_ERROR( err );
00533     CPPUNIT_ASSERT( r_lower <= logger.min() );
00534     CPPUNIT_ASSERT( r_upper >= logger.min() );
00535     // allow some freedom in choice of range, but not more than 30%
00536     CPPUNIT_ASSERT( ( 0.30 * ( r_upper - r_lower ) ) <= ( logger.max() - logger.min() ) );
00537     // range must contain all metric values
00538     CPPUNIT_ASSERT_EQUAL( 0, r_counts.front() );
00539     CPPUNIT_ASSERT_EQUAL( 0, r_counts.back() );
00540     CPPUNIT_ASSERT_EQUAL( intervals + 2, (int)r_counts.size() );
00541 
00542     // calculate expected histogram
00543     std::vector< int > counts( intervals, 0 );
00544     for( vector< double >::iterator i = logger.mValues.begin(); i != logger.mValues.end(); ++i )
00545     {
00546         double fract = ( *i - r_lower ) / ( r_upper - r_lower );
00547         int bucket;
00548         if( fabs( fract - 1.0 ) < DBL_EPSILON )
00549             bucket = intervals - 1;
00550         else
00551             bucket = (int)( intervals * fract );
00552         CPPUNIT_ASSERT( bucket >= 0 );
00553         CPPUNIT_ASSERT( bucket < intervals );
00554         ++counts[bucket];
00555     }
00556 
00557     // QA should have evaluated metric twice, so adjust values
00558     CPPUNIT_ASSERT_EQUAL( 12, (int)logger.mValues.size() );
00559     for( vector< int >::iterator j = counts.begin(); j != counts.end(); ++j )
00560         *j /= 2;
00561 
00562     switch( intervals )
00563     {
00564         default:
00565             for( int i = 10; i < intervals; ++i )
00566                 CPPUNIT_ASSERT_EQUAL( counts[i], r_counts[i + 1] );
00567         case 10:
00568             CPPUNIT_ASSERT_EQUAL( counts[9], r_counts[10] );
00569         case 9:
00570             CPPUNIT_ASSERT_EQUAL( counts[8], r_counts[9] );
00571         case 8:
00572             CPPUNIT_ASSERT_EQUAL( counts[7], r_counts[8] );
00573         case 7:
00574             CPPUNIT_ASSERT_EQUAL( counts[6], r_counts[7] );
00575         case 6:
00576             CPPUNIT_ASSERT_EQUAL( counts[5], r_counts[6] );
00577         case 5:
00578             CPPUNIT_ASSERT_EQUAL( counts[4], r_counts[5] );
00579         case 4:
00580             CPPUNIT_ASSERT_EQUAL( counts[3], r_counts[4] );
00581         case 3:
00582             CPPUNIT_ASSERT_EQUAL( counts[2], r_counts[3] );
00583         case 2:
00584             CPPUNIT_ASSERT_EQUAL( counts[1], r_counts[2] );
00585         case 1:
00586             CPPUNIT_ASSERT_EQUAL( counts[0], r_counts[1] );
00587         case 0:;
00588     }
00589 }
00590 
00591 void QualityAssessorTest::test_power_mean()
00592 {
00593     const double P1 = 3.0, P2 = 0.5;
00594     AspectRatioGammaQualityMetric metric1;
00595     LocalSizeQualityMetric metric2;
00596     ConditionNumberQualityMetric metric3;
00597     MetricLogger logger1( &metric1 ), logger2( &metric2 );
00598     QualityAssessor qa( &metric3 );
00599     qa.add_quality_assessment( &logger1, 0, P1 );
00600     qa.add_quality_assessment( &logger2, 0, P2 );
00601     qa.disable_printing_results();
00602 
00603     MsqError err;
00604     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00605     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00606     ASSERT_NO_ERROR( err );
00607 
00608     // get results
00609     const QualityAssessor::Assessor *result1, *result2, *result3;
00610     result1 = qa.get_results( &logger1 );
00611     result2 = qa.get_results( &logger2 );
00612     result3 = qa.get_results( &metric3 );
00613     CPPUNIT_ASSERT( NULL != result1 );
00614     CPPUNIT_ASSERT( NULL != result2 );
00615     CPPUNIT_ASSERT( NULL != result3 );
00616 
00617     // check flags
00618     CPPUNIT_ASSERT( result1->have_power_mean() );
00619     CPPUNIT_ASSERT( result2->have_power_mean() );
00620     CPPUNIT_ASSERT( !result3->have_power_mean() );
00621 
00622     // check values
00623     CPPUNIT_ASSERT_DOUBLES_EQUAL( P1, result1->get_power(), DBL_EPSILON );
00624     CPPUNIT_ASSERT_DOUBLES_EQUAL( P2, result2->get_power(), DBL_EPSILON );
00625     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger1.pmean( P1 ), result1->get_power_mean(), logger1.pmean( P1 ) * DBL_EPSILON );
00626     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger2.pmean( P2 ), result2->get_power_mean(), logger2.pmean( P2 ) * DBL_EPSILON );
00627 }
00628 
00629 void QualityAssessorTest::test_invalid_count()
00630 {
00631     MsqError err;
00632     ConditionNumberQualityMetric metric;
00633     QualityAssessor qa( &metric );
00634     qa.measure_free_samples_only( false );
00635     qa.disable_printing_results();
00636     const QualityAssessor::Assessor* results = qa.get_results( &metric );
00637     CPPUNIT_ASSERT( NULL != results );
00638 
00639     // try mesh with only valid elements
00640     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00641     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00642     ASSERT_NO_ERROR( err );
00643     CPPUNIT_ASSERT_EQUAL( 0, results->get_invalid_element_count() );
00644 
00645     // try mesh with one inverted element
00646     MeshDomainAssoc mesh_and_domain2 = MeshDomainAssoc( &invertedMesh, &myDomain );
00647     qa.loop_over_mesh( &mesh_and_domain2, &mySettings, err );
00648     ASSERT_NO_ERROR( err );
00649     CPPUNIT_ASSERT_EQUAL( 1, results->get_invalid_element_count() );
00650 }
00651 
00652 void QualityAssessorTest::test_inverted_count()
00653 {
00654     int inverted, samples;
00655     MsqError err;
00656     QualityAssessor qa;
00657     qa.measure_free_samples_only( false );
00658     qa.disable_printing_results();
00659 
00660     // try mesh with only valid elements
00661     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00662     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00663     ASSERT_NO_ERROR( err );
00664     inverted = samples = -1;
00665     qa.get_inverted_element_count( inverted, samples, err );
00666     ASSERT_NO_ERROR( err );
00667     CPPUNIT_ASSERT_EQUAL( 0, inverted );
00668     CPPUNIT_ASSERT_EQUAL( 0, samples );
00669 
00670     // try mesh with one inverted element
00671     MeshDomainAssoc mesh_and_domain2 = MeshDomainAssoc( &invertedMesh, &myDomain );
00672     qa.loop_over_mesh( &mesh_and_domain2, &mySettings, err );
00673     ASSERT_NO_ERROR( err );
00674     inverted = samples = -1;
00675     qa.get_inverted_element_count( inverted, samples, err );
00676     ASSERT_NO_ERROR( err );
00677     CPPUNIT_ASSERT_EQUAL( 1, inverted );
00678     CPPUNIT_ASSERT_EQUAL( 1, samples );
00679 }
00680 
00681 // Define class to handle redirecting output streams
00682 // (e.g. std::cout) to a stringstream instance.  We
00683 // use this to catch errors where QualityAssessor is
00684 // writing to a stream it should not be.
00685 class StreamRedirector
00686 {
00687   private:
00688     streambuf *outbuf, *errbuf, *logbuf;
00689 
00690   public:
00691     ostringstream outstr, errstr, logstr;
00692     StreamRedirector()
00693     {
00694         outbuf = cout.rdbuf();
00695         errbuf = cerr.rdbuf();
00696         logbuf = clog.rdbuf();
00697     }
00698 
00699     ~StreamRedirector()
00700     {
00701         restore();
00702     }
00703 
00704     void redirect()
00705     {
00706         clear();
00707         cout.rdbuf( outstr.rdbuf() );
00708         cerr.rdbuf( errstr.rdbuf() );
00709         clog.rdbuf( logstr.rdbuf() );
00710     }
00711 
00712     void restore()
00713     {
00714         cout.rdbuf( outbuf );
00715         cerr.rdbuf( errbuf );
00716         clog.rdbuf( logbuf );
00717     }
00718     // reset all streams to empty
00719     void clear()
00720     {
00721         string empty;
00722         outstr.str( empty );
00723         errstr.str( empty );
00724         logstr.str( empty );
00725     }
00726 
00727     void print_stream_bytes( const char* prefix, string bytes )
00728     {
00729         // fix output stream so we can use it to print results
00730         streambuf* sbuf = cout.rdbuf();
00731         cout.rdbuf( outbuf );
00732         cout << prefix << ": " << bytes.size() << " characters: ";
00733         for( string::iterator i = bytes.begin(); i != bytes.end(); ++i )
00734             if( isprint( *i ) || *i == ' ' )
00735                 cout << *i;
00736             else
00737                 switch( *i )
00738                 {
00739                     case '\a':
00740                         cout << "\\a";
00741                         break;
00742                     case '\b':
00743                         cout << "\\b";
00744                         break;
00745                     case '\t':
00746                         cout << "\\t";
00747                         break;
00748                     case '\n':
00749                         cout << "\\n";
00750                         break;
00751                     case '\v':
00752                         cout << "\\v";
00753                         break;
00754                     case '\f':
00755                         cout << "\\f";
00756                         break;
00757                     case '\r':
00758                         cout << "\\r";
00759                         break;
00760                     case '\0':
00761                         cout << "\\0";
00762                         break;
00763                     default:
00764                         cout << '\\' << setw( 3 ) << setfill( '0' ) << oct << *i;
00765                         break;
00766                 }
00767         cout << endl;
00768         // restore cout to incomming state
00769         cout.rdbuf( sbuf );
00770     }
00771 
00772     void check_stream( const char* name, ostringstream& str, bool& result )
00773     {
00774         string s = str.str();
00775         if( s.empty() ) return;
00776         result = true;
00777         print_stream_bytes( name, s );
00778     }
00779 
00780     // check if any stream was written to (may clear streams, depending on platform)
00781     bool have_data()
00782     {
00783         bool result = false;
00784         check_stream( "std::cout", outstr, result );
00785         check_stream( "std::cerr", errstr, result );
00786         check_stream( "std::clog", logstr, result );
00787         return result;
00788     }
00789 };
00790 
00791 void QualityAssessorTest::test_output_control()
00792 {
00793     // Redirect std::cout, std:cerr, and std::clog
00794     // so we know when they're written to.
00795     StreamRedirector redir;
00796 
00797     MsqError err;
00798     ConditionNumberQualityMetric metric;
00799 
00800     // disable output from constructor
00801     QualityAssessor qa1( &metric, 0, 0, false, 0, false );
00802     redir.redirect();
00803     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &invertedMesh, &myDomain );
00804     qa1.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00805     redir.restore();
00806     // make sure nothing was written to output streams
00807     CPPUNIT_ASSERT( !redir.have_data() );
00808 
00809     // disable output from function
00810     QualityAssessor qa2( &metric );
00811     qa2.disable_printing_results();
00812     redir.redirect();
00813     qa2.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00814     redir.restore();
00815     // make sure nothing was written to output streams
00816     CPPUNIT_ASSERT( !redir.have_data() );
00817 
00818     // send output to alternate stream
00819     stringstream deststr;
00820     QualityAssessor qa3( deststr, &metric );
00821     redir.redirect();
00822     qa3.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00823     redir.restore();
00824     // make sure nothing was written to output streams
00825     CPPUNIT_ASSERT( !redir.have_data() );
00826     // check that some data was written to specified stream
00827     CPPUNIT_ASSERT( !deststr.str().empty() );
00828 
00829     // redir destructor should now restore normal output streams
00830 }
00831 
00832 void QualityAssessorTest::test_tag_element()
00833 {
00834     const char* tag_name = "xxxTEST_TAGxxx";
00835     ConditionNumberQualityMetric metric;
00836     MetricLogger logger( &metric );
00837     QualityAssessor qa( &logger, 0, 0, false, tag_name );
00838     qa.disable_printing_results();
00839 
00840     MsqError err;
00841     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00842     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00843     ASSERT_NO_ERROR( err );
00844 
00845     TagHandle tag;
00846     tag = myMesh.tag_get( tag_name, err );
00847     ASSERT_NO_ERROR( err );
00848 
00849     string name;
00850     Mesh::TagType type;
00851     unsigned length;
00852     myMesh.tag_properties( tag, name, type, length, err );
00853     ASSERT_NO_ERROR( err );
00854     CPPUNIT_ASSERT_EQUAL( name, string( tag_name ) );
00855     CPPUNIT_ASSERT_EQUAL( Mesh::DOUBLE, type );
00856     CPPUNIT_ASSERT_EQUAL( 1u, length );
00857 
00858     CPPUNIT_ASSERT_EQUAL( logger.mValues.size(), logger.mHandles.size() );
00859     CPPUNIT_ASSERT( !logger.mValues.empty() );
00860     vector< double > tag_values( logger.mValues.size() );
00861     myMesh.tag_get_element_data( tag, logger.mHandles.size(), &logger.mHandles[0], arrptr( tag_values ), err );
00862     ASSERT_NO_ERROR( err );
00863 
00864     for( unsigned i = 0; i < logger.mValues.size(); ++i )
00865         CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.mValues[1], tag_values[1], DBL_EPSILON );
00866 }
00867 
00868 void QualityAssessorTest::test_tag_vertex()
00869 {
00870     const char* tag_name = "vertex-test-tag";
00871     LocalSizeQualityMetric metric;
00872     MetricLogger logger( &metric );
00873     QualityAssessor qa( &logger, 0, 0, false, tag_name );
00874     qa.disable_printing_results();
00875 
00876     MsqError err;
00877     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00878     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00879     ASSERT_NO_ERROR( err );
00880 
00881     TagHandle tag;
00882     tag = myMesh.tag_get( tag_name, err );
00883     ASSERT_NO_ERROR( err );
00884 
00885     string name;
00886     Mesh::TagType type;
00887     unsigned length;
00888     myMesh.tag_properties( tag, name, type, length, err );
00889     ASSERT_NO_ERROR( err );
00890     CPPUNIT_ASSERT_EQUAL( name, string( tag_name ) );
00891     CPPUNIT_ASSERT_EQUAL( Mesh::DOUBLE, type );
00892     CPPUNIT_ASSERT_EQUAL( 1u, length );
00893 
00894     CPPUNIT_ASSERT_EQUAL( logger.mValues.size(), logger.mHandles.size() );
00895     vector< double > tag_values( logger.mValues.size() );
00896     myMesh.tag_get_vertex_data( tag, logger.mHandles.size(), &logger.mHandles[0], arrptr( tag_values ), err );
00897     ASSERT_NO_ERROR( err );
00898 
00899     for( unsigned i = 0; i < logger.mValues.size(); ++i )
00900         CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.mValues[1], tag_values[1], DBL_EPSILON );
00901 }
00902 
00903 void QualityAssessorTest::test_tag_inverted()
00904 {
00905     const char* tag_name = "123INVERT456";
00906     QualityAssessor qa( false, false, tag_name );
00907 
00908     MsqError err;
00909     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &invertedMesh, &myDomain );
00910     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00911     ASSERT_NO_ERROR( err );
00912 
00913     TagHandle tag;
00914     tag = invertedMesh.tag_get( tag_name, err );
00915     ASSERT_NO_ERROR( err );
00916 
00917     string name;
00918     Mesh::TagType type;
00919     unsigned length;
00920     invertedMesh.tag_properties( tag, name, type, length, err );
00921     ASSERT_NO_ERROR( err );
00922     CPPUNIT_ASSERT_EQUAL( name, string( tag_name ) );
00923     CPPUNIT_ASSERT_EQUAL( Mesh::INT, type );
00924     CPPUNIT_ASSERT_EQUAL( 1u, length );
00925 
00926     vector< Mesh::ElementHandle > elements;
00927     invertedMesh.get_all_elements( elements, err );
00928     ASSERT_NO_ERROR( err );
00929 
00930     // expect two elements, one inverted
00931     CPPUNIT_ASSERT_EQUAL( (size_t)2, elements.size() );
00932     int data[2];
00933     invertedMesh.tag_get_element_data( tag, 2, arrptr( elements ), data, err );
00934     ASSERT_NO_ERROR( err );
00935 
00936     if( data[0] )
00937     {
00938         CPPUNIT_ASSERT_EQUAL( 0, data[1] );
00939     }
00940     else
00941     {
00942         CPPUNIT_ASSERT_EQUAL( 0, data[0] );
00943         CPPUNIT_ASSERT_EQUAL( 1, data[1] );
00944     }
00945 }
00946 
00947 void QualityAssessorTest::test_print_inverted()
00948 {
00949     MsqError err;
00950     stringstream str;
00951     QualityAssessor qa( str );
00952     qa.measure_free_samples_only( false );
00953     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &invertedMesh, &myDomain );
00954     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00955     ASSERT_NO_ERROR( err );
00956 
00957     // get inverted count from QA
00958     int expected = 0, other = 0;
00959     qa.get_inverted_element_count( expected, other, err );
00960     ASSERT_NO_ERROR( err );
00961 
00962     // At some point, expect "n INVERTED" to appear in output,
00963     // where 'n' is the count of inverted elements.
00964     int value;
00965     string prev, curr;
00966     for( ; str >> curr && curr != "INVERTED"; prev = curr )
00967         ;
00968     str.str( prev );
00969     str >> value;
00970 
00971     CPPUNIT_ASSERT_EQUAL( expected, value );
00972 }
00973 
00974 void QualityAssessorTest::test_print_stats()
00975 {
00976     MsqError err;
00977     stringstream str;
00978     ConditionNumberQualityMetric metric;
00979     QualityAssessor qa( str, &metric );
00980     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00981     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00982     ASSERT_NO_ERROR( err );
00983 
00984     // get results
00985     const QualityAssessor::Assessor* results = qa.get_results( &metric );
00986     CPPUNIT_ASSERT( NULL != results );
00987 
00988     // Expect metric name followed by 5 numerical values for statistics,
00989     // at some point in output.
00990 
00991     // Find occurance of metric name in output
00992     for( ;; )
00993     {
00994         stringstream name( metric.get_name() );
00995         string s, n;
00996         name >> n;
00997         while( str >> s && s != n )
00998             ;
00999         CPPUNIT_ASSERT( str );
01000         while( name >> n && str >> s && s == n )
01001             ;
01002         if( !name ) break;
01003     }
01004 
01005     // Read stats
01006     double min_s, max_s, avg_s, rms_s, dev_s;
01007     str >> min_s >> avg_s >> rms_s >> max_s >> dev_s;
01008 
01009     // The following commented out because they no longer pass due to a change
01010     //  in the QA Summary format
01011 
01012     // compare results
01013     //  CPPUNIT_ASSERT_DOUBLES_EQUAL( results->get_minimum(), min_s, min_s * 0.01 );
01014     //  CPPUNIT_ASSERT_DOUBLES_EQUAL( results->get_average(), avg_s, avg_s * 0.01 );
01015     //  CPPUNIT_ASSERT_DOUBLES_EQUAL( results->get_rms    (), rms_s, rms_s * 0.01 );
01016     //  CPPUNIT_ASSERT_DOUBLES_EQUAL( results->get_maximum(), max_s, max_s * 0.01 );
01017     //  CPPUNIT_ASSERT_DOUBLES_EQUAL( results->get_stddev (), dev_s, dev_s * 0.01 );
01018 }
01019 
01020 void QualityAssessorTest::test_print_name()
01021 {
01022     const char* NAME = "test_print_name_NAME";
01023 
01024     // expect QualityAssessor Name to be printed somewhere in output.
01025     MsqError err;
01026     stringstream str;
01027     ConditionNumberQualityMetric metric;
01028     QualityAssessor qa( str, &metric, 0, 0, false, 0, 0, NAME );
01029     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
01030     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
01031     ASSERT_NO_ERROR( err );
01032 
01033     // seach output for first occurance of name
01034     string s;
01035     while( str >> s && s != NAME )
01036         ;
01037     CPPUNIT_ASSERT_EQUAL( s, string( NAME ) );
01038 }
01039 
01040 // Test that adding a metric more than once changes
01041 // properties for a single occurance of the metric,
01042 // rather than introducing multiple copies.
01043 void QualityAssessorTest::test_modify_metric()
01044 {
01045     MsqError err;
01046     ConditionNumberQualityMetric metric;
01047     QualityAssessor qa( &metric );
01048     const QualityAssessor::Assessor* results = qa.get_results( &metric );
01049     CPPUNIT_ASSERT( NULL != results );
01050 
01051     CPPUNIT_ASSERT( !results->have_histogram() );
01052     qa.add_histogram_assessment( &metric, 0.0, 1.0, 10 );
01053     CPPUNIT_ASSERT( results->have_histogram() );
01054 
01055     CPPUNIT_ASSERT( !results->have_power_mean() );
01056     qa.add_quality_assessment( &metric, 0, 3.0 );
01057     CPPUNIT_ASSERT( results->have_power_mean() );
01058 }
01059 
01060 void QualityAssessorTest::test_free_only()
01061 {
01062     /*
01063      *   +--+--------o-----o
01064      *   |   \       |     |    + - fixed
01065      *   |    \      |     |    o - free
01066      *   +-----+-----o-----o
01067      */
01068 
01069     double coords[]      = { 0, 0, 0, 1, 0, 0, 0.5, 1, 0, 0, 1, 0, 2, 0, 0, 3, 0, 0, 3, 1, 0, 2, 1, 0 };
01070     int fixed[]          = { true, true, true, true, false, false, false, false };
01071     unsigned long conn[] = { 0, 1, 2, 3, 1, 4, 7, 2, 4, 5, 6, 7 };
01072 
01073     MsqError err;
01074     ArrayMesh mesh( 3, 8, coords, fixed, 3, QUADRILATERAL, conn );
01075     ConditionNumberQualityMetric metric;
01076     QualityAssessor qa_all( &metric, 0, 0.0, false, 0, false );
01077     QualityAssessor qa_free( &metric, 0, 0.0, true, 0, false );
01078     InstructionQueue q;
01079     q.add_quality_assessor( &qa_all, err );
01080     q.add_quality_assessor( &qa_free, err );
01081     PlanarDomain xy( PlanarDomain::XY );
01082     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &xy );
01083     q.run_instructions( &mesh_and_domain, err );
01084     ASSERT_NO_ERROR( err );
01085 
01086     const QualityAssessor::Assessor* data;
01087     data = qa_all.get_results( &metric );
01088     CPPUNIT_ASSERT( 0 != data );
01089     CPPUNIT_ASSERT_EQUAL( 3, data->get_count() );
01090 
01091     data = qa_free.get_results( &metric );
01092     CPPUNIT_ASSERT( 0 != data );
01093     CPPUNIT_ASSERT_EQUAL( 2, data->get_count() );
01094 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines