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