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