MOAB: Mesh Oriented datABase  (version 5.2.1)
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&, size_t, double&, vector< size_t >&, vector< Vector3D >&,
00308                                 vector< Matrix3D >&, MsqError& err )
00309     {  // shouldn't be called by QualityAssessor
00310         CPPUNIT_ASSERT( false );
00311         MSQ_SETERR( err )( MsqError::NOT_IMPLEMENTED );
00312         return false;
00313     }
00314 
00315   private:
00316     QualityMetric* mMetric;
00317 };
00318 
00319 void QualityAssessorTest::test_basic_stats_element()
00320 {
00321     AspectRatioGammaQualityMetric metric;
00322     MetricLogger logger( &metric );
00323     QualityAssessor qa( &logger );
00324     qa.disable_printing_results();
00325 
00326     MsqError err;
00327     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00328     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00329     ASSERT_NO_ERROR( err );
00330 
00331     // check didn't evaluate any element more than once
00332     CPPUNIT_ASSERT( logger.no_duplicate_evals() );
00333 
00334     // check one eval for each element
00335     vector< Mesh::ElementHandle > elems;
00336     myMesh.get_all_elements( elems, err );
00337     ASSERT_NO_ERROR( err );
00338     CPPUNIT_ASSERT_EQUAL( elems.size(), logger.mValues.size() );
00339 
00340     // check values
00341     const QualityAssessor::Assessor* results = qa.get_results( &logger );
00342     CPPUNIT_ASSERT( results != NULL );
00343     CPPUNIT_ASSERT_EQUAL( (int)elems.size(), results->get_count() );
00344     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.min(), results->get_minimum(), 1e-6 );
00345     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.avg(), results->get_average(), 1e-6 );
00346     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.rms(), results->get_rms(), 1e-6 );
00347     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.max(), results->get_maximum(), 1e-6 );
00348     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.dev(), results->get_stddev(), 1e-6 );
00349 }
00350 
00351 void QualityAssessorTest::test_basic_stats_vertex()
00352 {
00353     LocalSizeQualityMetric metric;
00354     MetricLogger logger( &metric );
00355     QualityAssessor qa( &logger );
00356     qa.measure_free_samples_only( false );
00357     qa.disable_printing_results();
00358 
00359     MsqError err;
00360     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00361     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00362     ASSERT_NO_ERROR( err );
00363 
00364     // check didn't evaluate any vertex more than once
00365     CPPUNIT_ASSERT( logger.no_duplicate_evals() );
00366 
00367     // check one eval for each element
00368     vector< Mesh::VertexHandle > verts;
00369     myMesh.get_all_vertices( verts, err );
00370     ASSERT_NO_ERROR( err );
00371     CPPUNIT_ASSERT_EQUAL( verts.size(), logger.mValues.size() );
00372 
00373     // check values
00374     const QualityAssessor::Assessor* results = qa.get_results( &logger );
00375     CPPUNIT_ASSERT( results != NULL );
00376     CPPUNIT_ASSERT_EQUAL( (int)verts.size(), results->get_count() );
00377     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.min(), results->get_minimum(), 1e-6 );
00378     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.avg(), results->get_average(), 1e-6 );
00379     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.rms(), results->get_rms(), 1e-6 );
00380     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.max(), results->get_maximum(), 1e-6 );
00381     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.dev(), results->get_stddev(), 1e-6 );
00382 }
00383 
00384 void QualityAssessorTest::test_basic_stats_sample()
00385 {
00386     TShapeNB1 tm;
00387     IdealShapeTarget tc;
00388     TQualityMetric metric( &tc, &tm );
00389     MetricLogger logger( &metric );
00390     QualityAssessor qa( &logger );
00391     qa.disable_printing_results();
00392 
00393     MsqError err;
00394     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00395     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00396     ASSERT_NO_ERROR( err );
00397 
00398     // check didn't evaluate any sample more than once
00399     CPPUNIT_ASSERT( logger.no_duplicate_evals() );
00400 
00401     // count total number of sample points in mesh
00402     vector< Mesh::ElementHandle > elems;
00403     myMesh.get_all_elements( elems, err );
00404     ASSERT_NO_ERROR( err );
00405 
00406     vector< Mesh::VertexHandle > free_verts;
00407     PatchData global_patch;
00408     global_patch.set_mesh( &myMesh );
00409     global_patch.set_domain( &myDomain );
00410     global_patch.attach_settings( &mySettings );
00411     global_patch.set_mesh_entities( elems, free_verts, err );
00412     ASSERT_NO_ERROR( err );
00413 
00414     size_t num_samples = 0;
00415     for( size_t i = 0; i < global_patch.num_elements(); ++i )
00416         num_samples += global_patch.get_samples( i ).num_nodes();
00417 
00418     // check correct number of metric evaluations
00419     CPPUNIT_ASSERT_EQUAL( num_samples, logger.mValues.size() );
00420 
00421     // check values
00422     const QualityAssessor::Assessor* results = qa.get_results( &logger );
00423     CPPUNIT_ASSERT( results != NULL );
00424     CPPUNIT_ASSERT_EQUAL( num_samples, (size_t)results->get_count() );
00425     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.min(), results->get_minimum(), 1e-6 );
00426     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.avg(), results->get_average(), 1e-6 );
00427     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.rms(), results->get_rms(), 1e-6 );
00428     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.max(), results->get_maximum(), 1e-6 );
00429     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.dev(), results->get_stddev(), 1e-6 );
00430 }
00431 
00432 void QualityAssessorTest::test_histogram_known_range()
00433 {
00434     const double lower = 1.0, upper = 3.0;
00435     const int intervals = 5;
00436 
00437     AspectRatioGammaQualityMetric metric;
00438     MetricLogger logger( &metric );
00439     QualityAssessor qa;
00440     qa.add_histogram_assessment( &logger, lower, upper, intervals );
00441     qa.disable_printing_results();
00442 
00443     MsqError err;
00444     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00445     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00446     ASSERT_NO_ERROR( err );
00447 
00448     // calculate expected histogram
00449     std::vector< int > counts( intervals + 2, 0 );
00450     for( vector< double >::iterator i = logger.mValues.begin(); i != logger.mValues.end(); ++i )
00451     {
00452         int bucket = (int)( ( *i - lower ) * intervals / ( upper - lower ) );
00453         if( bucket < 0 )
00454             ++counts.front();
00455         else if( bucket >= intervals )
00456             ++counts.back();
00457         else
00458             ++counts[bucket + 1];
00459     }
00460 
00461     // check values
00462     const QualityAssessor::Assessor* results = qa.get_results( &logger );
00463     CPPUNIT_ASSERT( results != NULL );
00464     CPPUNIT_ASSERT( results->have_histogram() );
00465 
00466     double r_lower, r_upper;
00467     vector< int > r_counts;
00468     results->get_histogram( r_lower, r_upper, r_counts, err );
00469     ASSERT_NO_ERROR( err );
00470     CPPUNIT_ASSERT_DOUBLES_EQUAL( lower, r_lower, DBL_EPSILON );
00471     CPPUNIT_ASSERT_DOUBLES_EQUAL( upper, r_upper, DBL_EPSILON );
00472     CPPUNIT_ASSERT_EQUAL( intervals + 2, (int)r_counts.size() );
00473 
00474     CPPUNIT_ASSERT_EQUAL( counts.front(), r_counts.front() );
00475     CPPUNIT_ASSERT_EQUAL( counts.back(), r_counts.back() );
00476     switch( intervals )
00477     {
00478         default:
00479             for( unsigned i = 11; i < intervals + 1; ++i )
00480                 CPPUNIT_ASSERT_EQUAL( counts[i], r_counts[i] );
00481         case 10:
00482             CPPUNIT_ASSERT_EQUAL( counts[10], r_counts[10] );
00483         case 9:
00484             CPPUNIT_ASSERT_EQUAL( counts[9], r_counts[9] );
00485         case 8:
00486             CPPUNIT_ASSERT_EQUAL( counts[8], r_counts[8] );
00487         case 7:
00488             CPPUNIT_ASSERT_EQUAL( counts[7], r_counts[7] );
00489         case 6:
00490             CPPUNIT_ASSERT_EQUAL( counts[6], r_counts[6] );
00491         case 5:
00492             CPPUNIT_ASSERT_EQUAL( counts[5], r_counts[5] );
00493         case 4:
00494             CPPUNIT_ASSERT_EQUAL( counts[4], r_counts[4] );
00495         case 3:
00496             CPPUNIT_ASSERT_EQUAL( counts[3], r_counts[3] );
00497         case 2:
00498             CPPUNIT_ASSERT_EQUAL( counts[2], r_counts[2] );
00499         case 1:
00500             CPPUNIT_ASSERT_EQUAL( counts[1], r_counts[1] );
00501         case 0:;
00502     }
00503 }
00504 
00505 void QualityAssessorTest::test_histogram_unknown_range()
00506 {
00507     const int intervals = 10;
00508 
00509     AspectRatioGammaQualityMetric metric;
00510     MetricLogger logger( &metric );
00511     QualityAssessor qa( &logger, intervals );
00512     qa.disable_printing_results();
00513 
00514     MsqError err;
00515     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00516     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00517     ASSERT_NO_ERROR( err );
00518 
00519     // check values
00520     const QualityAssessor::Assessor* results = qa.get_results( &logger );
00521     CPPUNIT_ASSERT( results != NULL );
00522     CPPUNIT_ASSERT( results->have_histogram() );
00523 
00524     double r_lower, r_upper;
00525     vector< int > r_counts;
00526     results->get_histogram( r_lower, r_upper, r_counts, err );
00527     ASSERT_NO_ERROR( err );
00528     CPPUNIT_ASSERT( r_lower <= logger.min() );
00529     CPPUNIT_ASSERT( r_upper >= logger.min() );
00530     // allow some freedom in choice of range, but not more than 30%
00531     CPPUNIT_ASSERT( ( 0.30 * ( r_upper - r_lower ) ) <= ( logger.max() - logger.min() ) );
00532     // range must contain all metric values
00533     CPPUNIT_ASSERT_EQUAL( 0, r_counts.front() );
00534     CPPUNIT_ASSERT_EQUAL( 0, r_counts.back() );
00535     CPPUNIT_ASSERT_EQUAL( intervals + 2, (int)r_counts.size() );
00536 
00537     // calculate expected histogram
00538     std::vector< int > counts( intervals, 0 );
00539     for( vector< double >::iterator i = logger.mValues.begin(); i != logger.mValues.end(); ++i )
00540     {
00541         double fract = ( *i - r_lower ) / ( r_upper - r_lower );
00542         int bucket;
00543         if( fabs( fract - 1.0 ) < DBL_EPSILON )
00544             bucket = intervals - 1;
00545         else
00546             bucket = (int)( intervals * fract );
00547         CPPUNIT_ASSERT( bucket >= 0 );
00548         CPPUNIT_ASSERT( bucket < intervals );
00549         ++counts[bucket];
00550     }
00551 
00552     // QA should have evaluated metric twice, so adjust values
00553     CPPUNIT_ASSERT_EQUAL( 12, (int)logger.mValues.size() );
00554     for( vector< int >::iterator j = counts.begin(); j != counts.end(); ++j )
00555         *j /= 2;
00556 
00557     switch( intervals )
00558     {
00559         default:
00560             for( int i = 10; i < intervals; ++i )
00561                 CPPUNIT_ASSERT_EQUAL( counts[i], r_counts[i + 1] );
00562         case 10:
00563             CPPUNIT_ASSERT_EQUAL( counts[9], r_counts[10] );
00564         case 9:
00565             CPPUNIT_ASSERT_EQUAL( counts[8], r_counts[9] );
00566         case 8:
00567             CPPUNIT_ASSERT_EQUAL( counts[7], r_counts[8] );
00568         case 7:
00569             CPPUNIT_ASSERT_EQUAL( counts[6], r_counts[7] );
00570         case 6:
00571             CPPUNIT_ASSERT_EQUAL( counts[5], r_counts[6] );
00572         case 5:
00573             CPPUNIT_ASSERT_EQUAL( counts[4], r_counts[5] );
00574         case 4:
00575             CPPUNIT_ASSERT_EQUAL( counts[3], r_counts[4] );
00576         case 3:
00577             CPPUNIT_ASSERT_EQUAL( counts[2], r_counts[3] );
00578         case 2:
00579             CPPUNIT_ASSERT_EQUAL( counts[1], r_counts[2] );
00580         case 1:
00581             CPPUNIT_ASSERT_EQUAL( counts[0], r_counts[1] );
00582         case 0:;
00583     }
00584 }
00585 
00586 void QualityAssessorTest::test_power_mean()
00587 {
00588     const double P1 = 3.0, P2 = 0.5;
00589     AspectRatioGammaQualityMetric metric1;
00590     LocalSizeQualityMetric metric2;
00591     ConditionNumberQualityMetric metric3;
00592     MetricLogger logger1( &metric1 ), logger2( &metric2 );
00593     QualityAssessor qa( &metric3 );
00594     qa.add_quality_assessment( &logger1, 0, P1 );
00595     qa.add_quality_assessment( &logger2, 0, P2 );
00596     qa.disable_printing_results();
00597 
00598     MsqError err;
00599     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00600     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00601     ASSERT_NO_ERROR( err );
00602 
00603     // get results
00604     const QualityAssessor::Assessor *result1, *result2, *result3;
00605     result1 = qa.get_results( &logger1 );
00606     result2 = qa.get_results( &logger2 );
00607     result3 = qa.get_results( &metric3 );
00608     CPPUNIT_ASSERT( NULL != result1 );
00609     CPPUNIT_ASSERT( NULL != result2 );
00610     CPPUNIT_ASSERT( NULL != result3 );
00611 
00612     // check flags
00613     CPPUNIT_ASSERT( result1->have_power_mean() );
00614     CPPUNIT_ASSERT( result2->have_power_mean() );
00615     CPPUNIT_ASSERT( !result3->have_power_mean() );
00616 
00617     // check values
00618     CPPUNIT_ASSERT_DOUBLES_EQUAL( P1, result1->get_power(), DBL_EPSILON );
00619     CPPUNIT_ASSERT_DOUBLES_EQUAL( P2, result2->get_power(), DBL_EPSILON );
00620     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger1.pmean( P1 ), result1->get_power_mean(), logger1.pmean( P1 ) * DBL_EPSILON );
00621     CPPUNIT_ASSERT_DOUBLES_EQUAL( logger2.pmean( P2 ), result2->get_power_mean(), logger2.pmean( P2 ) * DBL_EPSILON );
00622 }
00623 
00624 void QualityAssessorTest::test_invalid_count()
00625 {
00626     MsqError err;
00627     ConditionNumberQualityMetric metric;
00628     QualityAssessor qa( &metric );
00629     qa.measure_free_samples_only( false );
00630     qa.disable_printing_results();
00631     const QualityAssessor::Assessor* results = qa.get_results( &metric );
00632     CPPUNIT_ASSERT( NULL != results );
00633 
00634     // try mesh with only valid elements
00635     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00636     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00637     ASSERT_NO_ERROR( err );
00638     CPPUNIT_ASSERT_EQUAL( 0, results->get_invalid_element_count() );
00639 
00640     // try mesh with one inverted element
00641     MeshDomainAssoc mesh_and_domain2 = MeshDomainAssoc( &invertedMesh, &myDomain );
00642     qa.loop_over_mesh( &mesh_and_domain2, &mySettings, err );
00643     ASSERT_NO_ERROR( err );
00644     CPPUNIT_ASSERT_EQUAL( 1, results->get_invalid_element_count() );
00645 }
00646 
00647 void QualityAssessorTest::test_inverted_count()
00648 {
00649     int inverted, samples;
00650     MsqError err;
00651     QualityAssessor qa;
00652     qa.measure_free_samples_only( false );
00653     qa.disable_printing_results();
00654 
00655     // try mesh with only valid elements
00656     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00657     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00658     ASSERT_NO_ERROR( err );
00659     inverted = samples = -1;
00660     qa.get_inverted_element_count( inverted, samples, err );
00661     ASSERT_NO_ERROR( err );
00662     CPPUNIT_ASSERT_EQUAL( 0, inverted );
00663     CPPUNIT_ASSERT_EQUAL( 0, samples );
00664 
00665     // try mesh with one inverted element
00666     MeshDomainAssoc mesh_and_domain2 = MeshDomainAssoc( &invertedMesh, &myDomain );
00667     qa.loop_over_mesh( &mesh_and_domain2, &mySettings, err );
00668     ASSERT_NO_ERROR( err );
00669     inverted = samples = -1;
00670     qa.get_inverted_element_count( inverted, samples, err );
00671     ASSERT_NO_ERROR( err );
00672     CPPUNIT_ASSERT_EQUAL( 1, inverted );
00673     CPPUNIT_ASSERT_EQUAL( 1, samples );
00674 }
00675 
00676 // Define class to handle redirecting output streams
00677 // (e.g. std::cout) to a stringstream instance.  We
00678 // use this to catch errors where QualityAssessor is
00679 // writing to a stream it should not be.
00680 class StreamRedirector
00681 {
00682   private:
00683     streambuf *outbuf, *errbuf, *logbuf;
00684 
00685   public:
00686     ostringstream outstr, errstr, logstr;
00687     StreamRedirector()
00688     {
00689         outbuf = cout.rdbuf();
00690         errbuf = cerr.rdbuf();
00691         logbuf = clog.rdbuf();
00692     }
00693 
00694     ~StreamRedirector()
00695     {
00696         restore();
00697     }
00698 
00699     void redirect()
00700     {
00701         clear();
00702         cout.rdbuf( outstr.rdbuf() );
00703         cerr.rdbuf( errstr.rdbuf() );
00704         clog.rdbuf( logstr.rdbuf() );
00705     }
00706 
00707     void restore()
00708     {
00709         cout.rdbuf( outbuf );
00710         cerr.rdbuf( errbuf );
00711         clog.rdbuf( logbuf );
00712     }
00713     // reset all streams to empty
00714     void clear()
00715     {
00716         string empty;
00717         outstr.str( empty );
00718         errstr.str( empty );
00719         logstr.str( empty );
00720     }
00721 
00722     void print_stream_bytes( const char* prefix, string bytes )
00723     {
00724         // fix output stream so we can use it to print results
00725         streambuf* sbuf = cout.rdbuf();
00726         cout.rdbuf( outbuf );
00727         cout << prefix << ": " << bytes.size() << " characters: ";
00728         for( string::iterator i = bytes.begin(); i != bytes.end(); ++i )
00729             if( isprint( *i ) || *i == ' ' )
00730                 cout << *i;
00731             else
00732                 switch( *i )
00733                 {
00734                     case '\a':
00735                         cout << "\\a";
00736                         break;
00737                     case '\b':
00738                         cout << "\\b";
00739                         break;
00740                     case '\t':
00741                         cout << "\\t";
00742                         break;
00743                     case '\n':
00744                         cout << "\\n";
00745                         break;
00746                     case '\v':
00747                         cout << "\\v";
00748                         break;
00749                     case '\f':
00750                         cout << "\\f";
00751                         break;
00752                     case '\r':
00753                         cout << "\\r";
00754                         break;
00755                     case '\0':
00756                         cout << "\\0";
00757                         break;
00758                     default:
00759                         cout << '\\' << setw( 3 ) << setfill( '0' ) << oct << *i;
00760                         break;
00761                 }
00762         cout << endl;
00763         // restore cout to incomming state
00764         cout.rdbuf( sbuf );
00765     }
00766 
00767     void check_stream( const char* name, ostringstream& str, bool& result )
00768     {
00769         string s = str.str();
00770         if( s.empty() ) return;
00771         result = true;
00772         print_stream_bytes( name, s );
00773     }
00774 
00775     // check if any stream was written to (may clear streams, depending on platform)
00776     bool have_data()
00777     {
00778         bool result = false;
00779         check_stream( "std::cout", outstr, result );
00780         check_stream( "std::cerr", errstr, result );
00781         check_stream( "std::clog", logstr, result );
00782         return result;
00783     }
00784 };
00785 
00786 void QualityAssessorTest::test_output_control()
00787 {
00788     // Redirect std::cout, std:cerr, and std::clog
00789     // so we know when they're written to.
00790     StreamRedirector redir;
00791 
00792     MsqError err;
00793     ConditionNumberQualityMetric metric;
00794 
00795     // disable output from constructor
00796     QualityAssessor qa1( &metric, 0, 0, false, 0, false );
00797     redir.redirect();
00798     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &invertedMesh, &myDomain );
00799     qa1.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00800     redir.restore();
00801     // make sure nothing was written to output streams
00802     CPPUNIT_ASSERT( !redir.have_data() );
00803 
00804     // disable output from function
00805     QualityAssessor qa2( &metric );
00806     qa2.disable_printing_results();
00807     redir.redirect();
00808     qa2.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00809     redir.restore();
00810     // make sure nothing was written to output streams
00811     CPPUNIT_ASSERT( !redir.have_data() );
00812 
00813     // send output to alternate stream
00814     stringstream deststr;
00815     QualityAssessor qa3( deststr, &metric );
00816     redir.redirect();
00817     qa3.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00818     redir.restore();
00819     // make sure nothing was written to output streams
00820     CPPUNIT_ASSERT( !redir.have_data() );
00821     // check that some data was written to specified stream
00822     CPPUNIT_ASSERT( !deststr.str().empty() );
00823 
00824     // redir destructor should now restore normal output streams
00825 }
00826 
00827 void QualityAssessorTest::test_tag_element()
00828 {
00829     const char* tag_name = "xxxTEST_TAGxxx";
00830     ConditionNumberQualityMetric metric;
00831     MetricLogger logger( &metric );
00832     QualityAssessor qa( &logger, 0, 0, false, tag_name );
00833     qa.disable_printing_results();
00834 
00835     MsqError err;
00836     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00837     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00838     ASSERT_NO_ERROR( err );
00839 
00840     TagHandle tag;
00841     tag = myMesh.tag_get( tag_name, err );
00842     ASSERT_NO_ERROR( err );
00843 
00844     string name;
00845     Mesh::TagType type;
00846     unsigned length;
00847     myMesh.tag_properties( tag, name, type, length, err );
00848     ASSERT_NO_ERROR( err );
00849     CPPUNIT_ASSERT_EQUAL( name, string( tag_name ) );
00850     CPPUNIT_ASSERT_EQUAL( Mesh::DOUBLE, type );
00851     CPPUNIT_ASSERT_EQUAL( 1u, length );
00852 
00853     CPPUNIT_ASSERT_EQUAL( logger.mValues.size(), logger.mHandles.size() );
00854     CPPUNIT_ASSERT( !logger.mValues.empty() );
00855     vector< double > tag_values( logger.mValues.size() );
00856     myMesh.tag_get_element_data( tag, logger.mHandles.size(), &logger.mHandles[0], arrptr( tag_values ), err );
00857     ASSERT_NO_ERROR( err );
00858 
00859     for( unsigned i = 0; i < logger.mValues.size(); ++i )
00860         CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.mValues[1], tag_values[1], DBL_EPSILON );
00861 }
00862 
00863 void QualityAssessorTest::test_tag_vertex()
00864 {
00865     const char* tag_name = "vertex-test-tag";
00866     LocalSizeQualityMetric metric;
00867     MetricLogger logger( &metric );
00868     QualityAssessor qa( &logger, 0, 0, false, tag_name );
00869     qa.disable_printing_results();
00870 
00871     MsqError err;
00872     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00873     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00874     ASSERT_NO_ERROR( err );
00875 
00876     TagHandle tag;
00877     tag = myMesh.tag_get( tag_name, err );
00878     ASSERT_NO_ERROR( err );
00879 
00880     string name;
00881     Mesh::TagType type;
00882     unsigned length;
00883     myMesh.tag_properties( tag, name, type, length, err );
00884     ASSERT_NO_ERROR( err );
00885     CPPUNIT_ASSERT_EQUAL( name, string( tag_name ) );
00886     CPPUNIT_ASSERT_EQUAL( Mesh::DOUBLE, type );
00887     CPPUNIT_ASSERT_EQUAL( 1u, length );
00888 
00889     CPPUNIT_ASSERT_EQUAL( logger.mValues.size(), logger.mHandles.size() );
00890     vector< double > tag_values( logger.mValues.size() );
00891     myMesh.tag_get_vertex_data( tag, logger.mHandles.size(), &logger.mHandles[0], arrptr( tag_values ), err );
00892     ASSERT_NO_ERROR( err );
00893 
00894     for( unsigned i = 0; i < logger.mValues.size(); ++i )
00895         CPPUNIT_ASSERT_DOUBLES_EQUAL( logger.mValues[1], tag_values[1], DBL_EPSILON );
00896 }
00897 
00898 void QualityAssessorTest::test_tag_inverted()
00899 {
00900     const char* tag_name = "123INVERT456";
00901     QualityAssessor qa( false, false, tag_name );
00902 
00903     MsqError err;
00904     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &invertedMesh, &myDomain );
00905     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00906     ASSERT_NO_ERROR( err );
00907 
00908     TagHandle tag;
00909     tag = invertedMesh.tag_get( tag_name, err );
00910     ASSERT_NO_ERROR( err );
00911 
00912     string name;
00913     Mesh::TagType type;
00914     unsigned length;
00915     invertedMesh.tag_properties( tag, name, type, length, err );
00916     ASSERT_NO_ERROR( err );
00917     CPPUNIT_ASSERT_EQUAL( name, string( tag_name ) );
00918     CPPUNIT_ASSERT_EQUAL( Mesh::INT, type );
00919     CPPUNIT_ASSERT_EQUAL( 1u, length );
00920 
00921     vector< Mesh::ElementHandle > elements;
00922     invertedMesh.get_all_elements( elements, err );
00923     ASSERT_NO_ERROR( err );
00924 
00925     // expect two elements, one inverted
00926     CPPUNIT_ASSERT_EQUAL( (size_t)2, elements.size() );
00927     int data[2];
00928     invertedMesh.tag_get_element_data( tag, 2, arrptr( elements ), data, err );
00929     ASSERT_NO_ERROR( err );
00930 
00931     if( data[0] ) { CPPUNIT_ASSERT_EQUAL( 0, data[1] ); }
00932     else
00933     {
00934         CPPUNIT_ASSERT_EQUAL( 0, data[0] );
00935         CPPUNIT_ASSERT_EQUAL( 1, data[1] );
00936     }
00937 }
00938 
00939 void QualityAssessorTest::test_print_inverted()
00940 {
00941     MsqError err;
00942     stringstream str;
00943     QualityAssessor qa( str );
00944     qa.measure_free_samples_only( false );
00945     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &invertedMesh, &myDomain );
00946     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00947     ASSERT_NO_ERROR( err );
00948 
00949     // get inverted count from QA
00950     int expected = 0, other = 0;
00951     qa.get_inverted_element_count( expected, other, err );
00952     ASSERT_NO_ERROR( err );
00953 
00954     // At some point, expect "n INVERTED" to appear in output,
00955     // where 'n' is the count of inverted elements.
00956     int value;
00957     string prev, curr;
00958     for( ; str >> curr && curr != "INVERTED"; prev = curr )
00959         ;
00960     str.str( prev );
00961     str >> value;
00962 
00963     CPPUNIT_ASSERT_EQUAL( expected, value );
00964 }
00965 
00966 void QualityAssessorTest::test_print_stats()
00967 {
00968     MsqError err;
00969     stringstream str;
00970     ConditionNumberQualityMetric metric;
00971     QualityAssessor qa( str, &metric );
00972     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
00973     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
00974     ASSERT_NO_ERROR( err );
00975 
00976     // get results
00977     const QualityAssessor::Assessor* results = qa.get_results( &metric );
00978     CPPUNIT_ASSERT( NULL != results );
00979 
00980     // Expect metric name followed by 5 numerical values for statistics,
00981     // at some point in output.
00982 
00983     // Find occurance of metric name in output
00984     for( ;; )
00985     {
00986         stringstream name( metric.get_name() );
00987         string s, n;
00988         name >> n;
00989         while( str >> s && s != n )
00990             ;
00991         CPPUNIT_ASSERT( str );
00992         while( name >> n && str >> s && s == n )
00993             ;
00994         if( !name ) break;
00995     }
00996 
00997     // Read stats
00998     double min_s, max_s, avg_s, rms_s, dev_s;
00999     str >> min_s >> avg_s >> rms_s >> max_s >> dev_s;
01000 
01001     // The following commented out because they no longer pass due to a change
01002     //  in the QA Summary format
01003 
01004     // compare results
01005     //  CPPUNIT_ASSERT_DOUBLES_EQUAL( results->get_minimum(), min_s, min_s * 0.01 );
01006     //  CPPUNIT_ASSERT_DOUBLES_EQUAL( results->get_average(), avg_s, avg_s * 0.01 );
01007     //  CPPUNIT_ASSERT_DOUBLES_EQUAL( results->get_rms    (), rms_s, rms_s * 0.01 );
01008     //  CPPUNIT_ASSERT_DOUBLES_EQUAL( results->get_maximum(), max_s, max_s * 0.01 );
01009     //  CPPUNIT_ASSERT_DOUBLES_EQUAL( results->get_stddev (), dev_s, dev_s * 0.01 );
01010 }
01011 
01012 void QualityAssessorTest::test_print_name()
01013 {
01014     const char* NAME = "test_print_name_NAME";
01015 
01016     // expect QualityAssessor Name to be printed somewhere in output.
01017     MsqError err;
01018     stringstream str;
01019     ConditionNumberQualityMetric metric;
01020     QualityAssessor qa( str, &metric, 0, 0, false, 0, 0, NAME );
01021     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &myMesh, &myDomain );
01022     qa.loop_over_mesh( &mesh_and_domain, &mySettings, err );
01023     ASSERT_NO_ERROR( err );
01024 
01025     // seach output for first occurance of name
01026     string s;
01027     while( str >> s && s != NAME )
01028         ;
01029     CPPUNIT_ASSERT_EQUAL( s, string( NAME ) );
01030 }
01031 
01032 // Test that adding a metric more than once changes
01033 // properties for a single occurance of the metric,
01034 // rather than introducing multiple copies.
01035 void QualityAssessorTest::test_modify_metric()
01036 {
01037     MsqError err;
01038     ConditionNumberQualityMetric metric;
01039     QualityAssessor qa( &metric );
01040     const QualityAssessor::Assessor* results = qa.get_results( &metric );
01041     CPPUNIT_ASSERT( NULL != results );
01042 
01043     CPPUNIT_ASSERT( !results->have_histogram() );
01044     qa.add_histogram_assessment( &metric, 0.0, 1.0, 10 );
01045     CPPUNIT_ASSERT( results->have_histogram() );
01046 
01047     CPPUNIT_ASSERT( !results->have_power_mean() );
01048     qa.add_quality_assessment( &metric, 0, 3.0 );
01049     CPPUNIT_ASSERT( results->have_power_mean() );
01050 }
01051 
01052 void QualityAssessorTest::test_free_only()
01053 {
01054     /*
01055      *   +--+--------o-----o
01056      *   |   \       |     |    + - fixed
01057      *   |    \      |     |    o - free
01058      *   +-----+-----o-----o
01059      */
01060 
01061     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 };
01062     int fixed[]          = { true, true, true, true, false, false, false, false };
01063     unsigned long conn[] = { 0, 1, 2, 3, 1, 4, 7, 2, 4, 5, 6, 7 };
01064 
01065     MsqError err;
01066     ArrayMesh mesh( 3, 8, coords, fixed, 3, QUADRILATERAL, conn );
01067     ConditionNumberQualityMetric metric;
01068     QualityAssessor qa_all( &metric, 0, 0.0, false, 0, false );
01069     QualityAssessor qa_free( &metric, 0, 0.0, true, 0, false );
01070     InstructionQueue q;
01071     q.add_quality_assessor( &qa_all, err );
01072     q.add_quality_assessor( &qa_free, err );
01073     PlanarDomain xy( PlanarDomain::XY );
01074     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &xy );
01075     q.run_instructions( &mesh_and_domain, err );
01076     ASSERT_NO_ERROR( err );
01077 
01078     const QualityAssessor::Assessor* data;
01079     data = qa_all.get_results( &metric );
01080     CPPUNIT_ASSERT( 0 != data );
01081     CPPUNIT_ASSERT_EQUAL( 3, data->get_count() );
01082 
01083     data = qa_free.get_results( &metric );
01084     CPPUNIT_ASSERT( 0 != data );
01085     CPPUNIT_ASSERT_EQUAL( 2, data->get_count() );
01086 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines