MOAB: Mesh Oriented datABase  (version 5.2.1)
TMPQualityMetricTest.hpp
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     retian 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     (2010) kraftche@cae.wisc.edu
00025 
00026   ***************************************************************** */
00027 
00028 /*! \file TMPQualityMetricTest.cpp
00029 
00030 Unit testing for the TMPQualityMetric class
00031 \author Jasno Kraftcheck
00032 */
00033 
00034 #include "IdealShapeTarget.hpp"
00035 #include "MsqMatrix.hpp"
00036 #include "QualityMetricTester.hpp"
00037 #include "Settings.hpp"
00038 #include "UnitUtil.hpp"
00039 #include "PlanarDomain.hpp"
00040 #include "PatchData.hpp"
00041 #include "WeightCalculator.hpp"
00042 #include "ElementPMeanP.hpp"
00043 #include "ElemSampleQM.hpp"
00044 
00045 #include <iostream>
00046 
00047 using namespace MBMesquite;
00048 using std::cerr;
00049 using std::cout;
00050 using std::endl;
00051 
00052 /** Target metric (templatized by dimension) for use in misc. tests.
00053  *  'evaluate' method records input values and returns a constant.
00054  */
00055 template < typename B >
00056 class FauxMetric : public B
00057 {
00058   public:
00059     int count;
00060     double value;
00061     bool rval;
00062     MsqMatrix< 2, 2 > last_A_2D;
00063     MsqMatrix< 3, 3 > last_A_3D;
00064 
00065     FauxMetric( double v ) : count( 0 ), value( v ), rval( true ) {}
00066 
00067     std::string get_name() const
00068     {
00069         return "Faux";
00070     }
00071 
00072     bool evaluate( const MsqMatrix< 2, 2 >& A, const MsqMatrix< 2, 2 >& W, double& result, MsqError& )
00073     {
00074         last_A_2D = A;
00075         result    = value;
00076         ++count;
00077         return rval;
00078     }
00079 
00080     bool evaluate( const MsqMatrix< 3, 3 >& A, const MsqMatrix< 3, 3 >& W, double& result, MsqError& )
00081     {
00082         last_A_3D = A;
00083         result    = value;
00084         ++count;
00085         return rval;
00086     }
00087 
00088     bool evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& )
00089     {
00090         last_A_2D = T;
00091         result    = value;
00092         ++count;
00093         return rval;
00094     }
00095 
00096     bool evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& )
00097     {
00098         last_A_3D = T;
00099         result    = value;
00100         ++count;
00101         return rval;
00102     }
00103 };
00104 
00105 /** Weight calculator used for testing.  Returns constant weight. */
00106 class ScaleWeight : public WeightCalculator
00107 {
00108   public:
00109     ScaleWeight( double s ) : value( s ) {}
00110     double get_weight( PatchData&, size_t, Sample, MsqError& )
00111     {
00112         return value;
00113     }
00114     double value;
00115 };
00116 
00117 /** wrapper class to force numeric approximation of derivatives */
00118 template < class Base >
00119 class NumericalMetric : public Base
00120 {
00121   public:
00122     NumericalMetric( Base* real_metric ) : mMetric( real_metric ) {}
00123 
00124     ~NumericalMetric() {}
00125 
00126     std::string get_name() const
00127     {
00128         return "Numerical " + mMetric->get_name();
00129     }
00130 
00131     bool evaluate( const MsqMatrix< 2, 2 >& A, const MsqMatrix< 2, 2 >& W, double& result, MsqError& err )
00132     {
00133         return mMetric->evaluate( A, W, result, err );
00134     }
00135 
00136     bool evaluate( const MsqMatrix< 3, 3 >& A, const MsqMatrix< 3, 3 >& W, double& result, MsqError& err )
00137     {
00138         return mMetric->evaluate( A, W, result, err );
00139     }
00140 
00141     bool evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& err )
00142     {
00143         return mMetric->evaluate( T, result, err );
00144     }
00145 
00146     bool evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& err )
00147     {
00148         return mMetric->evaluate( T, result, err );
00149     }
00150 
00151   private:
00152     Base* mMetric;
00153 };
00154 
00155 /** Simple target metric for testing first partial derivatives.
00156  *  \f$\mu(A,W) = |A|^2\f$
00157  *  \f$\frac{\partial\mu}{\partial \A} = 2 A \f$
00158  */
00159 template < class Base >
00160 class TestGradTargetMetric : public Base
00161 {
00162   public:
00163     std::string get_name() const
00164     {
00165         return "TestGrad";
00166     }
00167 
00168     bool evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& err )
00169     {
00170         result = sqr_Frobenius( T );
00171         return true;
00172     }
00173 
00174     bool evaluate( const MsqMatrix< 2, 2 >& A, const MsqMatrix< 2, 2 >&, double& result, MsqError& err )
00175     {
00176         return evaluate( A, result, err );
00177     }
00178 
00179     bool evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& d, MsqError& err )
00180     {
00181         result = sqr_Frobenius( T );
00182         d      = 2 * T;
00183         return true;
00184     }
00185 
00186     bool evaluate_with_grad( const MsqMatrix< 2, 2 >& A, const MsqMatrix< 2, 2 >&, double& result, MsqMatrix< 2, 2 >& d,
00187                              MsqError& err )
00188     {
00189         return evaluate_with_grad( A, result, d, err );
00190     }
00191 
00192     bool evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& err )
00193     {
00194         result = sqr_Frobenius( T );
00195         return true;
00196     }
00197 
00198     bool evaluate( const MsqMatrix< 3, 3 >& A, const MsqMatrix< 3, 3 >&, double& result, MsqError& err )
00199     {
00200         return evaluate( A, result, err );
00201     }
00202 
00203     bool evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& d, MsqError& err )
00204     {
00205         result = sqr_Frobenius( T );
00206         d      = 2 * T;
00207         return true;
00208     }
00209 
00210     bool evaluate_with_grad( const MsqMatrix< 3, 3 >& A, const MsqMatrix< 3, 3 >&, double& result, MsqMatrix< 3, 3 >& d,
00211                              MsqError& err )
00212     {
00213         return evaluate_with_grad( A, result, d, err );
00214     }
00215 };
00216 
00217 /* class to force evaluation of mapping function only at element center
00218  * so that we can re-use tests from QualityMetricTester that work only
00219  * for element-based metrics (TMP metric is "element based" if only one
00220  * sample in each element.)
00221  */
00222 class CenterMF2D : public MappingFunction2D
00223 {
00224   public:
00225     CenterMF2D( const MappingFunction2D* real_mf ) : myFunc( real_mf ){};
00226     EntityTopology element_topology() const
00227     {
00228         return myFunc->element_topology();
00229     }
00230     int num_nodes() const
00231     {
00232         return myFunc->num_nodes();
00233     }
00234     NodeSet sample_points( NodeSet ) const
00235     {
00236         NodeSet s;
00237         s.set_mid_face_node( 0 );
00238         return s;
00239     }
00240     void coefficients( Sample l, NodeSet s, double* c, size_t* i, size_t& n, MsqError& e ) const
00241     {
00242         myFunc->coefficients( l, s, c, i, n, e );
00243     }
00244     void derivatives( Sample l, NodeSet s, size_t* i, MsqVector< 2 >* c, size_t& n, MsqError& e ) const
00245     {
00246         myFunc->derivatives( l, s, i, c, n, e );
00247     }
00248 
00249   private:
00250     const MappingFunction2D* myFunc;
00251 };
00252 class CenterMF3D : public MappingFunction3D
00253 {
00254   public:
00255     CenterMF3D( const MappingFunction3D* real_mf ) : myFunc( real_mf ){};
00256     EntityTopology element_topology() const
00257     {
00258         return myFunc->element_topology();
00259     }
00260     int num_nodes() const
00261     {
00262         return myFunc->num_nodes();
00263     }
00264     NodeSet sample_points( NodeSet ) const
00265     {
00266         NodeSet s;
00267         s.set_mid_region_node();
00268         return s;
00269     }
00270     void coefficients( Sample l, NodeSet s, double* c, size_t* i, size_t& n, MsqError& e ) const
00271     {
00272         myFunc->coefficients( l, s, c, i, n, e );
00273     }
00274     void derivatives( Sample l, NodeSet s, size_t* i, MsqVector< 3 >* c, size_t& n, MsqError& e ) const
00275     {
00276         myFunc->derivatives( l, s, i, c, n, e );
00277     }
00278 
00279   private:
00280     const MappingFunction3D* myFunc;
00281 };
00282 
00283 // Define a target calculator that returns targets for
00284 // ideally shaped elements, but also includes orientation
00285 // information (aligning surface elements to the xy plane
00286 // with the first column of the jacobian in the x direction).
00287 class IdealShapeXY : public IdealShapeTarget
00288 {
00289   public:
00290     bool have_surface_orient() const
00291     {
00292         return true;
00293     }
00294     bool get_surface_target( PatchData& pd, size_t element, Sample sample, MsqMatrix< 3, 2 >& W_out, MsqError& err )
00295     {
00296         MsqMatrix< 2, 2 > W;
00297         bool rval = get_2D_target( pd, element, sample, W, err );
00298         W_out.set_row( 0, W.row( 0 ) );
00299         W_out.set_row( 1, W.row( 1 ) );
00300         W_out.set_row( 2, MsqMatrix< 1, 2 >( 0.0 ) );
00301         return rval;
00302     }
00303 };
00304 
00305 #define REGISTER_TMP_TESTS                                      \
00306                                                                 \
00307     CPPUNIT_TEST( test_negate_flag );                           \
00308     CPPUNIT_TEST( test_supported_types );                       \
00309     CPPUNIT_TEST( test_get_evaluations );                       \
00310     CPPUNIT_TEST( test_get_element_evaluations );               \
00311                                                                 \
00312     CPPUNIT_TEST( test_evaluate_2D );                           \
00313     CPPUNIT_TEST( test_evaluate_surface );                      \
00314     CPPUNIT_TEST( test_evaluate_3D );                           \
00315     CPPUNIT_TEST( test_evaluate_2D_weight );                    \
00316     CPPUNIT_TEST( test_evaluate_surface_weight );               \
00317     CPPUNIT_TEST( test_evaluate_3D_weight );                    \
00318     CPPUNIT_TEST( test_2d_eval_ortho_quad );                    \
00319     CPPUNIT_TEST( test_surf_eval_ortho_quad );                  \
00320     CPPUNIT_TEST( test_3d_eval_ortho_hex );                     \
00321                                                                 \
00322     CPPUNIT_TEST( test_sample_indices );                        \
00323     CPPUNIT_TEST( test_evaluate_with_indices );                 \
00324     CPPUNIT_TEST( test_evaluate_fixed_indices );                \
00325                                                                 \
00326     CPPUNIT_TEST( test_gradient_2D );                           \
00327     CPPUNIT_TEST( test_gradient_surface );                      \
00328     CPPUNIT_TEST( test_gradient_3D );                           \
00329     CPPUNIT_TEST( compare_indices_and_gradient );               \
00330     CPPUNIT_TEST( test_ideal_element_gradient );                \
00331     CPPUNIT_TEST( compare_analytical_and_numerical_gradient );  \
00332     CPPUNIT_TEST( test_weighted_gradients );                    \
00333     CPPUNIT_TEST( test_gradient_with_fixed_vertices );          \
00334                                                                 \
00335     CPPUNIT_TEST( compare_indices_and_hessian );                \
00336     CPPUNIT_TEST( compare_gradient_and_hessian );               \
00337     CPPUNIT_TEST( compare_analytical_and_numerical_hessians );  \
00338     CPPUNIT_TEST( test_symmetric_hessian_diagonal );            \
00339     CPPUNIT_TEST( test_weighted_hessians );                     \
00340     CPPUNIT_TEST( test_hessian_with_fixed_vertices );           \
00341                                                                 \
00342     CPPUNIT_TEST( compare_indices_and_diagonal );               \
00343     CPPUNIT_TEST( compare_gradient_and_diagonal );              \
00344     CPPUNIT_TEST( compare_analytical_and_numerical_diagonals ); \
00345     CPPUNIT_TEST( test_weighted_diagonals );                    \
00346     CPPUNIT_TEST( test_diagonal_with_fixed_vertices );
00347 
00348 static double col_dot_prod( MsqMatrix< 2, 2 >& m )
00349 {
00350     return m( 0, 0 ) * m( 0, 1 ) + m( 1, 0 ) * m( 1, 1 );
00351 }
00352 
00353 template < class QMType >
00354 class TMPTypes
00355 {
00356 };
00357 
00358 template < class QMType >
00359 class TMPQualityMetricTest : public CppUnit::TestFixture
00360 {
00361   protected:
00362     QualityMetricTester tester;
00363 
00364     Settings settings;
00365     IdealShapeTarget ideal;
00366     IdealShapeXY surf_target;
00367     ScaleWeight e_weight;
00368 
00369     FauxMetric< typename TMPTypes< QMType >::MetricType > faux_pi, faux_zero, faux_two;
00370     typename TMPTypes< QMType >::TestType test_metric;
00371     NumericalMetric< typename QMType::MetricType > num_metric;
00372     QMType test_qm, test_qm_surf, zero_qm, weight_qm, center_qm;
00373     Settings centerOnly;
00374     CenterMF2D triCenter, quadCenter;
00375     CenterMF3D tetCenter, pyrCenter, priCenter, hexCenter;
00376 
00377   public:
00378     TMPQualityMetricTest()
00379         : tester( QualityMetricTester::ALL_FE_EXCEPT_SEPTAHEDRON, &settings ), e_weight( 2.7182818284590451 ),
00380           faux_pi( 3.14159 ), faux_zero( 0.0 ), faux_two( 2.0 ), num_metric( &test_metric ),
00381           test_qm( &ideal, &num_metric ), test_qm_surf( &surf_target, &num_metric ), zero_qm( &ideal, &faux_zero ),
00382           weight_qm( &ideal, &e_weight, &test_metric ), center_qm( &ideal, &test_metric ),
00383           triCenter( centerOnly.get_mapping_function_2D( TRIANGLE ) ),
00384           quadCenter( centerOnly.get_mapping_function_2D( QUADRILATERAL ) ),
00385           tetCenter( centerOnly.get_mapping_function_3D( TETRAHEDRON ) ),
00386           pyrCenter( centerOnly.get_mapping_function_3D( PYRAMID ) ),
00387           priCenter( centerOnly.get_mapping_function_3D( PRISM ) ),
00388           hexCenter( centerOnly.get_mapping_function_3D( HEXAHEDRON ) )
00389     {
00390         centerOnly.set_mapping_function( &triCenter );
00391         centerOnly.set_mapping_function( &quadCenter );
00392         centerOnly.set_mapping_function( &tetCenter );
00393         centerOnly.set_mapping_function( &pyrCenter );
00394         centerOnly.set_mapping_function( &priCenter );
00395         centerOnly.set_mapping_function( &hexCenter );
00396         tester.ideal_pyramid_base_equals_height( true );
00397     }
00398 
00399     void test_negate_flag()
00400     {
00401         CPPUNIT_ASSERT_EQUAL( 1, zero_qm.get_negate_flag() );
00402     }
00403     void test_supported_types()
00404     {
00405         tester.test_supported_element_types( &zero_qm );
00406     }
00407     void test_get_evaluations()
00408     {
00409         QMType edge_metric( &ideal, &faux_zero );
00410         tester.test_get_sample_evaluations( &zero_qm );
00411         tester.test_get_sample_evaluations( &edge_metric );
00412     }
00413     void test_get_element_evaluations()
00414     {
00415         QMType edge_metric( &ideal, &faux_zero );
00416         tester.test_get_in_element_evaluations( &zero_qm );
00417         tester.test_get_in_element_evaluations( &edge_metric );
00418     }
00419 
00420     void test_evaluate_2D();
00421     void test_evaluate_surface();
00422     void test_evaluate_3D();
00423 
00424     void test_evaluate_2D_weight()
00425     {
00426         MsqPrintError err( cout );
00427         PatchData pd;
00428         bool rval;
00429         double value;
00430 
00431         QMType m( &ideal, &e_weight, &faux_pi );
00432         tester.get_ideal_element( TRIANGLE, true, pd );
00433         rval = m.evaluate( pd, 0, value, err );
00434         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00435         CPPUNIT_ASSERT( rval );
00436         CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value * e_weight.value, value, DBL_EPSILON );
00437     }
00438 
00439     void test_evaluate_surface_weight()
00440     {
00441         MsqPrintError err( cout );
00442         PatchData pd;
00443         bool rval;
00444         double value;
00445 
00446         QMType m( &surf_target, &e_weight, &faux_pi );
00447 
00448         tester.get_ideal_element( TRIANGLE, true, pd );
00449         rval = m.evaluate( pd, 0, value, err );
00450         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00451         CPPUNIT_ASSERT( rval );
00452         CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value * e_weight.value, value, DBL_EPSILON );
00453     }
00454 
00455     void test_evaluate_3D_weight()
00456     {
00457         MsqPrintError err( cout );
00458         PatchData pd;
00459         bool rval;
00460         double value;
00461 
00462         QMType m( &ideal, &e_weight, &faux_two );
00463 
00464         tester.get_ideal_element( PRISM, true, pd );
00465         rval = m.evaluate( pd, 0, value, err );
00466         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00467         CPPUNIT_ASSERT( rval );
00468         CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_two.value * e_weight.value, value, DBL_EPSILON );
00469     }
00470 
00471     void test_2d_eval_ortho_quad()
00472     {
00473         MsqPrintError err( cout );
00474         PatchData pd;
00475         bool rval;
00476         double value;
00477 
00478         QMType m( &ideal, &faux_zero );
00479         faux_zero.count = 0;
00480 
00481         tester.get_ideal_element( QUADRILATERAL, true, pd );
00482         rval = m.evaluate( pd, 0, value, err );
00483         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00484         CPPUNIT_ASSERT( rval );
00485         CPPUNIT_ASSERT_EQUAL( 1, faux_zero.count );
00486         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, col_dot_prod( faux_zero.last_A_2D ), DBL_EPSILON );
00487     }
00488 
00489     void test_surf_eval_ortho_quad()
00490     {
00491         MsqPrintError err( cout );
00492         PatchData pd;
00493         bool rval;
00494         double value;
00495 
00496         QMType m( &surf_target, &faux_zero );
00497         faux_zero.count = 0;
00498 
00499         tester.get_ideal_element( QUADRILATERAL, true, pd );
00500         rval = m.evaluate( pd, 0, value, err );
00501         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00502         CPPUNIT_ASSERT( rval );
00503         CPPUNIT_ASSERT_EQUAL( 1, faux_zero.count );
00504         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, col_dot_prod( faux_zero.last_A_2D ), DBL_EPSILON );
00505     }
00506 
00507     void test_3d_eval_ortho_hex()
00508     {
00509         MsqPrintError err( cout );
00510         PatchData pd;
00511         bool rval;
00512         double value;
00513 
00514         QMType m( &ideal, &faux_zero );
00515         faux_zero.count = 0;
00516 
00517         tester.get_ideal_element( HEXAHEDRON, true, pd );
00518         rval = m.evaluate( pd, 0, value, err );
00519         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00520         CPPUNIT_ASSERT( rval );
00521         CPPUNIT_ASSERT_EQUAL( 1, faux_zero.count );
00522 
00523         // test that columns are orthogonal for ideal hex element
00524         MsqMatrix< 3, 3 > A = faux_zero.last_A_3D;
00525         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 0 ) % A.column( 1 ), 1e-6 );
00526         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 0 ) % A.column( 2 ), 1e-6 );
00527         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 1 ) % A.column( 2 ), 1e-6 );
00528     }
00529 
00530     void test_gradient_common( TargetCalculator* tc );
00531     void test_gradient_2D()
00532     {
00533         test_gradient_common( &ideal );
00534     }
00535     void test_gradient_surface()
00536     {
00537         test_gradient_common( &surf_target );
00538     }
00539     void test_gradient_3D();
00540 
00541     void test_sample_indices()
00542     {
00543         tester.test_get_sample_indices( &zero_qm );
00544     }
00545     void test_evaluate_with_indices()
00546     {
00547         tester.compare_eval_and_eval_with_indices( &zero_qm );
00548     }
00549     void test_evaluate_fixed_indices()
00550     {
00551         tester.test_get_indices_fixed( &zero_qm );
00552     }
00553 
00554     void compare_indices_and_gradient()
00555     {
00556         tester.compare_eval_with_indices_and_eval_with_gradient( &test_qm );
00557         tester.compare_eval_with_indices_and_eval_with_gradient( &test_qm_surf );
00558     }
00559     void test_ideal_element_gradient()
00560     {
00561         tester.test_ideal_element_zero_gradient( &test_qm, false );
00562         tester.test_ideal_element_zero_gradient( &test_qm_surf, false );
00563     }
00564     void compare_analytical_and_numerical_gradient()
00565     {
00566         compare_analytical_and_numerical_gradients( &test_qm );
00567         compare_analytical_and_numerical_gradients( &test_qm_surf );
00568     }
00569     void test_weighted_gradients()
00570     {
00571         compare_analytical_and_numerical_gradients( &weight_qm );
00572     }
00573     void test_gradient_with_fixed_vertices()
00574     {
00575         tester.test_gradient_with_fixed_vertex( &center_qm, &centerOnly );
00576     }
00577 
00578     void compare_indices_and_hessian()
00579     {
00580         tester.compare_eval_with_indices_and_eval_with_hessian( &test_qm );
00581         tester.compare_eval_with_indices_and_eval_with_hessian( &test_qm_surf );
00582     }
00583     void compare_gradient_and_hessian()
00584     {
00585         tester.compare_eval_with_grad_and_eval_with_hessian( &test_qm );
00586         tester.compare_eval_with_grad_and_eval_with_hessian( &test_qm_surf );
00587     }
00588     void compare_analytical_and_numerical_hessians()
00589     {
00590         compare_analytical_and_numerical_hessians( &test_qm );
00591         compare_analytical_and_numerical_hessians( &test_qm_surf );
00592     }
00593     void test_symmetric_hessian_diagonal()
00594     {
00595         tester.test_symmetric_Hessian_diagonal_blocks( &test_qm );
00596         tester.test_symmetric_Hessian_diagonal_blocks( &test_qm_surf );
00597     }
00598     void test_weighted_hessians()
00599     {
00600         compare_analytical_and_numerical_hessians( &weight_qm );
00601     }
00602     void test_hessian_with_fixed_vertices()
00603     {
00604         tester.test_hessian_with_fixed_vertex( &center_qm, &centerOnly );
00605     }
00606 
00607     void compare_indices_and_diagonal()
00608     {
00609         tester.compare_eval_with_indices_and_eval_with_diagonal( &test_qm );
00610         tester.compare_eval_with_indices_and_eval_with_diagonal( &test_qm_surf );
00611     }
00612     void compare_gradient_and_diagonal()
00613     {
00614         tester.compare_eval_with_grad_and_eval_with_diagonal( &test_qm );
00615         tester.compare_eval_with_grad_and_eval_with_diagonal( &test_qm_surf );
00616     }
00617     void compare_analytical_and_numerical_diagonals()
00618     {
00619         compare_analytical_and_numerical_diagonals( &test_qm );
00620         compare_analytical_and_numerical_diagonals( &test_qm_surf );
00621     }
00622     void test_weighted_diagonals()
00623     {
00624         compare_analytical_and_numerical_diagonals( &weight_qm );
00625     }
00626     void test_diagonal_with_fixed_vertices()
00627     {
00628         tester.test_diagonal_with_fixed_vertex( &center_qm, &centerOnly );
00629     }
00630 
00631     // Delcare specialized versions of the functions from
00632     // QualityMetricTester because we surface elements must
00633     // be handled differently.  For a surface element in the XY plane,
00634     // the finite difference approximations of the derivatives will
00635     // have non-zero values for derivatives wrt Z coordinates while the
00636     // analytical derivative calculations will return all derivatives
00637     // wrt Z coordiantes as zero.
00638 
00639     void get_nonideal_element( EntityTopology type, PatchData& pd )
00640     {
00641         tester.get_nonideal_element( type, pd, true );
00642         // Callers assume surface elements are in XY plane.
00643         // Verify this assumption.
00644         if( TopologyInfo::dimension( type ) == 2 )
00645         {
00646             for( size_t i = 0; i < pd.num_nodes(); ++i )
00647             {
00648                 CPPUNIT_ASSERT_DOUBLES_EQUAL( pd.vertex_by_index( i )[2], 0.0, 1e-6 );
00649             }
00650         }
00651     }
00652 
00653     void compare_analytical_and_numerical_gradients( QualityMetric* qm )
00654     {
00655         PatchData pd;
00656         const EntityTopology types[] = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON };
00657         const int num_types          = sizeof( types ) / sizeof( types[0] );
00658         for( int i = 0; i < num_types; ++i )
00659         {
00660             get_nonideal_element( types[i], pd );
00661             compare_analytical_and_numerical_gradients( qm, pd, TopologyInfo::dimension( types[i] ) );
00662         }
00663     }
00664 
00665     void compare_analytical_and_numerical_hessians( QualityMetric* qm );
00666     void compare_analytical_and_numerical_diagonals( QualityMetric* qm );
00667     void compare_analytical_and_numerical_gradients( QualityMetric* qm, PatchData&, int dim );
00668 };
00669 
00670 // CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(TMPQualityMetricTest, "TMPQualityMetricTest");
00671 // CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(TMPQualityMetricTest, "Unit");
00672 
00673 template < class QMType >
00674 inline void TMPQualityMetricTest< QMType >::test_evaluate_2D()
00675 {
00676     MsqPrintError err( cout );
00677     PatchData pd;
00678     bool rval;
00679     double value;
00680 
00681     QMType m( &ideal, &faux_pi );
00682 
00683     // test with aligned elements
00684     faux_pi.count = 0;
00685     tester.get_ideal_element( QUADRILATERAL, true, pd );
00686     rval = m.evaluate( pd, 0, value, err );
00687     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00688     CPPUNIT_ASSERT( rval );
00689     CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
00690     CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );
00691 
00692     // test that columns are orthogonal for ideal quad element
00693     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, col_dot_prod( faux_pi.last_A_2D ), 1e-6 );
00694 
00695     // test with an element rotated about X-axis
00696     faux_pi.count = 0;
00697     tester.get_ideal_element( QUADRILATERAL, true, pd );
00698     // rotate by 90 degrees about X axis
00699     for( size_t i = 0; i < pd.num_nodes(); ++i )
00700     {
00701         Vector3D orig = pd.vertex_by_index( i );
00702         Vector3D newp( orig[0], -orig[2], orig[1] );
00703         pd.set_vertex_coordinates( newp, i, err );
00704     }
00705     rval = m.evaluate( pd, 0, value, err );
00706     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00707     CPPUNIT_ASSERT( rval );
00708     CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
00709     CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );
00710 
00711     // test with an element rotated about Y-axis
00712     faux_pi.count = 0;
00713     tester.get_ideal_element( TRIANGLE, true, pd );
00714     // rotate by -90 degrees about Y axis
00715     for( size_t i = 0; i < pd.num_nodes(); ++i )
00716     {
00717         Vector3D orig = pd.vertex_by_index( i );
00718         Vector3D newp( orig[2], orig[1], -orig[0] );
00719         pd.set_vertex_coordinates( newp, i, err );
00720     }
00721     rval = m.evaluate( pd, 0, value, err );
00722     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00723     CPPUNIT_ASSERT( rval );
00724     CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
00725     CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );
00726 }
00727 
00728 template < class QMType >
00729 inline void TMPQualityMetricTest< QMType >::test_evaluate_surface()
00730 {
00731     MsqPrintError err( cout );
00732     PatchData pd;
00733     bool rval;
00734     double value;
00735 
00736     QMType m( &surf_target, &faux_pi );
00737 
00738     // test with aligned elements
00739     faux_pi.count = 0;
00740     tester.get_ideal_element( QUADRILATERAL, true, pd );
00741     rval = m.evaluate( pd, 0, value, err );
00742     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00743     CPPUNIT_ASSERT( rval );
00744     CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
00745     CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );
00746 
00747     // test that columns are orthogonal for ideal quad element
00748     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, col_dot_prod( faux_pi.last_A_2D ), 1e-6 );
00749 
00750     // test with an element rotated about X-axis
00751     faux_pi.count = 0;
00752     tester.get_ideal_element( QUADRILATERAL, true, pd );
00753     // rotate by 90 degrees about X axis
00754     for( size_t i = 0; i < pd.num_nodes(); ++i )
00755     {
00756         Vector3D orig = pd.vertex_by_index( i );
00757         Vector3D newp( orig[0], -orig[2], orig[1] );
00758         pd.set_vertex_coordinates( newp, i, err );
00759     }
00760     rval = m.evaluate( pd, 0, value, err );
00761     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00762     CPPUNIT_ASSERT( rval );
00763     CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
00764     CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );
00765 
00766     // test with an element rotated about Y-axis
00767     faux_pi.count = 0;
00768     tester.get_ideal_element( TRIANGLE, true, pd );
00769     // rotate by -90 degrees about Y axis
00770     for( size_t i = 0; i < pd.num_nodes(); ++i )
00771     {
00772         Vector3D orig = pd.vertex_by_index( i );
00773         Vector3D newp( orig[2], orig[1], -orig[0] );
00774         pd.set_vertex_coordinates( newp, i, err );
00775     }
00776     rval = m.evaluate( pd, 0, value, err );
00777     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00778     CPPUNIT_ASSERT( rval );
00779     CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_pi.value, value, DBL_EPSILON );
00780     CPPUNIT_ASSERT_EQUAL( 1, faux_pi.count );
00781 }
00782 
00783 template < class QMType >
00784 inline void TMPQualityMetricTest< QMType >::test_evaluate_3D()
00785 {
00786     MsqPrintError err( cout );
00787     PatchData pd;
00788     bool rval;
00789     double value;
00790 
00791     QMType m( &ideal, &faux_two );
00792 
00793     // test with aligned elements
00794     faux_two.count = 0;
00795     tester.get_ideal_element( HEXAHEDRON, true, pd );
00796     rval = m.evaluate( pd, 0, value, err );
00797     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00798     CPPUNIT_ASSERT( rval );
00799     CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_two.value, value, DBL_EPSILON );
00800     CPPUNIT_ASSERT_EQUAL( 1, faux_two.count );
00801 
00802     // test that columns are orthogonal for ideal hex element
00803     MsqMatrix< 3, 3 > A = faux_two.last_A_3D;
00804     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 0 ) % A.column( 1 ), 1e-6 );
00805     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 0 ) % A.column( 2 ), 1e-6 );
00806     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, A.column( 1 ) % A.column( 2 ), 1e-6 );
00807 
00808     // test with rotated element
00809     faux_two.count = 0;
00810     tester.get_ideal_element( TETRAHEDRON, true, pd );
00811     // rotate by 90-degrees about X axis
00812     for( size_t i = 0; i < pd.num_nodes(); ++i )
00813     {
00814         Vector3D orig = pd.vertex_by_index( i );
00815         Vector3D newp( orig[0], -orig[2], orig[1] );
00816         pd.set_vertex_coordinates( newp, i, err );
00817     }
00818     rval = m.evaluate( pd, 0, value, err );
00819     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00820     CPPUNIT_ASSERT( rval );
00821     CPPUNIT_ASSERT_DOUBLES_EQUAL( faux_two.value, value, DBL_EPSILON );
00822     CPPUNIT_ASSERT_EQUAL( 1, faux_two.count );
00823 }
00824 
00825 template < class QMType >
00826 inline void TMPQualityMetricTest< QMType >::test_gradient_common( TargetCalculator* tc )
00827 {
00828     MsqPrintError err( std::cout );
00829 
00830     // check for expected value at center of flattened hex
00831 
00832     // construct flattened hex
00833     const double y          = 0.5;
00834     const double vertices[] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, y, 0.0, 0.0, y, 0.0 };
00835     size_t conn[8]          = { 0, 1, 2, 3 };
00836     PatchData pd;
00837     pd.fill( 4, vertices, 1, QUADRILATERAL, conn, 0, err );
00838     ASSERT_NO_ERROR( err );
00839 
00840     // calculate Jacobian matrix at element center
00841 
00842     // derivatives of bilinear map at quad center
00843     const double deriv_vals[] = { -0.5, -0.5, 0.5, -0.5, 0.5, 0.5, -0.5, 0.5 };
00844     MsqMatrix< 4, 2 > coeff_derivs( deriv_vals );
00845     MsqMatrix< 4, 3 > coords( vertices );
00846     MsqMatrix< 3, 2 > J = transpose( coords ) * coeff_derivs;
00847     // calculate expected metric value
00848     const double expt_val = sqr_Frobenius( J );
00849     // calculate derivative for each element vertex
00850     MsqVector< 3 > expt_grad[4];
00851     for( int v = 0; v < 4; ++v )
00852         expt_grad[v] = 2 * J * transpose( coeff_derivs.row( v ) );
00853 
00854     // construct metric
00855     pd.attach_settings( &settings );
00856     TestGradTargetMetric< typename TMPTypes< QMType >::MetricType > tm;
00857     // IdealShapeTarget tc;
00858     QMType m( tc, &tm );
00859     PlanarDomain plane( PlanarDomain::XY );
00860     pd.set_domain( &plane );
00861 
00862     // evaluate metric
00863     double act_val;
00864     std::vector< size_t > indices;
00865     std::vector< Vector3D > act_grad;
00866     size_t h = ElemSampleQM::handle( Sample( 2, 0 ), 0 );
00867     m.evaluate_with_gradient( pd, h, act_val, indices, act_grad, err );
00868     ASSERT_NO_ERROR( err );
00869 
00870     // compare values
00871     CPPUNIT_ASSERT_DOUBLES_EQUAL( expt_val, act_val, 1e-10 );
00872     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[0]].data() ), act_grad[0], 1e-10 );
00873     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[1]].data() ), act_grad[1], 1e-10 );
00874     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[2]].data() ), act_grad[2], 1e-10 );
00875     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[3]].data() ), act_grad[3], 1e-10 );
00876 
00877     // check numerical approx of gradient
00878     m.QualityMetric::evaluate_with_gradient( pd, h, act_val, indices, act_grad, err );
00879     ASSERT_NO_ERROR( err );
00880     CPPUNIT_ASSERT_DOUBLES_EQUAL( expt_val, act_val, 1e-10 );
00881     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[0]].data() ), act_grad[0], 1e-5 );
00882     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[1]].data() ), act_grad[1], 1e-5 );
00883     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[2]].data() ), act_grad[2], 1e-5 );
00884     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[3]].data() ), act_grad[3], 1e-5 );
00885 }
00886 
00887 template < class QMType >
00888 inline void TMPQualityMetricTest< QMType >::test_gradient_3D()
00889 {
00890     MsqPrintError err( std::cout );
00891 
00892     // check for expected value at center of flattened hex
00893 
00894     // construct flattened hex
00895     const double z          = 0.5;
00896     const double vertices[] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0,
00897                                 0.0, 0.0, z,   1.0, 0.0, z,   1.0, 1.0, z,   0.0, 1.0, z };
00898     size_t conn[8]          = { 0, 1, 2, 3, 4, 5, 6, 7 };
00899     PatchData pd;
00900     pd.fill( 8, vertices, 1, HEXAHEDRON, conn, 0, err );
00901     ASSERT_NO_ERROR( err );
00902 
00903     // calculate Jacobian matrix at element center
00904 
00905     // derivatives of trilinear map at hex center
00906     const double deriv_vals[8][3] = { { -0.25, -0.25, -0.25 }, { 0.25, -0.25, -0.25 }, { 0.25, 0.25, -0.25 },
00907                                       { -0.25, 0.25, -0.25 },  { -0.25, -0.25, 0.25 }, { 0.25, -0.25, 0.25 },
00908                                       { 0.25, 0.25, 0.25 },    { -0.25, 0.25, 0.25 } };
00909     MsqMatrix< 8, 3 > coeff_derivs( deriv_vals );
00910     MsqMatrix< 8, 3 > coords( vertices );
00911     MsqMatrix< 3, 3 > J = transpose( coords ) * coeff_derivs;
00912     // calculate expected metric value
00913     const double expt_val = sqr_Frobenius( J );
00914     // calculate derivative for each element vertex
00915     MsqVector< 3 > expt_grad[8];
00916     for( int v = 0; v < 8; ++v )
00917         expt_grad[v] = 2 * J * transpose( coeff_derivs.row( v ) );
00918 
00919     // construct metric
00920     pd.attach_settings( &settings );
00921     TestGradTargetMetric< typename TMPTypes< QMType >::MetricType > tm;
00922     IdealShapeTarget tc;
00923     QMType m( &tc, 0, &tm );
00924 
00925     // evaluate metric
00926     double act_val;
00927     std::vector< size_t > indices;
00928     std::vector< Vector3D > act_grad;
00929     size_t h = ElemSampleQM::handle( Sample( 3, 0 ), 0 );
00930     m.evaluate_with_gradient( pd, h, act_val, indices, act_grad, err );
00931     ASSERT_NO_ERROR( err );
00932 
00933     // compare values
00934     CPPUNIT_ASSERT_DOUBLES_EQUAL( expt_val, act_val, 1e-10 );
00935     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[0]].data() ), act_grad[0], 1e-10 );
00936     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[1]].data() ), act_grad[1], 1e-10 );
00937     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[2]].data() ), act_grad[2], 1e-10 );
00938     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[3]].data() ), act_grad[3], 1e-10 );
00939     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[4]].data() ), act_grad[4], 1e-10 );
00940     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[5]].data() ), act_grad[5], 1e-10 );
00941     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[6]].data() ), act_grad[6], 1e-10 );
00942     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[7]].data() ), act_grad[7], 1e-10 );
00943 
00944     // check numerical approx of gradient
00945     m.QualityMetric::evaluate_with_gradient( pd, h, act_val, indices, act_grad, err );
00946     ASSERT_NO_ERROR( err );
00947     CPPUNIT_ASSERT_DOUBLES_EQUAL( expt_val, act_val, 1e-10 );
00948     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[0]].data() ), act_grad[0], 1e-5 );
00949     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[1]].data() ), act_grad[1], 1e-5 );
00950     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[2]].data() ), act_grad[2], 1e-5 );
00951     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[3]].data() ), act_grad[3], 1e-5 );
00952     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[4]].data() ), act_grad[4], 1e-5 );
00953     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[5]].data() ), act_grad[5], 1e-5 );
00954     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[6]].data() ), act_grad[6], 1e-5 );
00955     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( expt_grad[indices[7]].data() ), act_grad[7], 1e-5 );
00956 }
00957 
00958 template < class QMType >
00959 inline void TMPQualityMetricTest< QMType >::compare_analytical_and_numerical_gradients( QualityMetric* qm,
00960                                                                                         PatchData& pd, int dim )
00961 {
00962     MsqPrintError err( std::cout );
00963 
00964     std::vector< size_t > handles, indices1, indices2;
00965     std::vector< Vector3D > grad1, grad2;
00966     double qm_val1, qm_val2;
00967     bool rval;
00968 
00969     qm->get_evaluations( pd, handles, false, err );
00970     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00971     CPPUNIT_ASSERT( !handles.empty() );
00972     for( size_t j = 0; j < handles.size(); ++j )
00973     {
00974         rval = qm->QualityMetric::evaluate_with_gradient( pd, handles[j], qm_val1, indices1, grad1, err );
00975         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00976         CPPUNIT_ASSERT( rval );
00977 
00978         // For analytical gradient of a 2D element in the XY plane,
00979         // we expect all Z terms to be zero.
00980         if( dim == 2 )
00981             for( size_t k = 0; k < grad1.size(); ++k )
00982                 grad1[k][2] = 0.0;
00983 
00984         rval = qm->evaluate_with_gradient( pd, handles[j], qm_val2, indices2, grad2, err );
00985         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00986         CPPUNIT_ASSERT( rval );
00987 
00988         CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
00989         CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
00990         CPPUNIT_ASSERT( !indices1.empty() );
00991 
00992         std::vector< size_t >::iterator it1, it2;
00993         for( it1 = indices1.begin(); it1 != indices1.end(); ++it1 )
00994         {
00995             it2 = std::find( indices2.begin(), indices2.end(), *it1 );
00996             CPPUNIT_ASSERT( it2 != indices2.end() );
00997 
00998             size_t idx1 = it1 - indices1.begin();
00999             size_t idx2 = it2 - indices2.begin();
01000             CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[idx1], grad2[idx2], 0.01 );
01001         }
01002     }
01003 }
01004 
01005 // Delcare specialized versions of the functions from
01006 // QualityMetricTester because we surface elements must
01007 // be handled differently.  For a surface element in the XY plane,
01008 // the finite difference approximations of the derivatives will
01009 // have non-zero values for derivatives wrt Z coordinates while the
01010 // analytical derivative calculations will return all derivatives
01011 // wrt Z coordiantes as zero.
01012 template < class QMType >
01013 inline void TMPQualityMetricTest< QMType >::compare_analytical_and_numerical_hessians( QualityMetric* qm )
01014 {
01015     MsqPrintError err( std::cout );
01016     PatchData pd;
01017     const EntityTopology types[] = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON };
01018     const int num_types          = sizeof( types ) / sizeof( types[0] );
01019     for( int i = 0; i < num_types; ++i )
01020     {
01021         get_nonideal_element( types[i], pd );
01022 
01023         std::vector< size_t > handles, indices1, indices2;
01024         std::vector< Vector3D > grad1, grad2;
01025         std::vector< Matrix3D > Hess1, Hess2;
01026         double qm_val1, qm_val2;
01027         bool rval;
01028 
01029         qm->get_evaluations( pd, handles, false, err );
01030         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01031         CPPUNIT_ASSERT( !handles.empty() );
01032         for( size_t j = 0; j < handles.size(); ++j )
01033         {
01034             rval = qm->QualityMetric::evaluate_with_Hessian( pd, handles[j], qm_val1, indices1, grad1, Hess1, err );
01035             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01036             CPPUNIT_ASSERT( rval );
01037 
01038             // For analytical gradient of a 2D element in the XY plane,
01039             // we expect all Z terms to be zero.
01040 #ifdef PLANAR_HESSIAN
01041             if( TopologyInfo::dimension( types[i] ) == 2 )
01042                 for( size_t k = 0; k < Hess1.size(); ++k )
01043                     Hess1[k]( 0, 2 ) = Hess1[k]( 1, 2 ) = Hess1[k]( 2, 0 ) = Hess1[k]( 2, 1 ) = Hess1[k]( 2, 2 ) = 0.0;
01044 #endif
01045 
01046             rval = qm->evaluate_with_Hessian( pd, handles[j], qm_val2, indices2, grad2, Hess2, err );
01047             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01048             CPPUNIT_ASSERT( rval );
01049 
01050             CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
01051             CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
01052             CPPUNIT_ASSERT( !indices1.empty() );
01053 
01054             std::vector< size_t >::iterator it;
01055             unsigned h = 0;
01056             for( unsigned r = 0; r < indices1.size(); ++r )
01057             {
01058                 it = std::find( indices2.begin(), indices2.end(), indices1[r] );
01059                 CPPUNIT_ASSERT( it != indices2.end() );
01060                 unsigned r2 = it - indices2.begin();
01061 
01062                 for( unsigned c = r; c < indices1.size(); ++c, ++h )
01063                 {
01064                     it = std::find( indices2.begin(), indices2.end(), indices1[c] );
01065                     CPPUNIT_ASSERT( it != indices2.end() );
01066                     unsigned c2 = it - indices2.begin();
01067 
01068                     unsigned h2;
01069                     if( r2 <= c2 )
01070                         h2 = indices2.size() * r - r * ( r + 1 ) / 2 + c;
01071                     else
01072                         h2 = indices2.size() * c - c * ( c + 1 ) / 2 + r;
01073 
01074                     // if (!utest_mat_equal(Hess1[h],Hess2[h2],0.001))
01075                     //  assert(false);
01076                     CPPUNIT_ASSERT_MATRICES_EQUAL( Hess1[h], Hess2[h2], 0.05 );
01077                 }
01078             }
01079         }
01080     }
01081 }
01082 
01083 // Delcare specialized versions of the functions from
01084 // QualityMetricTester because we surface elements must
01085 // be handled differently.  For a surface element in the XY plane,
01086 // the finite difference approximations of the derivatives will
01087 // have non-zero values for derivatives wrt Z coordinates while the
01088 // analytical derivative calculations will return all derivatives
01089 // wrt Z coordiantes as zero.
01090 template < class QMType >
01091 inline void TMPQualityMetricTest< QMType >::compare_analytical_and_numerical_diagonals( QualityMetric* qm )
01092 {
01093     MsqPrintError err( std::cout );
01094     PatchData pd;
01095     const EntityTopology types[] = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON };
01096     const int num_types          = sizeof( types ) / sizeof( types[0] );
01097     for( int i = 0; i < num_types; ++i )
01098     {
01099         get_nonideal_element( types[i], pd );
01100 
01101         std::vector< size_t > handles, indices1, indices2;
01102         std::vector< Vector3D > grad1, grad2;
01103         std::vector< Matrix3D > Hess1;
01104         std::vector< SymMatrix3D > Hess2;
01105         double qm_val1, qm_val2;
01106         bool rval;
01107 
01108         qm->get_evaluations( pd, handles, false, err );
01109         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01110         CPPUNIT_ASSERT( !handles.empty() );
01111         for( size_t j = 0; j < handles.size(); ++j )
01112         {
01113             rval = qm->QualityMetric::evaluate_with_Hessian( pd, handles[j], qm_val1, indices1, grad1, Hess1, err );
01114             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01115             CPPUNIT_ASSERT( rval );
01116 
01117             // For analytical gradient of a 2D element in the XY plane,
01118             // we expect all Z terms to be zero.
01119 #ifdef PLANAR_HESSIAN
01120             if( TopologyInfo::dimension( types[i] ) == 2 )
01121                 for( size_t k = 0; k < Hess1.size(); ++k )
01122                     Hess1[k]( 0, 2 ) = Hess1[k]( 1, 2 ) = Hess1[k]( 2, 0 ) = Hess1[k]( 2, 1 ) = Hess1[k]( 2, 2 ) = 0.0;
01123 #endif
01124 
01125             rval = qm->evaluate_with_Hessian_diagonal( pd, handles[j], qm_val2, indices2, grad2, Hess2, err );
01126             CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
01127             CPPUNIT_ASSERT( rval );
01128 
01129             CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
01130             CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
01131             CPPUNIT_ASSERT( !indices1.empty() );
01132             CPPUNIT_ASSERT_EQUAL( indices1.size() * ( indices1.size() + 1 ) / 2, Hess1.size() );
01133             CPPUNIT_ASSERT_EQUAL( indices2.size(), Hess2.size() );
01134 
01135             size_t h = 0;
01136             std::vector< size_t >::iterator it;
01137             for( unsigned r = 0; r < indices1.size(); ++r )
01138             {
01139                 it = std::find( indices2.begin(), indices2.end(), indices1[r] );
01140                 CPPUNIT_ASSERT( it != indices2.end() );
01141                 unsigned r2 = it - indices2.begin();
01142                 // if (!utest_mat_equal(Hess1[h],Hess2[r2],0.001))
01143                 //  assert(false);
01144                 CPPUNIT_ASSERT_MATRICES_EQUAL( Hess1[h], Hess2[r2], 0.05 );
01145                 h += indices1.size() - r;
01146             }
01147         }
01148     }
01149 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines