MOAB: Mesh Oriented datABase  (version 5.2.1)
TQualityMetricTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2010 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     (2010) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file TQualityMetricTest.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "TQualityMetric.hpp"
00034 #include "TInverseMeanRatio.hpp"
00035 #include "IdealWeightInverseMeanRatio.hpp"
00036 #include "TMPQualityMetricTest.hpp"
00037 #include "TShapeNB1.hpp"
00038 
00039 template <>
00040 class TMPTypes< TQualityMetric >
00041 {
00042   public:
00043     typedef TMetric MetricType;
00044     typedef TShapeNB1 TestType;
00045 };
00046 
00047 class TQualityMetricTest : public TMPQualityMetricTest< TQualityMetric >
00048 {
00049     CPPUNIT_TEST_SUITE( TQualityMetricTest );
00050 
00051     REGISTER_TMP_TESTS
00052 
00053     CPPUNIT_TEST( test_inverse_mean_ratio_grad );
00054     CPPUNIT_TEST( test_inverse_mean_ratio_hess );
00055     CPPUNIT_TEST( test_inverse_mean_ratio_hess_diag );
00056     CPPUNIT_TEST( regression_inverse_mean_ratio_grad );
00057     CPPUNIT_TEST( regression_inverse_mean_ratio_hess );
00058 
00059     CPPUNIT_TEST_SUITE_END();
00060 
00061   public:
00062     void test_inverse_mean_ratio_grad();
00063     void test_inverse_mean_ratio_hess();
00064     void test_inverse_mean_ratio_hess_diag();
00065     void regression_inverse_mean_ratio_grad();
00066     void regression_inverse_mean_ratio_hess();
00067 };
00068 
00069 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TQualityMetricTest, "TQualityMetricTest" );
00070 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TQualityMetricTest, "Unit" );
00071 
00072 void TQualityMetricTest::test_inverse_mean_ratio_grad()
00073 {
00074     TInverseMeanRatio tm;
00075     IdealShapeTarget target;
00076     TQualityMetric metric( &target, &tm );
00077     ElementPMeanP avg( 1.0, &metric );
00078 
00079     tester.test_gradient_reflects_quality( &metric );
00080     compare_analytical_and_numerical_gradients( &metric );
00081     tester.test_gradient_with_fixed_vertex( &avg );
00082 }
00083 
00084 void TQualityMetricTest::test_inverse_mean_ratio_hess()
00085 {
00086     TInverseMeanRatio tm;
00087     IdealShapeTarget target;
00088     TQualityMetric metric( &target, &tm );
00089     ElementPMeanP avg( 1.0, &metric );
00090 
00091     compare_analytical_and_numerical_hessians( &metric );
00092     tester.test_symmetric_Hessian_diagonal_blocks( &metric );
00093     tester.test_hessian_with_fixed_vertex( &avg );
00094 }
00095 
00096 void TQualityMetricTest::test_inverse_mean_ratio_hess_diag()
00097 {
00098     TInverseMeanRatio tm;
00099     IdealShapeTarget target;
00100     TQualityMetric metric( &target, &tm );
00101 
00102     compare_analytical_and_numerical_diagonals( &metric );
00103     tester.compare_eval_with_diag_and_eval_with_hessian( &metric );
00104 }
00105 
00106 void TQualityMetricTest::regression_inverse_mean_ratio_grad()
00107 {
00108     MsqError err;
00109     TInverseMeanRatio tm;
00110     IdealShapeTarget target;
00111     TQualityMetric metric( &target, &tm );
00112     const double coords[]  = { -0.80000000000000004, -0.80000000000000004, 0,
00113                               0.00000000000000000,  2.00000000000000000,  0,
00114                               -1.73205079999999990, 1.00000000000000000,  0 };
00115     const size_t indices[] = { 0, 1, 2 };
00116     PatchData pd;
00117     pd.fill( 3, coords, 1, TRIANGLE, indices, 0, err );
00118     pd.attach_settings( &settings );
00119     PlanarDomain dom( PlanarDomain::XY, coords[0] );
00120     pd.set_domain( &dom );
00121 
00122     IdealWeightInverseMeanRatio ref_metric;
00123 
00124     double exp_val, act_val;
00125     std::vector< size_t > exp_idx, act_idx, handles;
00126     std::vector< Vector3D > exp_grad, act_grad;
00127 
00128     handles.clear();
00129     ref_metric.get_evaluations( pd, handles, false, err );
00130     ASSERT_NO_ERROR( err );
00131     CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
00132     const size_t hand1 = handles.front();
00133     handles.clear();
00134     metric.get_evaluations( pd, handles, false, err );
00135     ASSERT_NO_ERROR( err );
00136     CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
00137     const size_t hand2 = handles.front();
00138 
00139     bool exp_rval, act_rval;
00140     exp_rval = ref_metric.evaluate_with_gradient( pd, hand1, exp_val, exp_idx, exp_grad, err );
00141     ASSERT_NO_ERROR( err );
00142     act_rval = metric.evaluate_with_gradient( pd, hand2, act_val, act_idx, act_grad, err );
00143     ASSERT_NO_ERROR( err );
00144 
00145     CPPUNIT_ASSERT( exp_rval );
00146     CPPUNIT_ASSERT( act_rval );
00147     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_val - 1.0, act_val, 1e-5 );
00148     CPPUNIT_ASSERT_EQUAL( (size_t)3, exp_idx.size() );
00149     CPPUNIT_ASSERT_EQUAL( (size_t)3, act_idx.size() );
00150 
00151     std::vector< size_t > sorted( exp_idx );
00152     std::sort( sorted.begin(), sorted.end() );
00153     CPPUNIT_ASSERT_EQUAL( (size_t)0, sorted[0] );
00154     CPPUNIT_ASSERT_EQUAL( (size_t)1, sorted[1] );
00155     CPPUNIT_ASSERT_EQUAL( (size_t)2, sorted[2] );
00156 
00157     sorted = act_idx;
00158     std::sort( sorted.begin(), sorted.end() );
00159     CPPUNIT_ASSERT_EQUAL( (size_t)0, sorted[0] );
00160     CPPUNIT_ASSERT_EQUAL( (size_t)1, sorted[1] );
00161     CPPUNIT_ASSERT_EQUAL( (size_t)2, sorted[2] );
00162 
00163     const size_t idx_map[] = { std::find( act_idx.begin(), act_idx.end(), exp_idx[0] ) - act_idx.begin(),
00164                                std::find( act_idx.begin(), act_idx.end(), exp_idx[1] ) - act_idx.begin(),
00165                                std::find( act_idx.begin(), act_idx.end(), exp_idx[2] ) - act_idx.begin() };
00166     CPPUNIT_ASSERT_VECTORS_EQUAL( exp_grad[0], act_grad[idx_map[0]], 1e-5 );
00167     CPPUNIT_ASSERT_VECTORS_EQUAL( exp_grad[1], act_grad[idx_map[1]], 1e-5 );
00168     CPPUNIT_ASSERT_VECTORS_EQUAL( exp_grad[2], act_grad[idx_map[2]], 1e-5 );
00169 }
00170 
00171 void TQualityMetricTest::regression_inverse_mean_ratio_hess()
00172 {
00173     MsqError err;
00174     TInverseMeanRatio tm;
00175     IdealShapeTarget target;
00176     TQualityMetric metric( &target, &tm );
00177     const double coords[] = { 4.158984727, 4.6570859130000004, 5, 4.51742825, 4.51742825, 5, 4.3103448279999999, 5, 5 };
00178     const bool fixed[]    = { false, false, true };
00179     const size_t indices[] = { 0, 1, 2 };
00180     PatchData pd;
00181     pd.fill( 3, coords, 1, TRIANGLE, indices, fixed, err );
00182     pd.attach_settings( &settings );
00183     PlanarDomain dom( PlanarDomain::XY, coords[2] );
00184     pd.set_domain( &dom );
00185 
00186     IdealWeightInverseMeanRatio ref_metric;
00187 
00188     double exp_val, act_val;
00189     std::vector< size_t > exp_idx, act_idx, handles;
00190     std::vector< Vector3D > exp_grad, act_grad;
00191     std::vector< Matrix3D > exp_hess, act_hess;
00192 
00193     handles.clear();
00194     ref_metric.get_evaluations( pd, handles, false, err );
00195     ASSERT_NO_ERROR( err );
00196     CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
00197     const size_t hand1 = handles.front();
00198     handles.clear();
00199     metric.get_evaluations( pd, handles, false, err );
00200     ASSERT_NO_ERROR( err );
00201     CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
00202     const size_t hand2 = handles.front();
00203 
00204     // first make sure that non-TMP IMR metric works correctly
00205     bool exp_rval, act_rval;
00206     exp_rval = ref_metric.evaluate_with_Hessian( pd, hand1, exp_val, exp_idx, exp_grad, exp_hess, err );
00207     ASSERT_NO_ERROR( err );
00208     act_rval = ref_metric.QualityMetric::evaluate_with_Hessian( pd, hand1, act_val, act_idx, act_grad, act_hess, err );
00209     CPPUNIT_ASSERT( exp_rval );
00210     CPPUNIT_ASSERT( act_rval );
00211     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_val, act_val, 1e-5 );
00212     CPPUNIT_ASSERT_EQUAL( (size_t)2, exp_idx.size() );
00213     CPPUNIT_ASSERT_EQUAL( (size_t)2, act_idx.size() );
00214 
00215     if( act_idx[0] == exp_idx[0] )
00216     {
00217         CPPUNIT_ASSERT_EQUAL( exp_idx[1], act_idx[1] );
00218         CPPUNIT_ASSERT_VECTORS_EQUAL( exp_grad[0], act_grad[0], 1e-5 );
00219         CPPUNIT_ASSERT_VECTORS_EQUAL( exp_grad[1], act_grad[1], 1e-5 );
00220         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[0], act_hess[0], 5e-3 );
00221         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[1], act_hess[1], 5e-3 );
00222         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[2], act_hess[2], 5e-3 );
00223     }
00224     else
00225     {
00226         CPPUNIT_ASSERT_EQUAL( exp_idx[0], act_idx[1] );
00227         CPPUNIT_ASSERT_EQUAL( exp_idx[1], act_idx[0] );
00228         CPPUNIT_ASSERT_VECTORS_EQUAL( exp_grad[0], act_grad[1], 1e-5 );
00229         CPPUNIT_ASSERT_VECTORS_EQUAL( exp_grad[1], act_grad[0], 1e-5 );
00230         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[0], act_hess[2], 5e-3 );
00231         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[1], transpose( act_hess[1] ), 5e-3 );
00232         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[2], act_hess[0], 5e-3 );
00233     }
00234 
00235     // now compare TMP metric with non-TMP metric
00236     act_rval = metric.evaluate_with_Hessian( pd, hand2, act_val, act_idx, act_grad, act_hess, err );
00237     ASSERT_NO_ERROR( err );
00238     CPPUNIT_ASSERT( act_rval );
00239     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_val - 1.0, act_val, 1e-5 );
00240     CPPUNIT_ASSERT_EQUAL( (size_t)2, act_idx.size() );
00241 
00242 #ifdef PLANAR_HESSIAN
00243     // zero derivatives with respect to Z
00244     for( int i = 0; i < 2; ++i )
00245         exp_grad[i][2] = 0.0;
00246     for( int i = 0; i < 3; ++i )
00247     {
00248         for( int j = 0; j < 3; ++j )
00249             exp_hess[i][j][2] = exp_hess[i][2][j] = 0.0;
00250     }
00251 #else
00252     // don't compare double out-of-plane term because metrics
00253     // make varying assumptions about behavior of Hessian
00254     for( int i = 0; i < 3; ++i )
00255         exp_hess[i][2][2] = act_hess[i][2][2] = 0.0;
00256 #endif
00257 
00258     if( act_idx[0] == exp_idx[0] )
00259     {
00260         CPPUNIT_ASSERT_EQUAL( exp_idx[1], act_idx[1] );
00261         CPPUNIT_ASSERT_VECTORS_EQUAL( exp_grad[0], act_grad[0], 1e-5 );
00262         CPPUNIT_ASSERT_VECTORS_EQUAL( exp_grad[1], act_grad[1], 1e-5 );
00263         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[0], act_hess[0], 5e-3 );
00264         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[1], act_hess[1], 5e-3 );
00265         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[2], act_hess[2], 5e-3 );
00266     }
00267     else
00268     {
00269         CPPUNIT_ASSERT_EQUAL( exp_idx[0], act_idx[1] );
00270         CPPUNIT_ASSERT_EQUAL( exp_idx[1], act_idx[0] );
00271         CPPUNIT_ASSERT_VECTORS_EQUAL( exp_grad[0], act_grad[1], 1e-5 );
00272         CPPUNIT_ASSERT_VECTORS_EQUAL( exp_grad[1], act_grad[0], 1e-5 );
00273         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[0], act_hess[2], 5e-3 );
00274         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[1], transpose( act_hess[1] ), 5e-3 );
00275         CPPUNIT_ASSERT_MATRICES_EQUAL( exp_hess[2], act_hess[0], 5e-3 );
00276     }
00277 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines