MOAB: Mesh Oriented datABase  (version 5.4.1)
QualityMetricTester.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 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     (2006) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 #include "QualityMetricTester.hpp"
00028 #include "PatchData.hpp"
00029 #include "QualityMetric.hpp"
00030 #include "IdealElements.hpp"
00031 #include "UnitUtil.hpp"
00032 #include "ElemSampleQM.hpp"
00033 #include "TopologyInfo.hpp"
00034 #include "EdgeQM.hpp"
00035 
00036 #include <cppunit/extensions/HelperMacros.h>
00037 
00038 static std::vector< EntityTopology > types_in_group( QualityMetricTester::ElemTypeGroup group )
00039 {
00040     std::vector< EntityTopology > types;
00041     switch( group )
00042     {
00043         default:
00044         case QualityMetricTester::ALL:
00045             types.push_back( POLYHEDRON );
00046             types.push_back( POLYGON );
00047         case QualityMetricTester::ALL_FE:
00048             types.push_back( SEPTAHEDRON );
00049         case QualityMetricTester::ALL_FE_EXCEPT_SEPTAHEDRON:
00050             types.push_back( PYRAMID );
00051             types.push_back( PRISM );
00052         case QualityMetricTester::NON_MIXED_FE:
00053             types.push_back( HEXAHEDRON );
00054             types.push_back( QUADRILATERAL );
00055         case QualityMetricTester::SIMPLICIES:
00056             types.push_back( TETRAHEDRON );
00057             types.push_back( TRIANGLE );
00058             break;
00059 
00060         case QualityMetricTester::THREE_D:
00061             types.push_back( POLYHEDRON );
00062         case QualityMetricTester::THREE_D_FE:
00063             types.push_back( SEPTAHEDRON );
00064         case QualityMetricTester::THREE_D_FE_EXCEPT_SEPTAHEDRON:
00065             types.push_back( PYRAMID );
00066             types.push_back( PRISM );
00067         case QualityMetricTester::THREE_D_NON_MIXED_FE:
00068             types.push_back( HEXAHEDRON );
00069             types.push_back( TETRAHEDRON );
00070             break;
00071 
00072         case QualityMetricTester::TWO_D:
00073             types.push_back( POLYGON );
00074         case QualityMetricTester::TWO_D_FE:
00075             types.push_back( TRIANGLE );
00076             types.push_back( QUADRILATERAL );
00077             break;
00078     }
00079     std::reverse( types.begin(), types.end() );
00080     return types;
00081 }
00082 
00083 QualityMetricTester::QualityMetricTester( const EntityTopology* supported_elem_types,
00084                                           size_t len,
00085                                           const Settings* settings )
00086     : degenHexPyramid( false ), types( len ), mSettings( settings ),
00087       geomPlane( Vector3D( 0, 0, 1 ), Vector3D( 0, 0, 0 ) )
00088 {
00089     std::copy( supported_elem_types, supported_elem_types + len, types.begin() );
00090 }
00091 
00092 QualityMetricTester::QualityMetricTester( ElemTypeGroup group, const Settings* settings )
00093     : degenHexPyramid( false ), types( types_in_group( group ) ), mSettings( settings ),
00094       geomPlane( Vector3D( 0, 0, 1 ), Vector3D( 0, 0, 0 ) )
00095 {
00096 }
00097 
00098 void QualityMetricTester::get_ideal_tris( PatchData& pd, bool unit_area )
00099 {
00100     //      6 ------- 5     .
00101     //       /\     /\      .
00102     //      /  \   /  \     .
00103     //     /    \ /    \    .
00104     //  1 <------X------> 4 .
00105     //     \    /0\    /    .
00106     //      \  /   \  /     .
00107     //       \/     \/      .
00108     //      2 ------- 3     .
00109     static const double y3         = MSQ_SQRT_THREE_DIV_TWO;
00110     double ideal_tri_verts[]       = { 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, -0.5, -y3,  0.0, 0.5, -y3,
00111                                  0.0, 1.0, 0.0, 0.0,  0.5, y3,  0.0,  -0.5, y3,  0.0 };
00112     const size_t ideal_tri_elems[] = { 0, 1, 2, 0, 2, 3, 0, 4, 5, 0, 4, 5, 0, 5, 6, 0, 6, 1 };
00113     const bool fixed[]             = { false, true, true, true, true, true, true };
00114 
00115     if( unit_area )
00116     {
00117         const double unit_tri_scale = 2.0 / std::sqrt( 6.0 );
00118         for( size_t i = 0; i < sizeof( ideal_tri_verts ) / sizeof( double ); ++i )
00119             ideal_tri_verts[i] *= unit_tri_scale;
00120     }
00121 
00122     pd.set_domain( &geomPlane );
00123 
00124     MsqPrintError err( std::cout );
00125     pd.fill( 7, ideal_tri_verts, 6, TRIANGLE, ideal_tri_elems, fixed, err );
00126     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00127 }
00128 
00129 void QualityMetricTester::get_ideal_quads( PatchData& pd )
00130 {
00131     //            7
00132     //   8 +------+------+ 6
00133     //     |      |      |
00134     //     |      |      |
00135     //     |      |0     |
00136     //   1 +------+------+ 5
00137     //     |      |      |
00138     //     |      |      |
00139     //     |      |      |
00140     //   2 +------+------+ 4
00141     //            3
00142     static const double ideal_quad_verts[] = { 0.0, 0.0,  0.0, -1.0, 0.0,  0.0, -1.0, -1.0, 0.0,
00143                                                0.0, -1.0, 0.0, 1.0,  -1.0, 0.0, 1.0,  0.0,  0.0,
00144                                                1.0, 1.0,  0.0, 0.0,  1.0,  0.0, -1.0, 1.0,  0.0 };
00145     static const size_t ideal_quad_elems[] = { 0, 1, 2, 3, 0, 3, 4, 5, 0, 5, 6, 7, 0, 7, 8, 1 };
00146     const bool fixed[]                     = { false, true, true, true, true, true, true, true, true };
00147 
00148     pd.set_domain( &geomPlane );
00149 
00150     MsqPrintError err( std::cout );
00151     pd.fill( 9, ideal_quad_verts, 4, QUADRILATERAL, ideal_quad_elems, fixed, err );
00152     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00153 }
00154 
00155 void QualityMetricTester::get_ideal_hexes( PatchData& pd )
00156 {
00157     static const double ideal_hex_verts[] = { 0.0,  0.0,  0.0,  -1.0, 0.0,  0.0,  -1.0, -1.0, 0.0,  0.0,  -1.0, 0.0,
00158                                               1.0,  -1.0, 0.0,  1.0,  0.0,  0.0,  1.0,  1.0,  0.0,  0.0,  1.0,  0.0,
00159                                               -1.0, 1.0,  0.0,  0.0,  0.0,  -1.0, -1.0, 0.0,  -1.0, -1.0, -1.0, -1.0,
00160                                               0.0,  -1.0, -1.0, 1.0,  -1.0, -1.0, 1.0,  0.0,  -1.0, 1.0,  1.0,  -1.0,
00161                                               0.0,  1.0,  -1.0, -1.0, 1.0,  -1.0, 0.0,  0.0,  1.0,  -1.0, 0.0,  1.0,
00162                                               -1.0, -1.0, 1.0,  0.0,  -1.0, 1.0,  1.0,  -1.0, 1.0,  1.0,  0.0,  1.0,
00163                                               1.0,  1.0,  1.0,  0.0,  1.0,  1.0,  -1.0, 1.0,  1.0 };
00164     static const size_t ideal_hex_elems[] = { 9, 10, 11, 12, 0,  1,  2,  3,  9, 12, 13, 14, 0,  3,  4,  5,
00165                                               9, 14, 15, 16, 0,  5,  6,  7,  9, 16, 17, 10, 0,  7,  8,  1,
00166                                               0, 1,  2,  3,  18, 19, 20, 21, 0, 3,  4,  5,  18, 21, 22, 23,
00167                                               0, 5,  6,  7,  18, 23, 24, 25, 0, 7,  8,  1,  18, 25, 26, 19 };
00168     bool fixed[27];
00169     fixed[0] = false;
00170     for( size_t f = 1; f < 27; ++f )
00171         fixed[f] = true;
00172 
00173     MsqPrintError err( std::cout );
00174     pd.fill( 27, ideal_hex_verts, 8, HEXAHEDRON, ideal_hex_elems, fixed, err );
00175     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00176 }
00177 
00178 // free_vertex_index:
00179 //  >= 0 : all vertice except specified one are fixed
00180 //    -1 : all vertices are free
00181 //    -2 : first vertex is fixed, all others are free
00182 void QualityMetricTester::get_ideal_element( EntityTopology type, bool unit_area, PatchData& pd, int free_vertex_index )
00183 {
00184     if( TopologyInfo::dimension( type ) == 2 ) pd.set_domain( &geomPlane );
00185 
00186     const size_t n = TopologyInfo::corners( type );
00187     const Vector3D* coords =
00188         unit_area ? unit_element( type, degenHexPyramid ) : unit_edge_element( type, degenHexPyramid );
00189     CPPUNIT_ASSERT( coords != 0 );
00190 
00191     CPPUNIT_ASSERT( sizeof( double ) * 3 == sizeof( Vector3D ) );
00192     const double* elem_verts = reinterpret_cast< const double* >( coords );
00193 
00194     bool* fixed = new bool[n];
00195     std::vector< size_t > conn( n );
00196     for( size_t i = 0; i < n; ++i )
00197     {
00198         conn[i]  = i;
00199         fixed[i] = ( free_vertex_index >= 0 );
00200     }
00201     if( free_vertex_index == -2 )
00202         fixed[0] = true;
00203     else if( free_vertex_index >= 0 )
00204         fixed[free_vertex_index] = false;
00205 
00206     MsqPrintError err( std::cout );
00207     pd.fill( n, elem_verts, 1, type, arrptr( conn ), fixed, err );
00208     delete[] fixed;
00209     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00210 }
00211 void QualityMetricTester::get_ideal_element( EntityTopology type,
00212                                              bool unit_area,
00213                                              PatchData& pd,
00214                                              bool first_vertex_fixed )
00215 {
00216     get_ideal_element( type, unit_area, pd, first_vertex_fixed ? -2 : -1 );
00217 }
00218 
00219 /** Begin with an ideal element, and move the first vertex one half
00220  * of the distnace to the centroid */
00221 void QualityMetricTester::get_nonideal_element( EntityTopology type, PatchData& pd, int free_vtx_type )
00222 {
00223     MsqError err;
00224     get_ideal_element( type, false, pd, free_vtx_type );
00225     const size_t* verts     = pd.element_by_index( 0 ).get_vertex_index_array();
00226     const MsqVertex* coords = pd.get_vertex_array( err );
00227     size_t n                = pd.element_by_index( 0 ).vertex_count();
00228 
00229     Vector3D sum( 0, 0, 0 );
00230     for( unsigned i = 0; i < n; ++i )
00231         sum += coords[verts[i]];
00232 
00233     sum /= n;
00234     sum += coords[verts[0]];
00235     sum *= 0.5;
00236     pd.set_vertex_coordinates( sum, verts[0], err );
00237 }
00238 
00239 void QualityMetricTester::get_nonideal_element( EntityTopology type, PatchData& pd, bool first_vertex_fixed )
00240 {
00241     get_nonideal_element( type, pd, first_vertex_fixed ? -2 : -1 );
00242 }
00243 
00244 /** Collapse first and second vertices of an ideal element */
00245 void QualityMetricTester::get_degenerate_element( EntityTopology type, PatchData& pd )
00246 {
00247     MsqError err;
00248     get_ideal_element( type, false, pd );
00249     const size_t* verts     = pd.element_by_index( 0 ).get_vertex_index_array();
00250     const MsqVertex* coords = pd.get_vertex_array( err );
00251     pd.set_vertex_coordinates( coords[verts[1]], verts[0], err );
00252 }
00253 
00254 /** Create element with zero area/volume
00255  * For quads and hexes, the results are also inverted elements.
00256  */
00257 void QualityMetricTester::get_zero_element( EntityTopology type, PatchData& pd )
00258 {
00259     MsqError err;
00260     get_ideal_element( type, false, pd );
00261     MsqMeshEntity& elem     = pd.element_by_index( 0 );
00262     const size_t* verts     = elem.get_vertex_index_array();
00263     const MsqVertex* coords = pd.get_vertex_array( err );
00264     unsigned i;
00265     Vector3D sum( 0, 0, 0 );
00266 
00267     switch( type )
00268     {
00269         case TRIANGLE:
00270             pd.set_vertex_coordinates( 0.5 * ( coords[verts[0]] + coords[verts[1]] ), verts[2], err );
00271             break;
00272         case QUADRILATERAL:
00273             pd.set_vertex_coordinates( coords[verts[0]], verts[1], err );
00274             pd.set_vertex_coordinates( coords[verts[3]], verts[2], err );
00275             break;
00276         case TETRAHEDRON:
00277             for( i = 0; i < 3; ++i )
00278                 sum += coords[verts[i]];
00279             pd.set_vertex_coordinates( sum / 3.0, verts[3], err );
00280             break;
00281         case PYRAMID:
00282             for( i = 0; i < 4; ++i )
00283                 sum += coords[verts[i]];
00284             pd.set_vertex_coordinates( 0.25 * sum, verts[4], err );
00285             break;
00286         case PRISM:
00287             pd.set_vertex_coordinates( 0.5 * ( coords[verts[0]] + coords[verts[1]] ), verts[2], err );
00288             pd.set_vertex_coordinates( 0.5 * ( coords[verts[3]] + coords[verts[4]] ), verts[5], err );
00289             break;
00290         case HEXAHEDRON:
00291             pd.set_vertex_coordinates( coords[verts[0]], verts[4], err );
00292             pd.set_vertex_coordinates( coords[verts[1]], verts[5], err );
00293             pd.set_vertex_coordinates( coords[verts[2]], verts[6], err );
00294             pd.set_vertex_coordinates( coords[verts[3]], verts[7], err );
00295             break;
00296         default:
00297             CPPUNIT_ASSERT( false );
00298             break;
00299     }
00300 }
00301 
00302 /** Create inverted elements.
00303  * For tri and tet elements, reflect one vertex about the opposite side.
00304  * For all other elements, introduce a concave vertex in one of the
00305  * quadrilateral faces.
00306  */
00307 void QualityMetricTester::get_inverted_element( EntityTopology type, PatchData& pd )
00308 {
00309     MsqError err;
00310     get_ideal_element( type, false, pd );
00311     MsqMeshEntity& elem     = pd.element_by_index( 0 );
00312     const size_t* verts     = elem.get_vertex_index_array();
00313     const MsqVertex* coords = pd.get_vertex_array( err );
00314     unsigned i;
00315     Vector3D sum( 0, 0, 0 );
00316 
00317     //  switch (type) {
00318     //    case TRIANGLE:
00319     //      coords[verts[2]] = coords[verts[0]] + coords[verts[1]] - coords[verts[2]];
00320     //      break;
00321     //    case QUADRILATERAL:
00322     //    case PYRAMID:
00323     //    case HEXAHEDRON:
00324     //      coords[verts[2]] += 0.75 * (coords[verts[0]] - coords[verts[2]]);
00325     //      break;
00326     //    case TETRAHEDRON:
00327     //      for (i = 0; i < 3; ++i)
00328     //        sum += coords[verts[i]];
00329     //      coords[verts[3]] += 0.5*sum - 1.5*coords[verts[3]];
00330     //      break;
00331     //    case PRISM:
00332     //      coords[verts[4]] += 0.75 * (coords[verts[0]] - coords[verts[4]]);
00333     //      break;
00334     //    default:
00335     //      CPPUNIT_ASSERT(false);
00336     //      break;
00337     //  }
00338     if( type == TRIANGLE )
00339         pd.set_vertex_coordinates( coords[verts[0]] + coords[verts[1]] - coords[verts[2]], verts[2], err );
00340     else if( type == QUADRILATERAL || type == PYRAMID || type == HEXAHEDRON )
00341         pd.set_vertex_coordinates( 0.75 * ( coords[verts[0]] - coords[verts[2]] ), verts[2], err );
00342     else if( type == TETRAHEDRON )
00343     {
00344         for( i = 0; i < 3; ++i )
00345             sum += coords[verts[i]];
00346         pd.set_vertex_coordinates( 0.5 * sum - 1.5 * coords[verts[3]], verts[3], err );
00347     }
00348     else if( type == PRISM )
00349         pd.set_vertex_coordinates( 0.75 * ( coords[verts[0]] - coords[verts[4]] ), verts[4], err );
00350     else
00351         CPPUNIT_ASSERT( false );
00352 
00353     if( TopologyInfo::dimension( type ) == 3 || pd.domain_set() )
00354     {
00355         int inverted, total;
00356         elem.check_element_orientation( pd, inverted, total, err );
00357         CPPUNIT_ASSERT( !err );
00358         CPPUNIT_ASSERT( total );
00359     }
00360 }
00361 
00362 void QualityMetricTester::test_type_is_supported( EntityTopology type, QualityMetric* qm )
00363 {
00364     MsqPrintError err( std::cout );
00365     bool rval;
00366     // create patch containing only elements of sepcified type
00367     PatchData pd;
00368     if( mSettings ) pd.attach_settings( mSettings );
00369     get_ideal_element( type, false, pd );
00370     // get list of evaluation locations in patch
00371     std::vector< size_t > handles;
00372     qm->get_evaluations( pd, handles, false, err );
00373     CPPUNIT_ASSERT( !err );
00374     CPPUNIT_ASSERT( !handles.empty() );
00375     // evaluate the metric at one location
00376     double value = -HUGE_VAL;
00377     rval         = qm->evaluate( pd, handles[0], value, err );
00378     CPPUNIT_ASSERT( !err );
00379     CPPUNIT_ASSERT( rval );
00380     // if metric returned input value for 'value' try again
00381     // with a different initial value to verify that the metric
00382     // is setting 'value' to something.
00383     if( value == -HUGE_VAL )
00384     {
00385         value = 0.0;
00386         rval  = qm->evaluate( pd, handles[0], value, err );
00387         CPPUNIT_ASSERT( !err && rval );
00388         CPPUNIT_ASSERT( value != 0.0 );
00389     }
00390     std::vector< size_t > indices;
00391     rval = qm->evaluate_with_indices( pd, handles[0], value, indices, err );
00392     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00393     CPPUNIT_ASSERT( rval );
00394     std::vector< Vector3D > grad;
00395     rval = qm->evaluate_with_gradient( pd, handles[0], value, indices, grad, err );
00396     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00397     CPPUNIT_ASSERT( rval );
00398     std::vector< Matrix3D > hess;
00399     rval = qm->evaluate_with_Hessian( pd, handles[0], value, indices, grad, hess, err );
00400     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00401     CPPUNIT_ASSERT( rval );
00402 }
00403 
00404 void QualityMetricTester::test_type_is_not_supported( EntityTopology type, QualityMetric* qm )
00405 {
00406     MsqError err;
00407     bool rval;
00408     // create patch containing only elements of sepcified type
00409     PatchData pd;
00410     if( mSettings ) pd.attach_settings( mSettings );
00411     get_ideal_element( type, false, pd );
00412     // get list of evaluation locations in patch
00413     std::vector< size_t > handles;
00414     qm->get_evaluations( pd, handles, false, err );
00415     if( err.error_code() == MsqError::UNSUPPORTED_ELEMENT ) return;
00416     CPPUNIT_ASSERT( !err );
00417     if( handles.empty() ) return;
00418     // evaluate the metric at one location
00419     double value;
00420     rval = qm->evaluate( pd, handles[0], value, err );
00421     CPPUNIT_ASSERT_EQUAL( MsqError::UNSUPPORTED_ELEMENT, err.error_code() );
00422     std::vector< size_t > indices;
00423     rval = qm->evaluate_with_indices( pd, handles[0], value, indices, err );
00424     CPPUNIT_ASSERT_EQUAL( MsqError::UNSUPPORTED_ELEMENT, err.error_code() );
00425     std::vector< Vector3D > grad;
00426     rval = qm->evaluate_with_gradient( pd, handles[0], value, indices, grad, err );
00427     CPPUNIT_ASSERT_EQUAL( MsqError::UNSUPPORTED_ELEMENT, err.error_code() );
00428     std::vector< Matrix3D > hess;
00429     rval = qm->evaluate_with_Hessian( pd, handles[0], value, indices, grad, hess, err );
00430     CPPUNIT_ASSERT_EQUAL( MsqError::UNSUPPORTED_ELEMENT, err.error_code() );
00431 }
00432 
00433 void QualityMetricTester::test_supported_element_types( QualityMetric* qm )
00434 {
00435     for( int i = TRIANGLE; i < MIXED; ++i )
00436     {
00437         if( i == POLYGON || i == POLYHEDRON || i == SEPTAHEDRON )
00438             continue;
00439         else if( type_is_supported( (EntityTopology)i ) )
00440             test_type_is_supported( (EntityTopology)i, qm );
00441         else
00442             test_type_is_not_supported( (EntityTopology)i, qm );
00443     }
00444 }
00445 
00446 static void test_evaluate( PatchData& pd, QualityMetric* qm, double value )
00447 {
00448     MsqPrintError err( std::cout );
00449     std::vector< size_t > handles;
00450     qm->get_evaluations( pd, handles, false, err );
00451     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00452     CPPUNIT_ASSERT( !handles.empty() );
00453     for( unsigned i = 0; i < handles.size(); ++i )
00454     {
00455         double qmval;
00456         bool rval = qm->evaluate( pd, handles[i], qmval, err );
00457         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00458         CPPUNIT_ASSERT( rval );
00459         CPPUNIT_ASSERT_DOUBLES_EQUAL( value, qmval, 1e-6 );
00460     }
00461 }
00462 
00463 void QualityMetricTester::test_evaluate_unit_element( QualityMetric* qm, EntityTopology type, double value )
00464 {
00465     PatchData pd;
00466     if( mSettings ) pd.attach_settings( mSettings );
00467     get_ideal_element( type, true, pd );
00468     test_evaluate( pd, qm, value );
00469 }
00470 
00471 void QualityMetricTester::test_evaluate_unit_edge_element( QualityMetric* qm, EntityTopology type, double value )
00472 {
00473     PatchData pd;
00474     if( mSettings ) pd.attach_settings( mSettings );
00475     get_ideal_element( type, false, pd );
00476     test_evaluate( pd, qm, value );
00477 }
00478 
00479 void QualityMetricTester::test_evaluate_unit_tris_about_vertex( QualityMetric* qm, double value )
00480 {
00481     PatchData pd;
00482     if( mSettings ) pd.attach_settings( mSettings );
00483     get_ideal_tris( pd, true );
00484     test_evaluate( pd, qm, value );
00485 }
00486 
00487 void QualityMetricTester::test_evaluate_unit_quads_about_vertex( QualityMetric* qm, double value )
00488 {
00489     PatchData pd;
00490     if( mSettings ) pd.attach_settings( mSettings );
00491     get_ideal_quads( pd );
00492     test_evaluate( pd, qm, value );
00493 }
00494 
00495 void QualityMetricTester::test_evaluate_unit_hexes_about_vertex( QualityMetric* qm, double value )
00496 {
00497     PatchData pd;
00498     if( mSettings ) pd.attach_settings( mSettings );
00499     get_ideal_hexes( pd );
00500     test_evaluate( pd, qm, value );
00501 }
00502 
00503 void QualityMetricTester::test_evaluate_unit_edge_tris_about_vertex( QualityMetric* qm, double value )
00504 {
00505     PatchData pd;
00506     if( mSettings ) pd.attach_settings( mSettings );
00507     get_ideal_tris( pd, false );
00508     test_evaluate( pd, qm, value );
00509 }
00510 
00511 void QualityMetricTester::test_get_element_evaluations( QualityMetric* qm )
00512 {
00513     MsqPrintError err( std::cout );
00514     PatchData pd;
00515     if( mSettings ) pd.attach_settings( mSettings );
00516     unsigned count = 0;
00517     std::vector< size_t > handles;
00518 
00519     if( type_is_supported( HEXAHEDRON ) )
00520     {
00521         get_ideal_hexes( pd );
00522         qm->get_evaluations( pd, handles, false, err );
00523         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00524         CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
00525         std::sort( handles.begin(), handles.end() );
00526         handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
00527         CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
00528         ++count;
00529     }
00530 
00531     if( type_is_supported( QUADRILATERAL ) )
00532     {
00533         get_ideal_quads( pd );
00534         qm->get_evaluations( pd, handles, false, err );
00535         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00536         CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
00537         std::sort( handles.begin(), handles.end() );
00538         handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
00539         CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
00540         ++count;
00541     }
00542 
00543     if( type_is_supported( TRIANGLE ) )
00544     {
00545         get_ideal_tris( pd, false );
00546         qm->get_evaluations( pd, handles, false, err );
00547         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00548         CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
00549         std::sort( handles.begin(), handles.end() );
00550         handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
00551         CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
00552         ++count;
00553     }
00554 
00555     CPPUNIT_ASSERT( count > 0 );
00556 }
00557 
00558 void QualityMetricTester::test_get_vertex_evaluations( QualityMetric* qm )
00559 {
00560     MsqPrintError err( std::cout );
00561     PatchData pd;
00562     if( mSettings ) pd.attach_settings( mSettings );
00563     unsigned count = 0;
00564     std::vector< size_t > handles;
00565 
00566     if( type_is_supported( HEXAHEDRON ) )
00567     {
00568         get_ideal_hexes( pd );
00569         qm->get_evaluations( pd, handles, false, err );
00570         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00571         CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
00572         ++count;
00573     }
00574 
00575     if( type_is_supported( TRIANGLE ) )
00576     {
00577         get_ideal_tris( pd, false );
00578         qm->get_evaluations( pd, handles, false, err );
00579         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00580         CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
00581         ++count;
00582     }
00583 
00584     if( type_is_supported( QUADRILATERAL ) )
00585     {
00586         get_ideal_quads( pd );
00587         qm->get_evaluations( pd, handles, false, err );
00588         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00589         CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
00590         ++count;
00591     }
00592 
00593     CPPUNIT_ASSERT( count > 0 );
00594 }
00595 
00596 void QualityMetricTester::test_get_sample_evaluations( QualityMetric* qm )
00597 {
00598     MsqPrintError err( std::cout );
00599     PatchData pd;
00600     if( mSettings ) pd.attach_settings( mSettings );
00601     unsigned count = 0;
00602     std::vector< size_t > handles;
00603 
00604     if( type_is_supported( HEXAHEDRON ) )
00605     {
00606         get_ideal_hexes( pd );
00607         size_t expected_evals = 0;
00608         for( size_t i = 0; i < pd.num_elements(); ++i )
00609             expected_evals += pd.get_samples( i ).num_nodes();
00610 
00611         qm->get_evaluations( pd, handles, false, err );
00612         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00613         CPPUNIT_ASSERT_EQUAL( expected_evals, handles.size() );
00614 
00615         std::sort( handles.begin(), handles.end() );
00616         handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
00617         CPPUNIT_ASSERT_EQUAL( expected_evals, handles.size() );
00618 
00619         ++count;
00620     }
00621 
00622     if( type_is_supported( TRIANGLE ) )
00623     {
00624         get_ideal_tris( pd, false );
00625         size_t expected_evals = 0;
00626         for( size_t i = 0; i < pd.num_elements(); ++i )
00627             expected_evals += pd.get_samples( i ).num_nodes();
00628 
00629         qm->get_evaluations( pd, handles, false, err );
00630         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00631         CPPUNIT_ASSERT_EQUAL( expected_evals, handles.size() );
00632 
00633         std::sort( handles.begin(), handles.end() );
00634         handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
00635         CPPUNIT_ASSERT_EQUAL( expected_evals, handles.size() );
00636 
00637         ++count;
00638     }
00639 
00640     if( type_is_supported( QUADRILATERAL ) )
00641     {
00642         get_ideal_quads( pd );
00643         size_t expected_evals = 0;
00644         for( size_t i = 0; i < pd.num_elements(); ++i )
00645             expected_evals += pd.get_samples( i ).num_nodes();
00646 
00647         qm->get_evaluations( pd, handles, false, err );
00648         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00649         CPPUNIT_ASSERT_EQUAL( expected_evals, handles.size() );
00650 
00651         std::sort( handles.begin(), handles.end() );
00652         handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
00653         CPPUNIT_ASSERT_EQUAL( expected_evals, handles.size() );
00654 
00655         ++count;
00656     }
00657 
00658     CPPUNIT_ASSERT( count > 0 );
00659 }
00660 void get_ideal_element( EntityTopology type, bool unit_area, PatchData& pd, int free_vertex_index );
00661 
00662 void QualityMetricTester::test_get_edge_evaluations( EdgeQM* qm )
00663 {
00664     MsqPrintError err( std::cout );
00665     PatchData pd;
00666     if( mSettings ) pd.attach_settings( mSettings );
00667     std::vector< size_t > handles;
00668 
00669     CPPUNIT_ASSERT( !types.empty() );
00670     for( size_t i = 0; i < types.size(); ++i )
00671     {
00672         get_ideal_element( types[i], true, pd );
00673         handles.clear();
00674         qm->get_evaluations( pd, handles, false, err );
00675         ASSERT_NO_ERROR( err );
00676         CPPUNIT_ASSERT_EQUAL( TopologyInfo::edges( types[i] ), (unsigned)handles.size() );
00677     }
00678 }
00679 
00680 void QualityMetricTester::test_get_in_element_evaluations( ElemSampleQM* qm )
00681 {
00682     MsqPrintError err( std::cout );
00683     PatchData pd;
00684     if( mSettings ) pd.attach_settings( mSettings );
00685     std::vector< size_t > handles1, handles2;
00686 
00687     CPPUNIT_ASSERT( !types.empty() );
00688 
00689     get_ideal_element( types[0], false, pd );
00690     CPPUNIT_ASSERT( pd.num_elements() == 1 );
00691     qm->get_evaluations( pd, handles1, false, err );
00692     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00693     CPPUNIT_ASSERT( !handles1.empty() );
00694     qm->get_element_evaluations( pd, 0, handles2, err );
00695     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00696     CPPUNIT_ASSERT( !handles2.empty() );
00697     std::sort( handles1.begin(), handles1.end() );
00698     std::sort( handles2.begin(), handles2.end() );
00699     CPPUNIT_ASSERT( handles1 == handles2 );
00700 
00701     get_ideal_element( types.back(), false, pd );
00702     CPPUNIT_ASSERT( pd.num_elements() == 1 );
00703     qm->get_evaluations( pd, handles1, false, err );
00704     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00705     CPPUNIT_ASSERT( !handles1.empty() );
00706     qm->get_element_evaluations( pd, 0, handles2, err );
00707     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00708     CPPUNIT_ASSERT( !handles2.empty() );
00709     std::sort( handles1.begin(), handles1.end() );
00710     std::sort( handles2.begin(), handles2.end() );
00711     CPPUNIT_ASSERT( handles1 == handles2 );
00712 }
00713 
00714 void QualityMetricTester::test_get_indices_fixed( QualityMetric* qm )
00715 {
00716     MsqPrintError err( std::cout );
00717     PatchData pd1, pd2;
00718     if( mSettings )
00719     {
00720         pd1.attach_settings( mSettings );
00721         pd2.attach_settings( mSettings );
00722     }
00723     std::vector< size_t > handles1, handles2, indices1, indices2;
00724     double qm_val1, qm_val2;
00725 
00726     CPPUNIT_ASSERT( !types.empty() );
00727     // get element with no fixed vertices
00728     get_ideal_element( types.back(), false, pd1, false );
00729     qm->get_evaluations( pd1, handles1, false, err );
00730     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00731     CPPUNIT_ASSERT( !handles1.empty() );
00732     // call evaluate with indices
00733     size_t handle = handles1.back();
00734     qm->evaluate_with_indices( pd1, handle, qm_val1, indices1, err );
00735     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00736     CPPUNIT_ASSERT( !indices1.empty() );
00737     // get element with one fixed vertex
00738     get_ideal_element( types.back(), false, pd2, true );
00739     qm->get_evaluations( pd2, handles2, false, err );
00740     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00741     CPPUNIT_ASSERT( !handles2.empty() );
00742     // If vertex-based metric, need to find corresponding handle
00743     // in second patch.  The vertex order changed as a result of
00744     // one of the vertices beging fixed in the second patch.
00745     if( qm->get_metric_type() == QualityMetric::VERTEX_BASED )
00746     {
00747         handle              = handles1.back();
00748         const size_t* conn1 = pd1.element_by_index( 0 ).get_vertex_index_array();
00749         const size_t* conn2 = pd2.element_by_index( 0 ).get_vertex_index_array();
00750         const size_t len    = TopologyInfo::corners( types.back() );
00751         const size_t* ptr   = std::find( conn1, conn1 + len, handle );
00752         handle              = conn2[ptr - conn1];
00753     }
00754     // call evaluate with indices
00755     qm->evaluate_with_indices( pd2, handle, qm_val2, indices2, err );
00756     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00757     CPPUNIT_ASSERT( !indices2.empty() );
00758     // make sure we got the same QM value
00759     // (shoudn't be affected by fixed/non-fixed)
00760     CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
00761     // should have gotten back one less index (for element-based metrics)
00762     if( handles1.size() == 1 )
00763     {
00764         CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() + 1 );
00765     }
00766     // indices2 shouldn't contain any fixed vertices
00767     std::sort( indices2.begin(), indices2.end() );
00768     CPPUNIT_ASSERT( indices2.back() < pd2.num_free_vertices() );
00769     // indices2 shouldn/t contain any duplicates
00770     CPPUNIT_ASSERT( std::unique( indices2.begin(), indices2.end() ) == indices2.end() );
00771 }
00772 
00773 void QualityMetricTester::test_get_element_indices( QualityMetric* qm )
00774 {
00775     MsqPrintError err( std::cout );
00776     PatchData pd;
00777     if( mSettings ) pd.attach_settings( mSettings );
00778     std::vector< size_t > handles, indices, verts;
00779     double qm_val;
00780 
00781     CPPUNIT_ASSERT( qm->get_metric_type() == QualityMetric::ELEMENT_BASED );
00782     CPPUNIT_ASSERT( !types.empty() );
00783 
00784     for( size_t i = 0; i < types.size(); ++i )
00785     {
00786         // construct patch w/ one fixed vertex
00787         get_ideal_element( types[i], false, pd, true );
00788         qm->get_evaluations( pd, handles, false, err );
00789         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00790         CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
00791         // get vertex list
00792         verts.clear();
00793         pd.element_by_index( 0 ).get_vertex_indices( verts );
00794         // remove fixed vertex
00795         for( size_t j = 0; j < verts.size(); ++j )
00796             if( verts[j] >= pd.num_free_vertices() )
00797             {
00798                 verts.erase( verts.begin() + j );
00799                 break;
00800             }
00801         // evaluate metric
00802         qm->evaluate_with_indices( pd, handles[0], qm_val, indices, err );
00803         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00804         // index list should match vertex list (element metric
00805         // should depend on all free vertices in element).
00806         CPPUNIT_ASSERT_EQUAL( verts.size(), indices.size() );
00807         std::sort( verts.begin(), verts.end() );
00808         std::sort( indices.begin(), indices.end() );
00809         CPPUNIT_ASSERT( verts == indices );
00810     }
00811 }
00812 
00813 void QualityMetricTester::test_get_vertex_indices( QualityMetric* qm )
00814 {
00815     MsqPrintError err( std::cout );
00816     PatchData pd;
00817     if( mSettings ) pd.attach_settings( mSettings );
00818     std::vector< size_t > indices, verts;
00819     double qm_val;
00820 
00821     CPPUNIT_ASSERT( qm->get_metric_type() == QualityMetric::VERTEX_BASED );
00822     CPPUNIT_ASSERT( !types.empty() );
00823 
00824     for( size_t i = 0; i < types.size(); ++i )
00825     {
00826         // construct patch w/ one fixed vertex
00827         get_ideal_element( types[i], false, pd, true );
00828         // get adjacent vertices
00829         verts.clear();
00830         pd.get_adjacent_vertex_indices( 0, verts, err );
00831         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00832         // remove fixed vertex
00833         for( size_t j = 0; j < verts.size(); ++j )
00834             if( verts[j] >= pd.num_free_vertices() )
00835             {
00836                 verts.erase( verts.begin() + j );
00837                 break;
00838             }
00839         verts.push_back( 0 );
00840         // evaluate metric
00841         qm->evaluate_with_indices( pd, 0, qm_val, indices, err );
00842         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00843         // vertex metric should depend on all adjacent free vertices
00844         CPPUNIT_ASSERT_EQUAL( verts.size(), indices.size() );
00845         std::sort( verts.begin(), verts.end() );
00846         std::sort( indices.begin(), indices.end() );
00847         CPPUNIT_ASSERT( verts == indices );
00848     }
00849 }
00850 
00851 void QualityMetricTester::test_get_edge_indices( EdgeQM* qm )
00852 {
00853     MsqPrintError err( std::cout );
00854     PatchData pd, pd2;
00855     if( mSettings )
00856     {
00857         pd.attach_settings( mSettings );
00858         pd2.attach_settings( mSettings );
00859     }
00860     std::vector< size_t > handles, indices, indices2, verts;
00861     std::vector< size_t >::iterator it;
00862     double qm_val;
00863 
00864     CPPUNIT_ASSERT( qm->get_metric_type() == QualityMetric::VERTEX_BASED );
00865     CPPUNIT_ASSERT( !types.empty() );
00866 
00867     for( size_t i = 0; i < types.size(); ++i )
00868     {
00869         // construct patch w/ no free vertices
00870         get_ideal_element( types[i], false, pd, false );
00871 
00872         qm->get_evaluations( pd, handles, false, err );
00873         ASSERT_NO_ERROR( err );
00874         CPPUNIT_ASSERT_EQUAL( (size_t)TopologyInfo::edges( types[i] ), handles.size() );
00875 
00876         for( size_t j = 0; j < handles.size(); ++j )
00877         {
00878             // evaluate metric
00879             qm->evaluate_with_indices( pd, handles[j], qm_val, indices, err );
00880             ASSERT_NO_ERROR( err );
00881             // evaluation at each edge should depend on at least the two end vertices
00882             CPPUNIT_ASSERT( indices.size() >= (size_t)2 );
00883         }
00884     }
00885 
00886     for( size_t i = 0; i < types.size(); ++i )
00887     {
00888         // construct patch w/ one free vertex
00889         get_ideal_element( types[i], false, pd, true );
00890         const size_t fixed_vertex = pd.num_free_vertices();
00891 
00892         qm->get_evaluations( pd, handles, false, err );
00893         ASSERT_NO_ERROR( err );
00894         CPPUNIT_ASSERT_EQUAL( (size_t)TopologyInfo::edges( types[i] ), handles.size() );
00895 
00896         for( size_t j = 0; j < handles.size(); ++j )
00897         {
00898             // evaluate metric
00899             qm->evaluate_with_indices( pd, handles[j], qm_val, indices, err );
00900             ASSERT_NO_ERROR( err );
00901             // evaluation at each edge should depend on at least the two end vertices
00902             CPPUNIT_ASSERT( !indices.empty() );
00903 
00904             // indices should never contain the index of a fixed vertex
00905             it = std::find( indices.begin(), indices.end(), fixed_vertex );
00906             CPPUNIT_ASSERT( it == indices.end() );
00907         }
00908     }
00909 }
00910 
00911 void QualityMetricTester::test_get_sample_indices( QualityMetric* qm )
00912 {
00913     MsqPrintError err( std::cout );
00914     PatchData pd;
00915     if( mSettings ) pd.attach_settings( mSettings );
00916     std::vector< size_t > handles, indices, verts;
00917     double qm_val;
00918 
00919     CPPUNIT_ASSERT( qm->get_metric_type() == QualityMetric::ELEMENT_BASED );
00920     CPPUNIT_ASSERT( !types.empty() );
00921 
00922     for( size_t i = 0; i < types.size(); ++i )
00923     {
00924         // get evaluation locations
00925         get_ideal_element( types[i], false, pd, true );
00926         const size_t count = pd.get_samples( 0 ).num_nodes();
00927         CPPUNIT_ASSERT( count > 0 );
00928         qm->get_evaluations( pd, handles, false, err );
00929         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00930         // should have one evaluation for each sample point
00931         CPPUNIT_ASSERT_EQUAL( count, handles.size() );
00932         // get vertex list
00933         verts.clear();
00934         pd.element_by_index( 0 ).get_vertex_indices( verts );
00935         // remove fixed vertex
00936         for( size_t j = 0; j < verts.size(); ++j )
00937             if( verts[j] >= pd.num_free_vertices() )
00938             {
00939                 verts.erase( verts.begin() + j );
00940                 break;
00941             }
00942         // evaluate metric
00943         qm->evaluate_with_indices( pd, handles[0], qm_val, indices, err );
00944         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00945         // only one fixed vertex, and presumably every possible evaluation
00946         // depends on more than one vertex, so should be at least one
00947         // other, non-fixed, vertex.
00948         CPPUNIT_ASSERT( indices.size() > 0 );
00949         // make sure no duplicates
00950         std::sort( indices.begin(), indices.end() );
00951         CPPUNIT_ASSERT( std::unique( indices.begin(), indices.end() ) == indices.end() );
00952         // check that all indices are in vertex list
00953         CPPUNIT_ASSERT( indices.size() <= verts.size() );
00954         for( std::vector< size_t >::iterator j = indices.begin(); j != indices.end(); ++j )
00955             CPPUNIT_ASSERT( std::find( verts.begin(), verts.end(), *j ) != verts.end() );
00956     }
00957 }
00958 
00959 void QualityMetricTester::compare_eval_and_eval_with_indices( QualityMetric* qm )
00960 {
00961     MsqPrintError err( std::cout );
00962     PatchData pd;
00963     if( mSettings ) pd.attach_settings( mSettings );
00964     std::vector< size_t > handles, indices;
00965     double qm_val1, qm_val2;
00966     bool rval;
00967 
00968     CPPUNIT_ASSERT( !types.empty() );
00969     for( size_t i = 0; i < types.size(); ++i )
00970     {
00971         get_ideal_element( types[i], false, pd );
00972         qm->get_evaluations( pd, handles, false, err );
00973         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00974         CPPUNIT_ASSERT( !handles.empty() );
00975         for( size_t j = 0; j < handles.size(); ++j )
00976         {
00977             rval = qm->evaluate( pd, handles[j], qm_val1, err );
00978             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00979             CPPUNIT_ASSERT( rval );
00980             rval = qm->evaluate_with_indices( pd, handles[j], qm_val2, indices, err );
00981             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00982             CPPUNIT_ASSERT( rval );
00983 
00984             CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
00985         }
00986     }
00987 }
00988 
00989 void QualityMetricTester::compare_eval_with_indices_and_eval_with_gradient( QualityMetric* qm, PatchData& pd )
00990 {
00991     MsqPrintError err( std::cout );
00992     std::vector< size_t > handles, indices1, indices2;
00993     std::vector< Vector3D > grad;
00994     double qm_val1, qm_val2;
00995     bool rval;
00996 
00997     qm->get_evaluations( pd, handles, false, err );
00998     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00999     CPPUNIT_ASSERT( !handles.empty() );
01000     for( size_t j = 0; j < handles.size(); ++j )
01001     {
01002         rval = qm->evaluate_with_indices( pd, handles[j], qm_val1, indices1, err );
01003         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01004         CPPUNIT_ASSERT( rval );
01005         rval = qm->evaluate_with_gradient( pd, handles[j], qm_val2, indices2, grad, err );
01006         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01007         CPPUNIT_ASSERT( rval );
01008 
01009         CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
01010         CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
01011 
01012         std::sort( indices1.begin(), indices1.end() );
01013         std::sort( indices2.begin(), indices2.end() );
01014         CPPUNIT_ASSERT( indices1 == indices2 );
01015         CPPUNIT_ASSERT( !indices1.empty() );
01016     }
01017 }
01018 
01019 void QualityMetricTester::compare_eval_with_indices_and_eval_with_gradient( QualityMetric* qm )
01020 {
01021     PatchData pd;
01022     if( mSettings ) pd.attach_settings( mSettings );
01023 
01024     CPPUNIT_ASSERT( !types.empty() );
01025     for( size_t i = 0; i < types.size(); ++i )
01026     {
01027         get_ideal_element( types[i], false, pd, true );
01028         compare_eval_with_indices_and_eval_with_gradient( qm, pd );
01029     }
01030 }
01031 
01032 void QualityMetricTester::compare_eval_with_indices_and_eval_with_hessian( QualityMetric* qm, PatchData& pd )
01033 {
01034     MsqPrintError err( std::cout );
01035     std::vector< size_t > handles, indices1, indices2;
01036     std::vector< Vector3D > grad;
01037     std::vector< Matrix3D > Hess;
01038     double qm_val1, qm_val2;
01039     bool rval;
01040 
01041     qm->get_evaluations( pd, handles, false, err );
01042     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01043     CPPUNIT_ASSERT( !handles.empty() );
01044     for( size_t j = 0; j < handles.size(); ++j )
01045     {
01046         rval = qm->evaluate_with_indices( pd, handles[j], qm_val1, indices1, err );
01047         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01048         CPPUNIT_ASSERT( rval );
01049         rval = qm->evaluate_with_Hessian( pd, handles[j], qm_val2, indices2, grad, Hess, err );
01050         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01051         CPPUNIT_ASSERT( rval );
01052 
01053         CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
01054         CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
01055 
01056         std::sort( indices1.begin(), indices1.end() );
01057         std::sort( indices2.begin(), indices2.end() );
01058         CPPUNIT_ASSERT( indices1 == indices2 );
01059         CPPUNIT_ASSERT( !indices1.empty() );
01060     }
01061 }
01062 
01063 void QualityMetricTester::compare_eval_with_indices_and_eval_with_hessian( QualityMetric* qm )
01064 {
01065     PatchData pd;
01066     if( mSettings ) pd.attach_settings( mSettings );
01067 
01068     CPPUNIT_ASSERT( !types.empty() );
01069     for( size_t i = 0; i < types.size(); ++i )
01070     {
01071         get_nonideal_element( types[i], pd, true );
01072         compare_eval_with_indices_and_eval_with_hessian( qm, pd );
01073     }
01074 }
01075 
01076 void QualityMetricTester::compare_eval_with_indices_and_eval_with_diagonal( QualityMetric* qm, PatchData& pd )
01077 {
01078     MsqPrintError err( std::cout );
01079     std::vector< size_t > handles, indices1, indices2;
01080     std::vector< Vector3D > grad;
01081     std::vector< SymMatrix3D > Hess;
01082     double qm_val1, qm_val2;
01083     bool rval;
01084 
01085     qm->get_evaluations( pd, handles, false, err );
01086     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01087     CPPUNIT_ASSERT( !handles.empty() );
01088     for( size_t j = 0; j < handles.size(); ++j )
01089     {
01090         rval = qm->evaluate_with_indices( pd, handles[j], qm_val1, indices1, err );
01091         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01092         CPPUNIT_ASSERT( rval );
01093         rval = qm->evaluate_with_Hessian_diagonal( pd, handles[j], qm_val2, indices2, grad, Hess, err );
01094         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01095         CPPUNIT_ASSERT( rval );
01096 
01097         CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
01098         CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
01099 
01100         std::sort( indices1.begin(), indices1.end() );
01101         std::sort( indices2.begin(), indices2.end() );
01102         CPPUNIT_ASSERT( indices1 == indices2 );
01103         CPPUNIT_ASSERT( !indices1.empty() );
01104     }
01105 }
01106 
01107 void QualityMetricTester::compare_eval_with_indices_and_eval_with_diagonal( QualityMetric* qm )
01108 {
01109     PatchData pd;
01110     if( mSettings ) pd.attach_settings( mSettings );
01111 
01112     CPPUNIT_ASSERT( !types.empty() );
01113     for( size_t i = 0; i < types.size(); ++i )
01114     {
01115         get_nonideal_element( types[i], pd, true );
01116         compare_eval_with_indices_and_eval_with_diagonal( qm, pd );
01117     }
01118 }
01119 
01120 void QualityMetricTester::compare_eval_with_grad_and_eval_with_hessian( QualityMetric* qm, PatchData& pd )
01121 {
01122     MsqPrintError err( std::cout );
01123     std::vector< size_t > handles, indices1, indices2;
01124     std::vector< Vector3D > grad1, grad2;
01125     std::vector< Matrix3D > Hess;
01126     double qm_val1, qm_val2;
01127     bool rval;
01128 
01129     qm->get_evaluations( pd, handles, false, err );
01130     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01131     CPPUNIT_ASSERT( !handles.empty() );
01132     for( size_t j = 0; j < handles.size(); ++j )
01133     {
01134         rval = qm->evaluate_with_gradient( pd, handles[j], qm_val1, indices1, grad1, err );
01135         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01136         CPPUNIT_ASSERT( rval );
01137         rval = qm->evaluate_with_Hessian( pd, handles[j], qm_val2, indices2, grad2, Hess, err );
01138         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01139         CPPUNIT_ASSERT( rval );
01140 
01141         CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
01142         CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
01143 
01144         for( size_t k = 0; k < indices1.size(); ++k )
01145         {
01146             std::vector< size_t >::iterator it = std::find( indices2.begin(), indices2.end(), indices1[k] );
01147             CPPUNIT_ASSERT( it != indices2.end() );
01148             size_t m = it - indices2.begin();
01149             CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[k], grad2[m], 1e-6 );
01150         }
01151     }
01152 }
01153 
01154 void QualityMetricTester::compare_eval_with_grad_and_eval_with_hessian( QualityMetric* qm )
01155 {
01156     PatchData pd;
01157     if( mSettings ) pd.attach_settings( mSettings );
01158 
01159     CPPUNIT_ASSERT( !types.empty() );
01160     for( size_t i = 0; i < types.size(); ++i )
01161     {
01162         get_nonideal_element( types[i], pd, true );
01163         compare_eval_with_grad_and_eval_with_hessian( qm, pd );
01164     }
01165 }
01166 
01167 void QualityMetricTester::compare_eval_with_grad_and_eval_with_diagonal( QualityMetric* qm, PatchData& pd )
01168 {
01169     MsqPrintError err( std::cout );
01170     std::vector< size_t > handles, indices1, indices2;
01171     std::vector< Vector3D > grad1, grad2;
01172     std::vector< SymMatrix3D > Hess;
01173     double qm_val1, qm_val2;
01174     bool rval;
01175 
01176     qm->get_evaluations( pd, handles, false, err );
01177     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01178     CPPUNIT_ASSERT( !handles.empty() );
01179     for( size_t j = 0; j < handles.size(); ++j )
01180     {
01181         rval = qm->evaluate_with_gradient( pd, handles[j], qm_val1, indices1, grad1, err );
01182         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01183         CPPUNIT_ASSERT( rval );
01184         rval = qm->evaluate_with_Hessian_diagonal( pd, handles[j], qm_val2, indices2, grad2, Hess, err );
01185         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01186         CPPUNIT_ASSERT( rval );
01187 
01188         CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
01189         CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
01190 
01191         for( size_t k = 0; k < indices1.size(); ++k )
01192         {
01193             std::vector< size_t >::iterator it = std::find( indices2.begin(), indices2.end(), indices1[k] );
01194             CPPUNIT_ASSERT( it != indices2.end() );
01195             size_t m = it - indices2.begin();
01196             CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[k], grad2[m], 1e-6 );
01197         }
01198     }
01199 }
01200 
01201 void QualityMetricTester::compare_eval_with_grad_and_eval_with_diagonal( QualityMetric* qm )
01202 {
01203     PatchData pd;
01204     if( mSettings ) pd.attach_settings( mSettings );
01205 
01206     CPPUNIT_ASSERT( !types.empty() );
01207     for( size_t i = 0; i < types.size(); ++i )
01208     {
01209         get_nonideal_element( types[i], pd, true );
01210         compare_eval_with_grad_and_eval_with_diagonal( qm, pd );
01211     }
01212 }
01213 
01214 void QualityMetricTester::compare_eval_with_diag_and_eval_with_hessian( QualityMetric* qm, PatchData& pd )
01215 {
01216     MsqPrintError err( std::cout );
01217     std::vector< size_t > handles, indices1, indices2;
01218     std::vector< Vector3D > grad1, grad2;
01219     std::vector< SymMatrix3D > hess1;
01220     std::vector< Matrix3D > hess2;
01221     double qm_val1, qm_val2;
01222     bool rval;
01223 
01224     qm->get_evaluations( pd, handles, false, err );
01225     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01226     CPPUNIT_ASSERT( !handles.empty() );
01227     for( size_t j = 0; j < handles.size(); ++j )
01228     {
01229         rval = qm->evaluate_with_Hessian_diagonal( pd, handles[j], qm_val1, indices1, grad1, hess1, err );
01230         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01231         CPPUNIT_ASSERT( rval );
01232         rval = qm->evaluate_with_Hessian( pd, handles[j], qm_val2, indices2, grad2, hess2, err );
01233         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01234         CPPUNIT_ASSERT( rval );
01235 
01236         CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
01237 
01238         CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
01239         CPPUNIT_ASSERT_EQUAL( grad1.size(), grad2.size() );
01240         CPPUNIT_ASSERT_EQUAL( hess2.size(), hess1.size() * ( hess1.size() + 1 ) / 2 );
01241         unsigned h2step                             = indices2.size();
01242         std::vector< Matrix3D >::const_iterator h2i = hess2.begin();
01243 
01244         for( size_t k = 0; k < indices2.size(); ++k )
01245         {
01246             std::vector< size_t >::iterator it = std::find( indices1.begin(), indices1.end(), indices2[k] );
01247             CPPUNIT_ASSERT( it != indices1.end() );
01248             size_t m = it - indices1.begin();
01249             CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[m], grad2[k], 5e-5 );
01250             CPPUNIT_ASSERT_MATRICES_EQUAL( Matrix3D( hess1[m] ), *h2i, 5e-5 );
01251             h2i += h2step--;
01252         }
01253     }
01254 }
01255 
01256 void QualityMetricTester::compare_eval_with_diag_and_eval_with_hessian( QualityMetric* qm )
01257 {
01258     PatchData pd;
01259     if( mSettings ) pd.attach_settings( mSettings );
01260 
01261     CPPUNIT_ASSERT( !types.empty() );
01262     for( size_t i = 0; i < types.size(); ++i )
01263     {
01264         get_nonideal_element( types[i], pd, true );
01265         compare_eval_with_diag_and_eval_with_hessian( qm, pd );
01266     }
01267 }
01268 
01269 void QualityMetricTester::compare_analytical_and_numerical_gradients( QualityMetric* qm, PatchData& pd )
01270 {
01271     MsqPrintError err( std::cout );
01272     std::vector< size_t > handles, indices1, indices2;
01273     std::vector< Vector3D > grad1, grad2;
01274     double qm_val1, qm_val2;
01275     bool rval;
01276 
01277     qm->get_evaluations( pd, handles, false, err );
01278     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01279     CPPUNIT_ASSERT( !handles.empty() );
01280     for( size_t j = 0; j < handles.size(); ++j )
01281     {
01282         rval = qm->QualityMetric::evaluate_with_gradient( pd, handles[j], qm_val1, indices1, grad1, err );
01283         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01284         CPPUNIT_ASSERT( rval );
01285         rval = qm->evaluate_with_gradient( pd, handles[j], qm_val2, indices2, grad2, err );
01286         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01287         CPPUNIT_ASSERT( rval );
01288 
01289         CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
01290         CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
01291         CPPUNIT_ASSERT( !indices1.empty() );
01292 
01293         std::vector< size_t >::iterator it1, it2;
01294         for( it1 = indices1.begin(); it1 != indices1.end(); ++it1 )
01295         {
01296             it2 = std::find( indices2.begin(), indices2.end(), *it1 );
01297             CPPUNIT_ASSERT( it2 != indices2.end() );
01298 
01299             size_t idx1 = it1 - indices1.begin();
01300             size_t idx2 = it2 - indices2.begin();
01301             CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[idx1], grad2[idx2], 0.01 );
01302         }
01303     }
01304 }
01305 
01306 void QualityMetricTester::compare_analytical_and_numerical_gradients( QualityMetric* qm )
01307 {
01308     PatchData pd;
01309     if( mSettings ) pd.attach_settings( mSettings );
01310 
01311     CPPUNIT_ASSERT( !types.empty() );
01312     for( size_t i = 0; i < types.size(); ++i )
01313     {
01314         get_nonideal_element( types[i], pd, true );
01315         compare_analytical_and_numerical_gradients( qm, pd );
01316     }
01317 }
01318 
01319 void QualityMetricTester::compare_analytical_and_numerical_hessians( QualityMetric* qm, PatchData& pd )
01320 {
01321     MsqPrintError err( std::cout );
01322     std::vector< size_t > handles, indices1, indices2;
01323     std::vector< Vector3D > grad1, grad2;
01324     std::vector< Matrix3D > Hess1, Hess2;
01325     double qm_val1, qm_val2;
01326     bool rval;
01327 
01328     qm->get_evaluations( pd, handles, false, err );
01329     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01330     CPPUNIT_ASSERT( !handles.empty() );
01331     for( size_t j = 0; j < handles.size(); ++j )
01332     {
01333         rval = qm->QualityMetric::evaluate_with_Hessian( pd, handles[j], qm_val1, indices1, grad1, Hess1, err );
01334         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01335         CPPUNIT_ASSERT( rval );
01336         rval = qm->evaluate_with_Hessian( pd, handles[j], qm_val2, indices2, grad2, Hess2, err );
01337         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01338         CPPUNIT_ASSERT( rval );
01339 
01340         CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
01341         CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
01342         CPPUNIT_ASSERT( !indices1.empty() );
01343 
01344         std::vector< size_t >::iterator it;
01345         unsigned h = 0;
01346         for( unsigned r = 0; r < indices1.size(); ++r )
01347         {
01348             it = std::find( indices2.begin(), indices2.end(), indices1[r] );
01349             CPPUNIT_ASSERT( it != indices2.end() );
01350             unsigned r2 = it - indices2.begin();
01351 
01352             for( unsigned c = r; c < indices1.size(); ++c, ++h )
01353             {
01354                 it = std::find( indices2.begin(), indices2.end(), indices1[c] );
01355                 CPPUNIT_ASSERT( it != indices2.end() );
01356                 unsigned c2 = it - indices2.begin();
01357 
01358                 unsigned h2;
01359                 if( r2 <= c2 )
01360                     h2 = indices2.size() * r - r * ( r + 1 ) / 2 + c;
01361                 else
01362                     h2 = indices2.size() * c - c * ( c + 1 ) / 2 + r;
01363 
01364                 // if (!utest_mat_equal(Hess1[h],Hess2[h2],0.001))
01365                 //  assert(false);
01366                 CPPUNIT_ASSERT_MATRICES_EQUAL( Hess1[h], Hess2[h2], 0.001 );
01367             }
01368         }
01369     }
01370 }
01371 
01372 void QualityMetricTester::compare_analytical_and_numerical_hessians( QualityMetric* qm )
01373 {
01374     PatchData pd;
01375     if( mSettings ) pd.attach_settings( mSettings );
01376 
01377     CPPUNIT_ASSERT( !types.empty() );
01378     for( size_t i = 0; i < types.size(); ++i )
01379     {
01380         get_nonideal_element( types[i], pd, true );
01381         compare_analytical_and_numerical_hessians( qm, pd );
01382     }
01383 }
01384 
01385 void QualityMetricTester::compare_analytical_and_numerical_diagonals( QualityMetric* qm, PatchData& pd )
01386 {
01387     MsqPrintError err( std::cout );
01388     std::vector< size_t > handles, indices1, indices2;
01389     std::vector< Vector3D > grad1, grad2;
01390     std::vector< Matrix3D > Hess1;
01391     std::vector< SymMatrix3D > Hess2;
01392     double qm_val1, qm_val2;
01393     bool rval;
01394 
01395     qm->get_evaluations( pd, handles, false, err );
01396     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01397     CPPUNIT_ASSERT( !handles.empty() );
01398     for( size_t j = 0; j < handles.size(); ++j )
01399     {
01400         rval = qm->QualityMetric::evaluate_with_Hessian( pd, handles[j], qm_val1, indices1, grad1, Hess1, err );
01401         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01402         CPPUNIT_ASSERT( rval );
01403         rval = qm->evaluate_with_Hessian_diagonal( pd, handles[j], qm_val2, indices2, grad2, Hess2, err );
01404         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01405         CPPUNIT_ASSERT( rval );
01406 
01407         CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
01408         CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
01409         CPPUNIT_ASSERT( !indices1.empty() );
01410         CPPUNIT_ASSERT_EQUAL( indices1.size() * ( indices1.size() + 1 ) / 2, Hess1.size() );
01411         CPPUNIT_ASSERT_EQUAL( indices2.size(), Hess2.size() );
01412 
01413         size_t h = 0;
01414         std::vector< size_t >::iterator it;
01415         for( unsigned r = 0; r < indices1.size(); ++r )
01416         {
01417             it = std::find( indices2.begin(), indices2.end(), indices1[r] );
01418             CPPUNIT_ASSERT( it != indices2.end() );
01419             unsigned r2 = it - indices2.begin();
01420             // if (!utest_mat_equal(Hess1[h],Hess2[r2],0.001))
01421             //  assert(false);
01422             CPPUNIT_ASSERT_MATRICES_EQUAL( Hess1[h], Hess2[r2], 0.001 );
01423             h += indices1.size() - r;
01424         }
01425     }
01426 }
01427 
01428 void QualityMetricTester::compare_analytical_and_numerical_diagonals( QualityMetric* qm )
01429 {
01430     PatchData pd;
01431     if( mSettings ) pd.attach_settings( mSettings );
01432     MsqError err;
01433 
01434     CPPUNIT_ASSERT( !types.empty() );
01435     for( size_t i = 0; i < types.size(); ++i )
01436     {
01437         get_nonideal_element( types[i], pd, true );
01438         compare_analytical_and_numerical_diagonals( qm, pd );
01439     }
01440 }
01441 
01442 class PatchTranslate : public QualityMetricTester::PatchXform
01443 {
01444     Vector3D offset;
01445 
01446   public:
01447     inline PatchTranslate( Vector3D s ) : offset( s ) {}
01448     virtual ~PatchTranslate();
01449     virtual void xform( PatchData& pd, PlanarDomain* dom );
01450     virtual void xform_grad( std::vector< Vector3D >& );
01451 };
01452 
01453 PatchTranslate::~PatchTranslate() {}
01454 void PatchTranslate::xform( PatchData& pd, PlanarDomain* dom )
01455 {
01456     // If patch has associated planar geometry, translate that too
01457     MeshDomain* geom = pd.get_domain();
01458     if( geom )
01459     {
01460         PlanarDomain* plane = dynamic_cast< PlanarDomain* >( geom );
01461         CPPUNIT_ASSERT( plane );
01462         CPPUNIT_ASSERT( dom );
01463 
01464         Vector3D norm  = plane->get_normal();
01465         Vector3D point = plane->get_origin() + offset;
01466         dom->set_plane( norm, point );
01467         pd.set_domain( dom );
01468     }
01469 
01470     // Translate mesh vertices
01471     MsqError err;
01472     for( size_t i = 0; i < pd.num_nodes(); ++i )
01473         pd.move_vertex( offset, i, err );
01474 }
01475 void PatchTranslate::xform_grad( std::vector< Vector3D >& ) {}
01476 
01477 class PatchScale : public QualityMetricTester::PatchXform
01478 {
01479     double scaleFactor;
01480 
01481   public:
01482     inline PatchScale( double s ) : scaleFactor( s ) {}
01483     virtual ~PatchScale();
01484     virtual void xform( PatchData& pd, PlanarDomain* dom );
01485     virtual void xform_grad( std::vector< Vector3D >& );
01486 };
01487 
01488 PatchScale::~PatchScale() {}
01489 void PatchScale::xform( PatchData& pd, PlanarDomain* dom )
01490 {
01491     // If patch has associated planar geometry, scale distance from origin
01492     MeshDomain* geom = pd.get_domain();
01493     if( geom )
01494     {
01495         PlanarDomain* plane = dynamic_cast< PlanarDomain* >( geom );
01496         CPPUNIT_ASSERT( plane );
01497         CPPUNIT_ASSERT( dom );
01498 
01499         Vector3D norm  = plane->get_normal();
01500         Vector3D point = scaleFactor * plane->get_origin();
01501         dom->set_plane( norm, point );
01502         pd.set_domain( dom );
01503     }
01504 
01505     // Scale about origin
01506     MsqError err;
01507     for( size_t i = 0; i < pd.num_nodes(); ++i )
01508         pd.set_vertex_coordinates( pd.vertex_by_index( i ) * scaleFactor, i, err );
01509 }
01510 void PatchScale::xform_grad( std::vector< Vector3D >& ) {}
01511 
01512 class PatchRotate : public QualityMetricTester::PatchXform
01513 {
01514     Matrix3D rotation;
01515 
01516   public:
01517     inline PatchRotate( Matrix3D s ) : rotation( s ) {}
01518     virtual ~PatchRotate();
01519     virtual void xform( PatchData& pd, PlanarDomain* dom );
01520     virtual void xform_grad( std::vector< Vector3D >& grad );
01521 };
01522 
01523 PatchRotate::~PatchRotate() {}
01524 void PatchRotate::xform( PatchData& pd, PlanarDomain* dom )
01525 {
01526     // If patch has associated planar geometry, rotate that also
01527     MeshDomain* geom = pd.get_domain();
01528     if( geom )
01529     {
01530         PlanarDomain* plane = dynamic_cast< PlanarDomain* >( geom );
01531         CPPUNIT_ASSERT( plane );
01532         CPPUNIT_ASSERT( dom );
01533 
01534         Vector3D norm  = rotation * plane->get_normal();
01535         Vector3D point = rotation * plane->get_origin();
01536         dom->set_plane( norm, point );
01537         pd.set_domain( dom );
01538     }
01539 
01540     // Scale about origin
01541     MsqError err;
01542     for( size_t i = 0; i < pd.num_nodes(); ++i )
01543         pd.set_vertex_coordinates( rotation * pd.vertex_by_index( i ), i, err );
01544 }
01545 void PatchRotate::xform_grad( std::vector< Vector3D >& grad )
01546 {
01547     for( unsigned i = 0; i < grad.size(); ++i )
01548         grad[i] = rotation * grad[i];
01549 }
01550 
01551 void QualityMetricTester::test_transform_invariant( QualityMetric* qm,
01552                                                     QualityMetricTester::PatchXform& xform,
01553                                                     bool untangle )
01554 {
01555     MsqPrintError err( std::cout );
01556     PatchData pd;
01557     if( mSettings ) pd.attach_settings( mSettings );
01558     std::vector< size_t > handles;
01559     double qmval;
01560     bool rval;
01561     PlanarDomain dom( Vector3D( 0, 0, 1 ), Vector3D( 0, 0, 0 ) );
01562     size_t i, j;
01563 
01564     CPPUNIT_ASSERT( !types.empty() );
01565     for( i = 0; i < types.size(); ++i )
01566     {
01567         if( untangle )
01568             get_inverted_element( types[i], pd );
01569         else
01570             get_nonideal_element( types[i], pd );
01571         qm->get_evaluations( pd, handles, false, err );
01572         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01573         CPPUNIT_ASSERT( !handles.empty() );
01574 
01575         // get value(s) for ideal element
01576         std::vector< double > values( handles.size() );
01577         for( j = 0; j < handles.size(); ++j )
01578         {
01579             rval = qm->evaluate( pd, handles[j], values[j], err );
01580             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01581             CPPUNIT_ASSERT( rval );
01582         }
01583 
01584         // compare values for transformed patch
01585         xform.xform( pd, &dom );
01586         for( j = 0; j < handles.size(); ++j )
01587         {
01588             rval = qm->evaluate( pd, handles[j], qmval, err );
01589             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01590             CPPUNIT_ASSERT( rval );
01591             CPPUNIT_ASSERT_DOUBLES_EQUAL( values[j], qmval, 1e-3 );
01592         }
01593     }
01594 }
01595 
01596 void QualityMetricTester::test_grad_transform_invariant( QualityMetric* qm,
01597                                                          QualityMetricTester::PatchXform& xform,
01598                                                          bool untangle )
01599 {
01600     MsqPrintError err( std::cout );
01601     PatchData pd;
01602     if( mSettings ) pd.attach_settings( mSettings );
01603     std::vector< size_t > handles, indices2;
01604     std::vector< Vector3D > grad2;
01605     double qmval;
01606     bool rval;
01607     PlanarDomain dom( Vector3D( 0, 0, 1 ), Vector3D( 0, 0, 0 ) );
01608     size_t i, j, k;
01609 
01610     CPPUNIT_ASSERT( !types.empty() );
01611     for( i = 0; i < types.size(); ++i )
01612     {
01613         if( untangle )
01614             get_inverted_element( types[i], pd );
01615         else
01616             get_nonideal_element( types[i], pd );
01617         qm->get_evaluations( pd, handles, false, err );
01618         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01619         CPPUNIT_ASSERT( !handles.empty() );
01620 
01621         // get value(s) for orignal patch
01622         std::vector< std::vector< Vector3D > > grad1( handles.size() );
01623         std::vector< std::vector< size_t > > indices1( handles.size() );
01624         for( j = 0; j < handles.size(); ++j )
01625         {
01626             rval = qm->evaluate_with_gradient( pd, handles[j], qmval, indices1[j], grad1[j], err );
01627             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01628             CPPUNIT_ASSERT( rval );
01629             CPPUNIT_ASSERT( !grad1[j].empty() );
01630         }
01631 
01632         // transform patch
01633         xform.xform( pd, &dom );
01634 
01635         // compare values
01636         for( j = 0; j < handles.size(); ++j )
01637         {
01638             rval = qm->evaluate_with_gradient( pd, handles[j], qmval, indices2, grad2, err );
01639             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01640             CPPUNIT_ASSERT( rval );
01641 
01642             CPPUNIT_ASSERT_EQUAL( indices1[j].size(), indices2.size() );
01643             CPPUNIT_ASSERT_EQUAL( grad1[j].size(), grad2.size() );
01644             CPPUNIT_ASSERT( indices1[j] == indices2 );
01645 
01646             xform.xform_grad( grad1[j] );
01647             for( k = 0; k < grad2.size(); ++k )
01648                 CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[j][k], grad2[k], 1e-3 );
01649         }
01650     }
01651 }
01652 
01653 void QualityMetricTester::test_hessian_transform_invariant( QualityMetric* qm,
01654                                                             QualityMetricTester::PatchXform& xform,
01655                                                             bool untangle )
01656 {
01657     MsqPrintError err( std::cout );
01658     PatchData pd;
01659     if( mSettings ) pd.attach_settings( mSettings );
01660     std::vector< size_t > handles, indices2;
01661     std::vector< Vector3D > grad1, grad2;
01662     std::vector< Matrix3D > hess2;
01663     double qmval;
01664     bool rval;
01665     PlanarDomain dom( Vector3D( 0, 0, 1 ), Vector3D( 0, 0, 0 ) );
01666     size_t i, j, k;
01667 
01668     CPPUNIT_ASSERT( !types.empty() );
01669     for( i = 0; i < types.size(); ++i )
01670     {
01671         if( untangle )
01672             get_inverted_element( types[i], pd );
01673         else
01674             get_nonideal_element( types[i], pd );
01675         qm->get_evaluations( pd, handles, false, err );
01676         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01677         CPPUNIT_ASSERT( !handles.empty() );
01678 
01679         std::vector< std::vector< size_t > > indices1( handles.size() );
01680         std::vector< std::vector< Matrix3D > > hess1( handles.size() );
01681         for( j = 0; j < handles.size(); ++j )
01682         {
01683             rval = qm->evaluate_with_Hessian( pd, handles[j], qmval, indices1[j], grad1, hess1[j], err );
01684             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01685             CPPUNIT_ASSERT( rval );
01686             CPPUNIT_ASSERT( !hess1[j].empty() );
01687         }
01688 
01689         xform.xform( pd, &dom );
01690 
01691         for( j = 0; j < handles.size(); ++j )
01692         {
01693             rval = qm->evaluate_with_Hessian( pd, handles[j], qmval, indices2, grad2, hess2, err );
01694             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01695             CPPUNIT_ASSERT( rval );
01696 
01697             CPPUNIT_ASSERT_EQUAL( indices1[j].size(), indices2.size() );
01698             CPPUNIT_ASSERT_EQUAL( hess1[j].size(), hess2.size() );
01699             CPPUNIT_ASSERT( indices1[j] == indices2 );
01700 
01701             for( k = 0; k < hess2.size(); ++k )
01702                 CPPUNIT_ASSERT_MATRICES_EQUAL( hess1[j][k], hess2[k], 1e-3 );
01703         }
01704     }
01705 }
01706 
01707 void QualityMetricTester::test_location_invariant( QualityMetric* qm, bool untangle )
01708 {
01709     PatchTranslate xform( Vector3D( 1, 1, -1 ) );
01710     test_transform_invariant( qm, xform, untangle );
01711 }
01712 
01713 void QualityMetricTester::test_scale_invariant( QualityMetric* qm, bool untangle )
01714 {
01715     PatchScale scale1( 0.5 ), scale2( 2.0 );
01716     test_transform_invariant( qm, scale1, untangle );
01717     test_transform_invariant( qm, scale2, untangle );
01718 }
01719 
01720 void QualityMetricTester::test_orient_invariant( QualityMetric* qm, bool untangle )
01721 {
01722     // rotate 45 degrees about x-axis
01723     const double f = std::sqrt( 2.0 ) / 2;
01724     Matrix3D rotation( 1, 0, 0, 0, f, f, 0, -f, f );
01725     PatchRotate xform( rotation );
01726     test_transform_invariant( qm, xform, untangle );
01727 }
01728 
01729 void QualityMetricTester::test_grad_location_invariant( QualityMetric* qm, bool untangle )
01730 {
01731     PatchTranslate xform( Vector3D( 1, 1, -1 ) );
01732     test_grad_transform_invariant( qm, xform, untangle );
01733 }
01734 
01735 void QualityMetricTester::test_hessian_location_invariant( QualityMetric* qm, bool untangle )
01736 {
01737     PatchTranslate xform( Vector3D( 1, 1, -1 ) );
01738     test_hessian_transform_invariant( qm, xform, untangle );
01739 }
01740 
01741 void QualityMetricTester::test_grad_orient_invariant( QualityMetric* qm, bool untangle )
01742 {
01743     // rotate 45 degrees about x-axis
01744     const double f = std::sqrt( 2.0 ) / 2;
01745     Matrix3D rotation( 1, 0, 0, 0, f, f, 0, -f, f );
01746     PatchRotate xform( rotation );
01747     test_grad_transform_invariant( qm, xform, untangle );
01748 }
01749 
01750 void QualityMetricTester::test_measures_transform( QualityMetric* qm,
01751                                                    QualityMetricTester::PatchXform& xform,
01752                                                    bool unit_area )
01753 {
01754     MsqPrintError err( std::cout );
01755     PatchData pd;
01756     if( mSettings ) pd.attach_settings( mSettings );
01757     std::vector< size_t > handles;
01758     double qmval1, qmval2;
01759     bool rval;
01760     PlanarDomain dom( Vector3D( 0, 0, 1 ), Vector3D( 0, 0, 0 ) );
01761     bool decreasing = false;
01762     if( qm->get_negate_flag() != 1 )
01763     {
01764         CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
01765         decreasing = true;
01766     }
01767 
01768     CPPUNIT_ASSERT( !types.empty() );
01769     for( size_t i = 0; i < types.size(); ++i )
01770     {
01771         get_ideal_element( types[i], unit_area, pd );
01772         qm->get_evaluations( pd, handles, false, err );
01773         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01774         CPPUNIT_ASSERT( !handles.empty() );
01775 
01776         rval = qm->evaluate( pd, handles[0], qmval1, err );
01777         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01778         CPPUNIT_ASSERT( rval );
01779 
01780         xform.xform( pd, &dom );
01781         rval = qm->evaluate( pd, handles[0], qmval2, err );
01782         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01783         CPPUNIT_ASSERT( rval );
01784         if( decreasing )
01785             CPPUNIT_ASSERT( qmval2 < qmval1 );
01786         else
01787             CPPUNIT_ASSERT( qmval2 > qmval1 );
01788     }
01789 }
01790 
01791 void QualityMetricTester::test_measures_size( QualityMetric* qm, bool unit_area )
01792 {
01793     PatchScale s1( 0.5 ), s2( 2.0 );
01794     test_measures_transform( qm, s1, unit_area );
01795     test_measures_transform( qm, s2, unit_area );
01796 }
01797 
01798 void QualityMetricTester::test_measures_in_plane_orientation( QualityMetric* qm )
01799 {
01800     // rotate 45 degrees about z-axis
01801     const double f = std::sqrt( 2.0 ) / 2;
01802     Matrix3D rotation( f, -f, 0, f, f, 0, 0, 0, 1 );
01803     PatchRotate xform( rotation );
01804     test_measures_transform( qm, xform, false );
01805 }
01806 
01807 void QualityMetricTester::test_measures_out_of_plane_orientation( QualityMetric* qm )
01808 {
01809     // rotate 45 degrees about z-axis
01810     const double f = std::sqrt( 2.0 ) / 2;
01811     Matrix3D rotation( 1, 0, 0, 0, f, f, 0, -f, f );
01812     PatchRotate xform( rotation );
01813     test_measures_transform( qm, xform, false );
01814 }
01815 
01816 typedef void ( QualityMetricTester::*patch_func_t )( EntityTopology, PatchData& );
01817 static void test_bad_element( QualityMetricTester* instance,
01818                               QualityMetric* qm,
01819                               patch_func_t type_func,
01820                               bool should_succeed,
01821                               const std::vector< EntityTopology >& types,
01822                               const Settings* mSettings )
01823 {
01824     MsqPrintError err( std::cout );
01825     PatchData pd;
01826     if( mSettings ) pd.attach_settings( mSettings );
01827     std::vector< size_t > handles;
01828     double qmval;
01829     bool rval;
01830 
01831     CPPUNIT_ASSERT( !types.empty() );
01832     for( size_t i = 0; i < types.size(); ++i )
01833     {
01834         ( instance->*type_func )( types[i], pd );
01835         qm->get_evaluations( pd, handles, false, err );
01836         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01837         CPPUNIT_ASSERT( !handles.empty() );
01838 
01839         rval = true;
01840         for( size_t j = 0; j < handles.size(); ++j )
01841         {
01842             rval = rval && qm->evaluate( pd, handles[j], qmval, err );
01843             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01844         }
01845         CPPUNIT_ASSERT_EQUAL( should_succeed, rval );
01846     }
01847 }
01848 
01849 void QualityMetricTester::test_evaluate_inverted_element( QualityMetric* qm, bool should_succeed )
01850 {
01851     test_bad_element( this, qm, &QualityMetricTester::get_inverted_element, should_succeed, types, mSettings );
01852 }
01853 void QualityMetricTester::test_evaluate_degenerate_element( QualityMetric* qm, bool should_succeed )
01854 {
01855     test_bad_element( this, qm, &QualityMetricTester::get_degenerate_element, should_succeed, types, mSettings );
01856 }
01857 void QualityMetricTester::test_evaluate_zero_element( QualityMetric* qm, bool should_succeed )
01858 {
01859     test_bad_element( this, qm, &QualityMetricTester::get_zero_element, should_succeed, types, mSettings );
01860 }
01861 
01862 void QualityMetricTester::test_measures_quality( QualityMetric* qm )
01863 {
01864     MsqPrintError err( std::cout );
01865     PatchData pd;
01866     if( mSettings ) pd.attach_settings( mSettings );
01867     std::vector< size_t > handles;
01868     double qmval1, qmval2;
01869     bool rval;
01870     bool decreasing = false;
01871     if( qm->get_negate_flag() != 1 )
01872     {
01873         CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
01874         decreasing = true;
01875     }
01876 
01877     CPPUNIT_ASSERT( !types.empty() );
01878     for( size_t i = 0; i < types.size(); ++i )
01879     {
01880         // get a patch
01881         get_ideal_element( types[i], false, pd );
01882         qm->get_evaluations( pd, handles, false, err );
01883         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01884         CPPUNIT_ASSERT( !handles.empty() );
01885         // evaluate for ideal element
01886         rval = qm->evaluate( pd, handles[0], qmval1, err );
01887         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01888         CPPUNIT_ASSERT( rval );
01889         // calculate element centroid
01890         const size_t* verts     = pd.element_by_index( 0 ).get_vertex_index_array();
01891         const MsqVertex* coords = pd.get_vertex_array( err );
01892         size_t n                = pd.element_by_index( 0 ).vertex_count();
01893         Vector3D cent( 0, 0, 0 );
01894         Vector3D coords0 = coords[verts[0]];
01895         for( unsigned j = 0; j < n; ++j )
01896             cent += coords[verts[j]];
01897         cent /= n;
01898         // evaluate for non-ideal element (decreased area)
01899         pd.set_vertex_coordinates( 0.5 * ( cent + coords0 ), verts[0], err );
01900         rval = qm->evaluate( pd, handles[0], qmval2, err );
01901         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01902         CPPUNIT_ASSERT( rval );
01903         // check values
01904         if( decreasing )
01905         {
01906             CPPUNIT_ASSERT( qmval2 < qmval1 );
01907         }
01908         else
01909         {
01910             CPPUNIT_ASSERT( qmval2 > qmval1 );
01911         }
01912         // evaluate for non-ideal element (increased area)
01913         pd.set_vertex_coordinates( 3 * coords0 - 2 * cent, verts[0], err );
01914         rval = qm->evaluate( pd, handles[0], qmval2, err );
01915         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01916         CPPUNIT_ASSERT( rval );
01917         // check values
01918         if( decreasing )
01919         {
01920             CPPUNIT_ASSERT( qmval2 < qmval1 );
01921         }
01922         else
01923         {
01924             CPPUNIT_ASSERT( qmval2 > qmval1 );
01925         }
01926     }
01927 }
01928 
01929 void QualityMetricTester::test_domain_deviation_quality( QualityMetric* qm )
01930 {
01931     MsqPrintError err( std::cout );
01932     PatchData pd;
01933     if( mSettings ) pd.attach_settings( mSettings );
01934     std::vector< size_t > handles, indices;
01935     size_t i;
01936     double qmval1, qmval2;
01937     bool rval;
01938     bool decreasing = false;
01939     if( qm->get_negate_flag() != 1 )
01940     {
01941         CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
01942         decreasing = true;
01943     }
01944 
01945     CPPUNIT_ASSERT( type_is_supported( TRIANGLE ) || type_is_supported( QUADRILATERAL ) );
01946 
01947     if( type_is_supported( TRIANGLE ) )
01948     {
01949         // get a patch
01950         get_ideal_tris( pd, false );
01951         qm->get_evaluations( pd, handles, false, err );
01952         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01953         CPPUNIT_ASSERT( !handles.empty() );
01954         // find an evaluation that depends on the free vertex
01955         for( i = 0; i < handles.size(); ++i )
01956         {
01957             rval = qm->evaluate_with_indices( pd, handles[i], qmval1, indices, err );
01958             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01959             CPPUNIT_ASSERT( rval );
01960             if( !indices.empty() ) break;
01961         }
01962         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );   // one free vertex in patch
01963         CPPUNIT_ASSERT_EQUAL( (size_t)0, indices.front() );  // one free vertex in patch
01964         const size_t handle = handles[i];
01965         // rotate patch about origin such that elements are no longer in the plane
01966         const Matrix3D rot( M_PI / 4.0, Vector3D( 1, 1, 0 ) );
01967         for( i = 0; i < pd.num_nodes(); ++i )
01968             pd.set_vertex_coordinates( rot * pd.vertex_by_index( i ), i, err );
01969         // evaluate for rotated patch
01970         rval = qm->evaluate( pd, handle, qmval2, err );
01971         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01972         CPPUNIT_ASSERT( rval );
01973         // check quality
01974         if( decreasing )
01975         {
01976             CPPUNIT_ASSERT( qmval1 > qmval2 );
01977         }
01978         else
01979         {
01980             CPPUNIT_ASSERT( qmval1 < qmval2 );
01981         }
01982     }
01983 
01984     if( type_is_supported( QUADRILATERAL ) )
01985     {
01986         // get a patch
01987         get_ideal_quads( pd );
01988         qm->get_evaluations( pd, handles, false, err );
01989         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01990         CPPUNIT_ASSERT( !handles.empty() );
01991         // find an evaluation that depends on the free vertex
01992         for( i = 0; i < handles.size(); ++i )
01993         {
01994             rval = qm->evaluate_with_indices( pd, handles[i], qmval1, indices, err );
01995             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01996             CPPUNIT_ASSERT( rval );
01997             if( !indices.empty() ) break;
01998         }
01999         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );   // one free vertex in patch
02000         CPPUNIT_ASSERT_EQUAL( (size_t)0, indices.front() );  // one free vertex in patch
02001         const size_t handle = handles[i];
02002         // rotate patch about origin such that elements are no longer in the plane
02003         const Matrix3D rot( M_PI / 3.0, Vector3D( 1, 0, 0 ) );
02004         for( i = 0; i < pd.num_nodes(); ++i )
02005             pd.set_vertex_coordinates( rot * pd.vertex_by_index( i ), i, err );
02006         // evaluate for rotated patch
02007         rval = qm->evaluate( pd, handle, qmval2, err );
02008         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02009         CPPUNIT_ASSERT( rval );
02010         // check quality
02011         if( decreasing )
02012         {
02013             CPPUNIT_ASSERT( qmval1 > qmval2 );
02014         }
02015         else
02016         {
02017             CPPUNIT_ASSERT( qmval1 < qmval2 );
02018         }
02019     }
02020 }
02021 
02022 void QualityMetricTester::test_domain_deviation_gradient( QualityMetric* qm )
02023 {
02024     MsqPrintError err( std::cout );
02025     PatchData pd;
02026     if( mSettings ) pd.attach_settings( mSettings );
02027     std::vector< size_t > handles, indices;
02028     std::vector< Vector3D > grad;
02029     size_t i;
02030     double qmval;
02031     bool rval;
02032 
02033     CPPUNIT_ASSERT( type_is_supported( TRIANGLE ) || type_is_supported( QUADRILATERAL ) );
02034 
02035     if( type_is_supported( TRIANGLE ) )
02036     {
02037         // get a patch
02038         get_ideal_tris( pd, false );
02039         qm->get_evaluations( pd, handles, false, err );
02040         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02041         CPPUNIT_ASSERT( !handles.empty() );
02042         // move the free vertex out of the plane
02043         pd.move_vertex( Vector3D( 0, 0, 0.2 ), 0, err );
02044         // find an evaluation that depends on the free vertex
02045         for( i = 0; i < handles.size(); ++i )
02046         {
02047             rval = qm->evaluate_with_gradient( pd, handles[i], qmval, indices, grad, err );
02048             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02049             CPPUNIT_ASSERT( rval );
02050             if( !grad.empty() ) break;
02051         }
02052         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );   // one free vertex in patch
02053         CPPUNIT_ASSERT_EQUAL( (size_t)0, indices.front() );  // one free vertex in patch
02054         CPPUNIT_ASSERT_EQUAL( (size_t)1, grad.size() );
02055 
02056         Vector3D v = qm->get_negate_flag() * Vector3D( 0, 0, 1 );
02057         CPPUNIT_ASSERT( v % grad[0] > DBL_EPSILON );
02058     }
02059 
02060     if( type_is_supported( QUADRILATERAL ) )
02061     {
02062         // get a patch
02063         get_ideal_quads( pd );
02064         qm->get_evaluations( pd, handles, false, err );
02065         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02066         CPPUNIT_ASSERT( !handles.empty() );
02067         // move the free vertex out of the plane
02068         pd.move_vertex( Vector3D( 0, 0, 0.2 ), 0, err );
02069         // find an evaluation that depends on the free vertex
02070         for( i = 0; i < handles.size(); ++i )
02071         {
02072             rval = qm->evaluate_with_gradient( pd, handles[i], qmval, indices, grad, err );
02073             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02074             CPPUNIT_ASSERT( rval );
02075             if( !grad.empty() ) break;
02076         }
02077         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );   // one free vertex in patch
02078         CPPUNIT_ASSERT_EQUAL( (size_t)0, indices.front() );  // one free vertex in patch
02079         CPPUNIT_ASSERT_EQUAL( (size_t)1, grad.size() );
02080 
02081         Vector3D v = qm->get_negate_flag() * Vector3D( 0, 0, 1 );
02082         CPPUNIT_ASSERT( v % grad[0] > DBL_EPSILON );
02083     }
02084 }
02085 
02086 void QualityMetricTester::test_gradient_reflects_quality( QualityMetric* qm )
02087 {
02088     MsqPrintError err( std::cout );
02089     PatchData pd;
02090     if( mSettings ) pd.attach_settings( mSettings );
02091     std::vector< size_t > handles, indices;
02092     std::vector< Vector3D > grad;
02093     double qmval;
02094     bool rval;
02095     bool decreasing = false;
02096     if( qm->get_negate_flag() != 1 )
02097     {
02098         CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
02099         decreasing = true;
02100     }
02101 
02102     CPPUNIT_ASSERT( !types.empty() );
02103     for( size_t i = 0; i < types.size(); ++i )
02104     {
02105         // get a patch
02106         get_ideal_element( types[i], false, pd );
02107         qm->get_evaluations( pd, handles, false, err );
02108         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02109         CPPUNIT_ASSERT( !handles.empty() );
02110         // make element non-ideal
02111         const size_t* verts     = pd.element_by_index( 0 ).get_vertex_index_array();
02112         const MsqVertex* coords = pd.get_vertex_array( err );
02113         size_t n                = pd.element_by_index( 0 ).vertex_count();
02114         Vector3D sum( 0, 0, 0 );
02115         for( unsigned j = 0; j < n; ++j )
02116             sum += coords[verts[j]];
02117         const Vector3D init_pos = coords[verts[0]];
02118         const Vector3D new_pos  = 0.5 * coords[verts[0]] + sum / ( 2 * n );
02119         pd.set_vertex_coordinates( new_pos, verts[0], err );
02120         // evaluate for non-ideal element
02121         rval = qm->evaluate_with_gradient( pd, handles[0], qmval, indices, grad, err );
02122         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02123         CPPUNIT_ASSERT( rval );
02124         // find gradient for vertex 0
02125         std::vector< size_t >::iterator it = std::find( indices.begin(), indices.end(), verts[0] );
02126         CPPUNIT_ASSERT( it != indices.end() );
02127         Vector3D g = grad[it - indices.begin()];
02128         // check gradient direction
02129         Vector3D v     = init_pos - coords[verts[0]];
02130         double dotprod = v % g;
02131         if( decreasing )
02132             CPPUNIT_ASSERT( dotprod > 1e-6 );
02133         else
02134             CPPUNIT_ASSERT( dotprod < -1e-6 );
02135     }
02136 }
02137 
02138 void QualityMetricTester::test_measures_vertex_quality( QualityMetric* qm )
02139 {
02140     MsqPrintError err( std::cout );
02141     PatchData pd;
02142     if( mSettings ) pd.attach_settings( mSettings );
02143     double qmval1, qmval2;
02144     bool rval;
02145     bool increasing = true;
02146     if( qm->get_negate_flag() != 1 )
02147     {
02148         CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
02149         increasing = true;
02150     }
02151 
02152     unsigned count = 0;
02153     for( size_t i = 0; i < 3; ++i )
02154     {
02155         // get a patch
02156         switch( i )
02157         {
02158             case 0:
02159                 if( !type_is_supported( TRIANGLE ) ) continue;
02160                 get_ideal_tris( pd, false );
02161                 break;
02162             case 1:
02163                 if( !type_is_supported( QUADRILATERAL ) ) continue;
02164                 get_ideal_quads( pd );
02165                 break;
02166             case 2:
02167                 if( !type_is_supported( HEXAHEDRON ) ) continue;
02168                 get_ideal_hexes( pd );
02169                 break;
02170         }
02171         ++count;
02172 
02173         // evaluate for ideal element
02174         rval = qm->evaluate( pd, 0, qmval1, err );
02175         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02176         CPPUNIT_ASSERT( rval );
02177         // move center vertex
02178         pd.move_vertex( Vector3D( 0.3, 0.3, 0.3 ), 0, err );
02179         // evaluate for non-ideal element
02180         rval = qm->evaluate( pd, 0, qmval2, err );
02181         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02182         CPPUNIT_ASSERT( rval );
02183         // check values
02184         if( increasing )
02185         {
02186             CPPUNIT_ASSERT( qmval2 > qmval1 );
02187         }
02188         else
02189         {
02190             CPPUNIT_ASSERT( qmval2 < qmval1 );
02191         }
02192     }
02193     // make sure we tested something
02194     CPPUNIT_ASSERT( count > 0 );
02195 }
02196 
02197 void QualityMetricTester::test_vertex_gradient_reflects_quality( QualityMetric* qm )
02198 {
02199     MsqPrintError err( std::cout );
02200     PatchData pd;
02201     if( mSettings ) pd.attach_settings( mSettings );
02202     std::vector< size_t > indices;
02203     std::vector< Vector3D > grad;
02204     double qmval;
02205     bool rval;
02206     bool increasing = true;
02207     if( qm->get_negate_flag() != 1 )
02208     {
02209         CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
02210         increasing = false;
02211     }
02212 
02213     unsigned count = 0;
02214     for( size_t i = 0; i < 3; ++i )
02215     {
02216         // get a patch
02217         switch( i )
02218         {
02219             case 0:
02220                 if( !type_is_supported( TRIANGLE ) ) continue;
02221                 get_ideal_tris( pd, false );
02222                 break;
02223             case 1:
02224                 if( !type_is_supported( QUADRILATERAL ) ) continue;
02225                 get_ideal_quads( pd );
02226                 break;
02227             case 2:
02228                 if( !type_is_supported( HEXAHEDRON ) ) continue;
02229                 get_ideal_hexes( pd );
02230                 break;
02231         }
02232 
02233         // move center vertex
02234         pd.move_vertex( Vector3D( 0.3, 0.3, 0.3 ), 0, err );
02235 
02236         // evaluate for ideal non-element
02237         rval = qm->evaluate_with_gradient( pd, 0, qmval, indices, grad, err );
02238         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02239         CPPUNIT_ASSERT( rval );
02240         // find gradient for vertex 0
02241         if( indices.empty() ) continue;
02242         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );  // only 1 free vertex
02243         CPPUNIT_ASSERT_EQUAL( (size_t)0, indices[0] );      // and that vertex is vertex 0
02244         ++count;
02245         Vector3D g = grad[0];
02246         // check gradient direction
02247         if( increasing )
02248             CPPUNIT_ASSERT( g[2] > 1e-6 );
02249         else
02250             CPPUNIT_ASSERT( g[2] < -1e-6 );
02251     }
02252     // make sure we tested something
02253     CPPUNIT_ASSERT( count > 0 );
02254 }
02255 
02256 void QualityMetricTester::test_ideal_element_zero_gradient( QualityMetric* qm, bool unit_area )
02257 {
02258     MsqPrintError err( std::cout );
02259     PatchData pd;
02260     if( mSettings ) pd.attach_settings( mSettings );
02261     std::vector< size_t > handles, indices;
02262     std::vector< Vector3D > grad;
02263     double qmval;
02264     bool rval;
02265 
02266     CPPUNIT_ASSERT( !types.empty() );
02267     for( size_t i = 0; i < types.size(); ++i )
02268     {
02269         // get a patch
02270         get_ideal_element( types[i], unit_area, pd );
02271         qm->get_evaluations( pd, handles, false, err );
02272         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02273         CPPUNIT_ASSERT( !handles.empty() );
02274         for( size_t j = 0; j < handles.size(); ++j )
02275         {
02276             // evaluate
02277             rval = qm->evaluate_with_gradient( pd, handles[j], qmval, indices, grad, err );
02278             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02279             CPPUNIT_ASSERT( rval );
02280             CPPUNIT_ASSERT_EQUAL( indices.size(), grad.size() );
02281             // check that all gradients are zero
02282             for( size_t k = 0; k < grad.size(); ++k )
02283             {
02284                 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, 0, 0 ), grad[k], 1e-4 );
02285             }  // for(k < grad.size())
02286         }      // for(j < handles.size())
02287     }          // for(i < types.size())
02288 }
02289 
02290 void QualityMetricTester::test_ideal_element_zero_vertex_gradient( QualityMetric* qm, bool unit_area )
02291 {
02292     MsqPrintError err( std::cout );
02293     PatchData pd;
02294     if( mSettings ) pd.attach_settings( mSettings );
02295     std::vector< size_t > handles, indices;
02296     std::vector< Vector3D > grad;
02297     double qmval;
02298     bool rval;
02299 
02300     unsigned count = 0;
02301     for( size_t i = 0; i < 3; ++i )
02302     {
02303         // get a patch
02304         switch( i )
02305         {
02306             case 0:
02307                 if( !type_is_supported( TRIANGLE ) ) continue;
02308                 get_ideal_tris( pd, unit_area );
02309                 break;
02310             case 1:
02311                 if( !type_is_supported( QUADRILATERAL ) ) continue;
02312                 get_ideal_quads( pd );
02313                 break;
02314             case 2:
02315                 if( !type_is_supported( HEXAHEDRON ) ) continue;
02316                 get_ideal_hexes( pd );
02317                 break;
02318         }
02319 
02320         qm->get_evaluations( pd, handles, false, err );
02321         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02322         CPPUNIT_ASSERT( !handles.empty() );
02323         for( size_t j = 0; j < handles.size(); ++j )
02324         {
02325             // evaluate
02326             rval = qm->evaluate_with_gradient( pd, handles[j], qmval, indices, grad, err );
02327             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02328             CPPUNIT_ASSERT( rval );
02329             CPPUNIT_ASSERT_EQUAL( indices.size(), grad.size() );
02330             if( grad.empty() ) continue;
02331             ++count;
02332             CPPUNIT_ASSERT_EQUAL( (size_t)1, grad.size() );  // only 1 free vertex in patch
02333                                                              // check that all gradients are zero
02334             CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, 0, 0 ), grad[0], 1e-4 );
02335         }  // for(j < handles.size())
02336     }      // for(i < 3)
02337            // make sure we tested something
02338     CPPUNIT_ASSERT( count > 0 );
02339 }
02340 
02341 static double** allocate_matrix( unsigned n )
02342 {
02343     unsigned num_ptr = n;
02344     if( num_ptr % 2 ) ++num_ptr;
02345 
02346     void* storage = malloc( n * n * sizeof( double ) + num_ptr * sizeof( double* ) );
02347     double** ptrs = (double**)storage;
02348     double* data  = (double*)( ptrs + num_ptr );
02349     for( unsigned i = 0; i < n; ++i )
02350         ptrs[i] = data + i * n;
02351     return ptrs;
02352 }
02353 
02354 static double value( const Matrix3D* blocks, unsigned n, unsigned i, unsigned j )
02355 {
02356     if( i > j ) std::swap( i, j );
02357 
02358     unsigned mi  = i / 3;
02359     unsigned mj  = j / 3;
02360     unsigned idx = n * mi - mi * ( mi + 1 ) / 2 + mj;
02361     return blocks[idx][i % 3][j % 3];
02362 }
02363 /*
02364 static bool positive_definite( const Matrix3D* blocks, unsigned n )
02365 {
02366     // Do Cholesky-Banachiewicz decompositon of the matrix.
02367     // If this results in any imaginary values for diagonal terms,
02368     // then the matrix is not positive-definite.
02369   const int N = 3*n;
02370   double** mat = allocate_matrix(N);
02371   int i, j, k;
02372   double s;
02373   for (i = 0; i < N; ++i)
02374   {
02375     for (j = 0; j < i; ++j)
02376     {
02377       s = value(blocks,n,i,j);
02378       for (k = 0; k < j - 1; ++k)
02379         s -= mat[i][k]*mat[j][k];
02380       mat[i][j] = s / mat[j][j];
02381     }
02382     s = value(blocks,n,i,i);
02383     for (k = 0; k < i - 1; ++k)
02384       s -= mat[i][k]*mat[i][k];
02385     if (s < 0.0)
02386     {
02387       free( mat );
02388       return false;
02389     }
02390     mat[i][i] = sqrt(s);
02391   }
02392   return true;
02393 }
02394 */
02395 
02396 /** Test that Hessian is positive-definite for ideal elements */
02397 void QualityMetricTester::test_ideal_element_positive_definite_Hessian( QualityMetric* qm, bool unit_area )
02398 {
02399     MsqPrintError err( std::cout );
02400     PatchData pd;
02401     if( mSettings ) pd.attach_settings( mSettings );
02402     std::vector< size_t > handles, indices;
02403     std::vector< Vector3D > grad;
02404     std::vector< Matrix3D > Hess;
02405     double qmval;
02406     bool rval;
02407 
02408     CPPUNIT_ASSERT( !types.empty() );
02409     for( size_t i = 0; i < types.size(); ++i )
02410     {
02411         // get a patch
02412         get_ideal_element( types[i], unit_area, pd );
02413         qm->get_evaluations( pd, handles, false, err );
02414         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02415         CPPUNIT_ASSERT( !handles.empty() );
02416         for( size_t j = 0; j < handles.size(); ++j )
02417         {
02418             // evaluate
02419             rval = qm->evaluate_with_Hessian( pd, handles[j], qmval, indices, grad, Hess, err );
02420             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02421             CPPUNIT_ASSERT( rval );
02422             const size_t n = indices.size();
02423             CPPUNIT_ASSERT_EQUAL( n, grad.size() );
02424             CPPUNIT_ASSERT_EQUAL( n * ( n + 1 ) / 2, Hess.size() );
02425             // if necessary, adjust Hessian for negate_flag();
02426             if( qm->get_negate_flag() == -1 )
02427                 for( size_t k = 0; k < Hess.size(); ++k )
02428                     Hess[k] *= -1.0;
02429             // check that all diagonal blocks are positive definite
02430             size_t h_idx = 0;
02431             for( size_t k = 0; k < n; h_idx += ( n - k ), ++k )
02432             {
02433                 CPPUNIT_ASSERT( Hess[h_idx].positive_definite() );
02434             }
02435             // check that the entire Hessian is positive definite
02436             // CPPUNIT_ASSERT( positive_definite( arrptr(Hess), n ) );
02437         }  // for(j < handles.size())
02438     }      // for(i < types.size())
02439 }
02440 
02441 /** Test that diagonal bocks of Hessian are symmetric */
02442 void QualityMetricTester::test_symmetric_Hessian_diagonal_blocks( QualityMetric* qm )
02443 {
02444     MsqPrintError err( std::cout );
02445     PatchData pd;
02446     if( mSettings ) pd.attach_settings( mSettings );
02447     std::vector< size_t > handles, indices;
02448     std::vector< Vector3D > grad;
02449     std::vector< Matrix3D > Hess;
02450     double qmval;
02451     bool rval;
02452 
02453     CPPUNIT_ASSERT( !types.empty() );
02454     for( size_t i = 0; i < types.size(); ++i )
02455     {
02456         // get a patch
02457         get_ideal_element( types[i], false, pd );
02458         qm->get_evaluations( pd, handles, false, err );
02459         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02460         CPPUNIT_ASSERT( !handles.empty() );
02461         for( size_t j = 0; j < handles.size(); ++j )
02462         {
02463             // evaluate
02464             rval = qm->evaluate_with_Hessian( pd, handles[j], qmval, indices, grad, Hess, err );
02465             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
02466             CPPUNIT_ASSERT( rval );
02467             const size_t n = indices.size();
02468             CPPUNIT_ASSERT_EQUAL( n, grad.size() );
02469             CPPUNIT_ASSERT_EQUAL( n * ( n + 1 ) / 2, Hess.size() );
02470             // check that diagonal Hessian blocks are symmetric
02471             size_t h_idx = 0;
02472             for( size_t k = 0; k < n; h_idx += ( n - k ), ++k )
02473             {
02474                 CPPUNIT_ASSERT_DOUBLES_EQUAL( Hess[h_idx][0][1], Hess[h_idx][1][0], 5e-3 );
02475                 CPPUNIT_ASSERT_DOUBLES_EQUAL( Hess[h_idx][0][2], Hess[h_idx][2][0], 5e-3 );
02476                 CPPUNIT_ASSERT_DOUBLES_EQUAL( Hess[h_idx][1][2], Hess[h_idx][2][1], 5e-3 );
02477             }  // for(k < Hess.size())
02478         }      // for(j < handles.size())
02479     }          // for(i < types.size())
02480 }
02481 
02482 void QualityMetricTester::test_gradient_with_fixed_vertex( QualityMetric* qm, const Settings* settings )
02483 {
02484     CPPUNIT_ASSERT( !types.empty() );
02485     for( unsigned i = 0; i < types.size(); ++i )
02486         test_gradient_with_fixed_vertex( types[i], qm, settings );
02487 }
02488 
02489 void QualityMetricTester::test_hessian_with_fixed_vertex( QualityMetric* qm, const Settings* settings )
02490 {
02491     CPPUNIT_ASSERT( !types.empty() );
02492     for( unsigned i = 0; i < types.size(); ++i )
02493         test_hessian_with_fixed_vertex( types[i], qm, settings );
02494 }
02495 
02496 void QualityMetricTester::test_diagonal_with_fixed_vertex( QualityMetric* qm, const Settings* settings )
02497 {
02498     CPPUNIT_ASSERT( !types.empty() );
02499     for( unsigned i = 0; i < types.size(); ++i )
02500         test_diagonal_with_fixed_vertex( types[i], qm, settings );
02501 }
02502 
02503 void QualityMetricTester::test_gradient_with_fixed_vertex( EntityTopology type,
02504                                                            QualityMetric* qm,
02505                                                            const Settings* settings )
02506 {
02507     MsqPrintError err( std::cout );
02508     PatchData pd;
02509     if( settings )
02510         pd.attach_settings( settings );
02511     else if( mSettings )
02512         pd.attach_settings( mSettings );
02513 
02514     std::vector< size_t > handles, indices, indices_fixed, conn;
02515     std::vector< Vector3D > grad, grad_fixed;
02516     double val, val_fixed;
02517     bool rval;
02518     const int n = TopologyInfo::corners( type );
02519 
02520     get_nonideal_element( type, pd );
02521     pd.element_by_index( 0 ).get_vertex_indices( conn );
02522 
02523     qm->get_evaluations( pd, handles, false, err );
02524     ASSERT_NO_ERROR( err );
02525     // For sample-based metrics, it is difficult to set up such that
02526     // both surface and volume elements have only one sample point.
02527     // Make things easier by skipping element types with no sample points.
02528     if( handles.empty() ) return;
02529     CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics", (size_t)1, handles.size() );
02530     size_t h = handles[0];
02531 
02532     rval = qm->evaluate_with_gradient( pd, h, val, indices, grad, err );
02533     ASSERT_NO_ERROR( err );
02534     CPPUNIT_ASSERT( rval );
02535     CPPUNIT_ASSERT_EQUAL( (size_t)n, indices.size() );
02536     CPPUNIT_ASSERT_EQUAL( (size_t)n, grad.size() );
02537 
02538     for( int i = 0; i < n; ++i )
02539     {
02540         get_nonideal_element( type, pd, i );
02541 
02542         qm->get_evaluations( pd, handles, false, err );
02543         ASSERT_NO_ERROR( err );
02544         CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics", (size_t)1, handles.size() );
02545         h = handles[0];
02546 
02547         rval = qm->evaluate_with_gradient( pd, h, val_fixed, indices_fixed, grad_fixed, err );
02548         ASSERT_NO_ERROR( err );
02549         CPPUNIT_ASSERT( rval );
02550         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices_fixed.size() );
02551         CPPUNIT_ASSERT_EQUAL( (size_t)1, grad_fixed.size() );
02552 
02553         CPPUNIT_ASSERT_DOUBLES_EQUAL( val, val_fixed, 1e-5 );
02554         CPPUNIT_ASSERT_VECTORS_EQUAL( grad[conn[i]], grad_fixed.front(), 1e-5 );
02555     }
02556 }
02557 
02558 void QualityMetricTester::test_hessian_with_fixed_vertex( EntityTopology type,
02559                                                           QualityMetric* qm,
02560                                                           const Settings* settings )
02561 {
02562     MsqPrintError err( std::cout );
02563     PatchData pd;
02564     if( settings )
02565         pd.attach_settings( settings );
02566     else if( mSettings )
02567         pd.attach_settings( mSettings );
02568     std::vector< size_t > handles, indices, indices_fixed, conn;
02569     std::vector< Vector3D > grad, grad_fixed;
02570     std::vector< Matrix3D > hess, hess_fixed;
02571     double val, val_fixed;
02572     bool rval;
02573     const int n = TopologyInfo::corners( type );
02574 
02575     get_nonideal_element( type, pd );
02576     pd.element_by_index( 0 ).get_vertex_indices( conn );
02577 
02578     qm->get_evaluations( pd, handles, false, err );
02579     ASSERT_NO_ERROR( err );
02580     // For sample-based metrics, it is difficult to set up such that
02581     // both surface and volume elements have only one sample point.
02582     // Make things easier by skipping element types with no sample points.
02583     if( handles.empty() ) return;
02584     CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics", (size_t)1, handles.size() );
02585     size_t h = handles[0];
02586 
02587     rval = qm->evaluate_with_Hessian( pd, h, val, indices, grad, hess, err );
02588     ASSERT_NO_ERROR( err );
02589     CPPUNIT_ASSERT( rval );
02590     CPPUNIT_ASSERT_EQUAL( (size_t)n, indices.size() );
02591     CPPUNIT_ASSERT_EQUAL( (size_t)n, grad.size() );
02592     CPPUNIT_ASSERT_EQUAL( (size_t)( n * ( n + 1 ) / 2 ), hess.size() );
02593 
02594     for( int i = 0; i < n; ++i )
02595     {
02596         get_nonideal_element( type, pd, i );
02597 
02598         qm->get_evaluations( pd, handles, false, err );
02599         ASSERT_NO_ERROR( err );
02600         CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics", (size_t)1, handles.size() );
02601         h = handles[0];
02602 
02603         rval = qm->evaluate_with_Hessian( pd, h, val_fixed, indices_fixed, grad_fixed, hess_fixed, err );
02604         ASSERT_NO_ERROR( err );
02605         CPPUNIT_ASSERT( rval );
02606         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices_fixed.size() );
02607         CPPUNIT_ASSERT_EQUAL( (size_t)1, grad_fixed.size() );
02608         CPPUNIT_ASSERT_EQUAL( (size_t)1, hess_fixed.size() );
02609 
02610         const size_t j = conn[i];
02611         CPPUNIT_ASSERT_DOUBLES_EQUAL( val, val_fixed, 1e-5 );
02612         CPPUNIT_ASSERT_VECTORS_EQUAL( grad[j], grad_fixed.front(), 1e-5 );
02613         CPPUNIT_ASSERT_MATRICES_EQUAL( hess[n * j - j * ( j - 1 ) / 2], hess_fixed.front(), 1e-4 );
02614     }
02615 }
02616 
02617 void QualityMetricTester::test_diagonal_with_fixed_vertex( EntityTopology type,
02618                                                            QualityMetric* qm,
02619                                                            const Settings* settings )
02620 {
02621     MsqPrintError err( std::cout );
02622     PatchData pd;
02623     if( settings )
02624         pd.attach_settings( settings );
02625     else if( mSettings )
02626         pd.attach_settings( mSettings );
02627     std::vector< size_t > handles, indices, indices_fixed, conn;
02628     std::vector< Vector3D > grad, grad_fixed;
02629     std::vector< SymMatrix3D > hess, hess_fixed;
02630     double val, val_fixed;
02631     bool rval;
02632     const int n = TopologyInfo::corners( type );
02633 
02634     get_nonideal_element( type, pd );
02635     pd.element_by_index( 0 ).get_vertex_indices( conn );
02636 
02637     qm->get_evaluations( pd, handles, false, err );
02638     ASSERT_NO_ERROR( err );
02639     // For sample-based metrics, it is difficult to set up such that
02640     // both surface and volume elements have only one sample point.
02641     // Make things easier by skipping element types with no sample points.
02642     if( handles.empty() ) return;
02643     CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics", (size_t)1, handles.size() );
02644     size_t h = handles[0];
02645 
02646     rval = qm->evaluate_with_Hessian_diagonal( pd, h, val, indices, grad, hess, err );
02647     ASSERT_NO_ERROR( err );
02648     CPPUNIT_ASSERT( rval );
02649     CPPUNIT_ASSERT_EQUAL( (size_t)n, indices.size() );
02650     CPPUNIT_ASSERT_EQUAL( (size_t)n, grad.size() );
02651     CPPUNIT_ASSERT_EQUAL( (size_t)n, hess.size() );
02652 
02653     for( int i = 0; i < n; ++i )
02654     {
02655         get_nonideal_element( type, pd, i );
02656 
02657         qm->get_evaluations( pd, handles, false, err );
02658         ASSERT_NO_ERROR( err );
02659         CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics", (size_t)1, handles.size() );
02660         h = handles[0];
02661 
02662         rval = qm->evaluate_with_Hessian_diagonal( pd, h, val_fixed, indices_fixed, grad_fixed, hess_fixed, err );
02663         ASSERT_NO_ERROR( err );
02664         CPPUNIT_ASSERT( rval );
02665         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices_fixed.size() );
02666         CPPUNIT_ASSERT_EQUAL( (size_t)1, grad_fixed.size() );
02667         CPPUNIT_ASSERT_EQUAL( (size_t)1, hess_fixed.size() );
02668 
02669         const size_t j = conn[i];
02670         CPPUNIT_ASSERT_DOUBLES_EQUAL( val, val_fixed, 1e-5 );
02671         CPPUNIT_ASSERT_VECTORS_EQUAL( grad[j], grad_fixed.front(), 1e-5 );
02672         CPPUNIT_ASSERT_MATRICES_EQUAL( hess[j], hess_fixed.front(), 1e-4 );
02673     }
02674 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines