MOAB: Mesh Oriented datABase  (version 5.4.0)
PMeanPMetricTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     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 
00025   ***************************************************************** */
00026 
00027 #include "ElementPMeanP.hpp"
00028 #include "VertexPMeanP.hpp"
00029 
00030 #include "cppunit/extensions/HelperMacros.h"
00031 #include "PatchDataInstances.hpp"
00032 #include "UnitUtil.hpp"
00033 #include "ElemSampleQM.hpp"
00034 
00035 using namespace MBMesquite;
00036 using std::cerr;
00037 using std::cout;
00038 using std::endl;
00039 
00040 class PMeanPMetricTest : public CppUnit::TestFixture
00041 {
00042 
00043   private:
00044     CPPUNIT_TEST_SUITE( PMeanPMetricTest );
00045     CPPUNIT_TEST( test_get_metric_type );
00046     CPPUNIT_TEST( test_get_element_evaluations );
00047     CPPUNIT_TEST( test_get_vertex_evaluations );
00048     CPPUNIT_TEST( test_element_evaluate );
00049     CPPUNIT_TEST( test_vertex_evaluate );
00050     CPPUNIT_TEST( test_indices );
00051     CPPUNIT_TEST( test_gradient );
00052     CPPUNIT_TEST( test_hessian );
00053     CPPUNIT_TEST( test_hessian_diagonal );
00054     CPPUNIT_TEST_SUITE_END();
00055 
00056     PatchData pd;
00057 
00058   public:
00059     void setUp();
00060 
00061     void test_get_metric_type();
00062     void test_get_element_evaluations();
00063     void test_get_vertex_evaluations();
00064     void test_element_evaluate();
00065     void test_vertex_evaluate();
00066     void test_indices();
00067     void test_gradient();
00068     void test_hessian();
00069     void test_hessian_diagonal();
00070 };
00071 
00072 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( PMeanPMetricTest, "PMeanPMetricTest" );
00073 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( PMeanPMetricTest, "Unit" );
00074 
00075 void PMeanPMetricTest::setUp()
00076 {
00077     MsqError err;
00078     create_four_quads_patch( pd, err );
00079     CPPUNIT_ASSERT( !err );
00080 }
00081 
00082 // bogus metric for testing.
00083 // returns vertex index value for metric value, {vertex index, 0, 1}
00084 // for gradint wrt vertex , and { (h1, h2, h1+j2), {h1, 1, h1+h2}, {2*h1, 2*h2, h2} }
00085 // for hessin where h1 and h2 are vertex indices.
00086 class FauxMetric : public ElemSampleQM
00087 {
00088   public:
00089     MetricType get_metric_type() const
00090     {
00091         return ELEMENT_BASED;
00092     }
00093     std::string get_name() const
00094     {
00095         return "FauxMetrix";
00096     }
00097     int get_negate_flag() const
00098     {
00099         return 1;
00100     }
00101     void get_evaluations( PatchData& pd, std::vector< size_t >& h, bool free, MsqError& );
00102     void get_element_evaluations( PatchData&, size_t, std::vector< size_t >&, MsqError& );
00103     bool evaluate( PatchData& pd, size_t h, double& v, MsqError& err );
00104     bool evaluate_with_indices( PatchData& pd, size_t h, double& v, std::vector< size_t >& indices, MsqError& err );
00105     bool evaluate_with_gradient( PatchData& pd,
00106                                  size_t h,
00107                                  double& v,
00108                                  std::vector< size_t >& indices,
00109                                  std::vector< Vector3D >& grads,
00110                                  MsqError& err );
00111     bool evaluate_with_Hessian( PatchData& pd,
00112                                 size_t h,
00113                                 double& v,
00114                                 std::vector< size_t >& indices,
00115                                 std::vector< Vector3D >& grad,
00116                                 std::vector< Matrix3D >& Hess,
00117                                 MsqError& err );
00118 };
00119 
00120 void FauxMetric::get_evaluations( PatchData& pd, std::vector< size_t >& h, bool free, MsqError& err )
00121 {
00122     h.clear();
00123     for( size_t i = 0; i < pd.num_elements(); ++i )
00124         get_element_evaluations( pd, i, h, err );
00125 }
00126 
00127 void FauxMetric::get_element_evaluations( PatchData& pd, size_t h, std::vector< size_t >& list, MsqError& )
00128 {
00129     MsqMeshEntity& elem = pd.element_by_index( h );
00130     for( unsigned i = 0; i < elem.corner_count(); ++i )
00131         list.push_back( handle( Sample( 0, i ), h ) );
00132 }
00133 
00134 bool FauxMetric::evaluate( PatchData& pd, size_t h, double& v, MsqError& )
00135 {
00136     size_t e      = ElemSampleQM::elem( h );
00137     Sample s      = ElemSampleQM::sample( h );
00138     size_t* verts = pd.element_by_index( e ).get_vertex_index_array();
00139     CPPUNIT_ASSERT_EQUAL( (unsigned short)0, s.dimension );
00140     v = (double)( verts[s.number] );
00141     return true;
00142 }
00143 
00144 bool FauxMetric::evaluate_with_indices( PatchData& pd,
00145                                         size_t h,
00146                                         double& v,
00147                                         std::vector< size_t >& indices,
00148                                         MsqError& err )
00149 {
00150     evaluate( pd, h, v, err );
00151     indices.resize( 3 );
00152     size_t e      = ElemSampleQM::elem( h );
00153     Sample s      = ElemSampleQM::sample( h );
00154     size_t* verts = pd.element_by_index( e ).get_vertex_index_array();
00155     size_t n      = pd.element_by_index( e ).vertex_count();
00156     CPPUNIT_ASSERT_EQUAL( (unsigned short)0, s.dimension );
00157     indices[0] = verts[s.number];
00158     indices[1] = verts[( s.number + 1 ) % n];
00159     indices[2] = verts[( s.number + n - 1 ) % n];
00160     return true;
00161 }
00162 
00163 bool FauxMetric::evaluate_with_gradient( PatchData& pd,
00164                                          size_t h,
00165                                          double& v,
00166                                          std::vector< size_t >& indices,
00167                                          std::vector< Vector3D >& grads,
00168                                          MsqError& err )
00169 {
00170     evaluate_with_indices( pd, h, v, indices, err );
00171     grads.clear();
00172     for( unsigned i = 0; i < indices.size(); ++i )
00173         grads.push_back( Vector3D( (double)indices[i], 0, 1 ) );
00174     return true;
00175 }
00176 
00177 bool FauxMetric::evaluate_with_Hessian( PatchData& pd,
00178                                         size_t h,
00179                                         double& v,
00180                                         std::vector< size_t >& indices,
00181                                         std::vector< Vector3D >& grad,
00182                                         std::vector< Matrix3D >& hess,
00183                                         MsqError& err )
00184 {
00185     evaluate_with_gradient( pd, h, v, indices, grad, err );
00186     hess.clear();
00187     for( unsigned r = 0; r < indices.size(); ++r )
00188         for( unsigned c = r; c < indices.size(); ++c )
00189             hess.push_back( Matrix3D( indices[r], indices[c], indices[r] + indices[c], indices[r], 1.0,
00190                                       indices[r] + indices[c], 2 * indices[r], 2 * indices[c], indices[c] ) );
00191     return true;
00192 }
00193 
00194 void PMeanPMetricTest::test_get_metric_type()
00195 {
00196     FauxMetric m;
00197     ElementPMeanP e( 1.0, &m );
00198     VertexPMeanP v( 1.0, &m );
00199     CPPUNIT_ASSERT( QualityMetric::ELEMENT_BASED == e.get_metric_type() );
00200     CPPUNIT_ASSERT( QualityMetric::VERTEX_BASED == v.get_metric_type() );
00201 }
00202 
00203 void PMeanPMetricTest::test_get_element_evaluations()
00204 {
00205     MsqError err;
00206     FauxMetric m;
00207     ElementPMeanP e( 1.0, &m );
00208     std::vector< size_t > handles;
00209     // test that handles array contains all elements
00210     e.get_evaluations( pd, handles, false, err );
00211     std::sort( handles.begin(), handles.end() );
00212     CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
00213     for( unsigned i = 1; i < handles.size(); ++i )
00214         CPPUNIT_ASSERT_EQUAL( handles[i - 1] + 1, handles[i] );
00215     // test that handles array contains all elements
00216     e.get_evaluations( pd, handles, true, err );
00217     std::sort( handles.begin(), handles.end() );
00218     CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
00219     for( unsigned i = 1; i < handles.size(); ++i )
00220         CPPUNIT_ASSERT_EQUAL( handles[i - 1] + 1, handles[i] );
00221 }
00222 
00223 void PMeanPMetricTest::test_get_vertex_evaluations()
00224 {
00225     MsqError err;
00226     FauxMetric m;
00227     VertexPMeanP e( 1.0, &m );
00228     std::vector< size_t > handles;
00229     // test that handles array contains all vertices
00230     e.get_evaluations( pd, handles, false, err );
00231     std::sort( handles.begin(), handles.end() );
00232     CPPUNIT_ASSERT_EQUAL( pd.num_nodes(), handles.size() );
00233     for( unsigned i = 1; i < handles.size(); ++i )
00234         CPPUNIT_ASSERT_EQUAL( handles[i - 1] + 1, handles[i] );
00235     // test that handles array contains all vertices
00236     e.get_evaluations( pd, handles, true, err );
00237     std::sort( handles.begin(), handles.end() );
00238     CPPUNIT_ASSERT_EQUAL( pd.num_nodes(), handles.size() );
00239     for( unsigned i = 1; i < handles.size(); ++i )
00240         CPPUNIT_ASSERT_EQUAL( handles[i - 1] + 1, handles[i] );
00241 }
00242 
00243 void PMeanPMetricTest::test_vertex_evaluate()
00244 {
00245     MsqError err;
00246     FauxMetric m;
00247     VertexPMeanP m1( 1.0, &m );
00248     VertexPMeanP m2( 2.0, &m );
00249 
00250     // evaluate average around vertex
00251     double v1, v2;
00252     m1.evaluate( pd, 0, v1, err );
00253     CPPUNIT_ASSERT( !err );
00254     m2.evaluate( pd, 0, v2, err );
00255     CPPUNIT_ASSERT( !err );
00256 
00257     // get elements adjacent to vertex
00258     size_t num_elem;
00259     const size_t* elems = pd.get_vertex_element_adjacencies( 0, num_elem, err );
00260     CPPUNIT_ASSERT( !err );
00261 
00262     // calculate expected values from underyling metric
00263     double ev1 = 0.0, ev2 = 0.0;
00264     for( unsigned i = 0; i < num_elem; ++i )
00265     {
00266         const MsqMeshEntity& elem = pd.element_by_index( elems[i] );
00267         const size_t* verts       = elem.get_vertex_index_array();
00268         const size_t* end         = verts + elem.node_count();
00269         const size_t* p           = std::find( verts, end, (size_t)0 );
00270         CPPUNIT_ASSERT( p < end );
00271         size_t h = ElemSampleQM::handle( Sample( 0, p - verts ), elems[i] );
00272 
00273         double v;
00274         m.evaluate( pd, h, v, err );
00275         CPPUNIT_ASSERT( !err );
00276         ev1 += v;
00277         ev2 += v * v;
00278     }
00279 
00280     ev1 /= num_elem;
00281     ev2 /= num_elem;
00282 
00283     CPPUNIT_ASSERT_DOUBLES_EQUAL( ev1, v1, 1e-6 );
00284     CPPUNIT_ASSERT_DOUBLES_EQUAL( ev2, v2, 1e-6 );
00285 }
00286 
00287 void PMeanPMetricTest::test_element_evaluate()
00288 {
00289     MsqError err;
00290     FauxMetric m;
00291     ElementPMeanP m1( 1.0, &m );
00292     ElementPMeanP m2( 2.0, &m );
00293 
00294     // evaluate average over element
00295     double v1, v2;
00296     m1.evaluate( pd, 0, v1, err );
00297     CPPUNIT_ASSERT( !err );
00298     m2.evaluate( pd, 0, v2, err );
00299     CPPUNIT_ASSERT( !err );
00300 
00301     // get sample points within element
00302     std::vector< size_t > handles;
00303     m.get_element_evaluations( pd, 0, handles, err );
00304     CPPUNIT_ASSERT( !err );
00305 
00306     // calculate expected values from underyling metric
00307     double ev1 = 0.0, ev2 = 0.0;
00308     for( unsigned i = 0; i < handles.size(); ++i )
00309     {
00310         double v;
00311         m.evaluate( pd, handles[i], v, err );
00312         CPPUNIT_ASSERT( !err );
00313         ev1 += v;
00314         ev2 += v * v;
00315     }
00316     ev1 /= handles.size();
00317     ev2 /= handles.size();
00318 
00319     CPPUNIT_ASSERT_DOUBLES_EQUAL( ev1, v1, 1e-6 );
00320     CPPUNIT_ASSERT_DOUBLES_EQUAL( ev2, v2, 1e-6 );
00321 }
00322 
00323 void PMeanPMetricTest::test_indices()
00324 {
00325     MsqError err;
00326     FauxMetric m;
00327     ElementPMeanP m2( 2.0, &m );
00328 
00329     double v1, v2;
00330     std::vector< size_t > indices;
00331     m2.evaluate_with_indices( pd, 0, v1, indices, err );
00332     CPPUNIT_ASSERT( !err );
00333     m2.evaluate( pd, 0, v2, err );
00334     CPPUNIT_ASSERT( !err );
00335     CPPUNIT_ASSERT_DOUBLES_EQUAL( v2, v1, 1e-6 );
00336 
00337     std::vector< size_t > vertices;
00338     pd.element_by_index( 0 ).get_vertex_indices( vertices );
00339 
00340     std::sort( vertices.begin(), vertices.end() );
00341     std::sort( indices.begin(), indices.end() );
00342     CPPUNIT_ASSERT( vertices == indices );
00343 }
00344 
00345 template < typename T >
00346 size_t index_of( const std::vector< T >& v, T a )
00347 {
00348     return std::find( v.begin(), v.end(), a ) - v.begin();
00349 }
00350 
00351 void PMeanPMetricTest::test_gradient()
00352 {
00353     MsqError err;
00354     FauxMetric m;
00355     ElementPMeanP m1( 1.0, &m );
00356     ElementPMeanP m2( 2.0, &m );
00357 
00358     // get vertices for later
00359     std::vector< size_t > vertices;
00360     pd.element_by_index( 0 ).get_vertex_indices( vertices );
00361 
00362     // evaluate without gradients
00363     double v1, v2, v3, v4;
00364     std::vector< size_t > indices1, indices2, indices3, indices4;
00365     std::vector< Vector3D > grads1, grads2;
00366     m1.evaluate_with_indices( pd, 0, v1, indices1, err );
00367     CPPUNIT_ASSERT( !err );
00368     m2.evaluate_with_indices( pd, 0, v2, indices2, err );
00369     CPPUNIT_ASSERT( !err );
00370 
00371     // evaluate with gradients
00372     m1.evaluate_with_gradient( pd, 0, v3, indices3, grads1, err );
00373     CPPUNIT_ASSERT( !err );
00374     m2.evaluate_with_gradient( pd, 0, v4, indices4, grads2, err );
00375     CPPUNIT_ASSERT( !err );
00376 
00377     // compare value and indices to eval w/out gradient
00378     CPPUNIT_ASSERT_DOUBLES_EQUAL( v1, v3, 1e-6 );
00379     CPPUNIT_ASSERT_DOUBLES_EQUAL( v2, v4, 1e-6 );
00380     std::sort( indices1.begin(), indices1.end() );
00381     std::vector< size_t > tm( indices3 );
00382     std::sort( tm.begin(), tm.end() );
00383     CPPUNIT_ASSERT( tm == indices1 );
00384     std::sort( indices2.begin(), indices2.end() );
00385     tm = indices4;
00386     std::sort( tm.begin(), tm.end() );
00387     CPPUNIT_ASSERT( tm == indices2 );
00388 
00389     // setup evaluation of underying metric
00390     std::vector< size_t > handles;
00391     m.get_element_evaluations( pd, 0, handles, err );
00392     CPPUNIT_ASSERT( !err );
00393 
00394     // calculate expected gradients
00395     std::vector< Vector3D > expected1, expected2, temp;
00396     expected1.resize( vertices.size(), Vector3D( 0, 0, 0 ) );
00397     expected2.resize( vertices.size(), Vector3D( 0, 0, 0 ) );
00398     for( unsigned i = 0; i < handles.size(); ++i )
00399     {
00400         double v;
00401         m.evaluate_with_gradient( pd, handles[i], v, indices3, temp, err );
00402         CPPUNIT_ASSERT( !err );
00403         for( unsigned k = 0; k < indices3.size(); ++k )
00404         {
00405             unsigned idx = index_of( vertices, indices3[k] );
00406             CPPUNIT_ASSERT( idx < vertices.size() );
00407             expected1[idx] += temp[k];
00408             expected2[idx] += 2 * v * temp[k];
00409         }
00410     }
00411     for( unsigned i = 0; i < vertices.size(); ++i )
00412     {
00413         expected1[i] /= handles.size();
00414         expected2[i] /= handles.size();
00415     }
00416 
00417     // compare gradients
00418     for( unsigned i = 0; i < indices1.size(); ++i )
00419     {
00420         unsigned k = index_of( vertices, indices1[i] );
00421         CPPUNIT_ASSERT( k < vertices.size() );
00422         CPPUNIT_ASSERT_VECTORS_EQUAL( expected1[k], grads1[i], 1e-6 );
00423         CPPUNIT_ASSERT_VECTORS_EQUAL( expected2[k], grads2[i], 1e-6 );
00424     }
00425 }
00426 
00427 void PMeanPMetricTest::test_hessian()
00428 {
00429     MsqError err;
00430     FauxMetric m;
00431     ElementPMeanP m1( 1.0, &m );
00432     ElementPMeanP m2( 2.0, &m );
00433 
00434     // get vertices for later
00435     std::vector< size_t > vertices;
00436     pd.element_by_index( 0 ).get_vertex_indices( vertices );
00437 
00438     // evaluate gradient
00439     double v1, v2, v3, v4;
00440     std::vector< size_t > indices1, indices2, indices3, indices4, tmpi;
00441     std::vector< Vector3D > grad1, grad2, grad3, grad4;
00442     std::vector< Matrix3D > hess3, hess4;
00443     m1.evaluate_with_gradient( pd, 0, v1, indices1, grad1, err );
00444     CPPUNIT_ASSERT( !err );
00445     m2.evaluate_with_gradient( pd, 0, v2, indices2, grad2, err );
00446     CPPUNIT_ASSERT( !err );
00447 
00448     // evaluate with Hessian
00449     m1.evaluate_with_Hessian( pd, 0, v3, indices3, grad3, hess3, err );
00450     CPPUNIT_ASSERT( !err );
00451     m2.evaluate_with_Hessian( pd, 0, v4, indices4, grad4, hess4, err );
00452     CPPUNIT_ASSERT( !err );
00453 
00454     // compare value and indices to eval w/out gradient
00455     CPPUNIT_ASSERT_DOUBLES_EQUAL( v1, v3, 1e-6 );
00456     CPPUNIT_ASSERT_DOUBLES_EQUAL( v2, v4, 1e-6 );
00457     // It isn't a requirement that the index order remain the same
00458     // for both eval_with_grad and eval_with_Hess, but assuming it
00459     // does simplifies a lot of stuff in this test.  Check that the
00460     // assumption remains valid.
00461     CPPUNIT_ASSERT( indices3 == indices1 );
00462     CPPUNIT_ASSERT( indices4 == indices2 );
00463     // It isn't a requirement that the index order remain the same
00464     // for any value of P, but assuming it does simplifies a lot
00465     // of stuff in this test, so check that the assumption is valid.
00466     CPPUNIT_ASSERT( indices1 == indices2 );
00467 
00468     // check that gradient values match
00469     for( size_t i = 0; i < indices1.size(); ++i )
00470     {
00471         CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[i], grad3[i], 1e-5 );
00472         CPPUNIT_ASSERT_VECTORS_EQUAL( grad2[i], grad4[i], 1e-5 );
00473     }
00474 
00475     // setup evaluation of underying metric
00476     std::vector< size_t > handles;
00477     m.get_element_evaluations( pd, 0, handles, err );
00478     CPPUNIT_ASSERT( !err );
00479 
00480     // calculate expected Hessians
00481     std::vector< Vector3D > g;
00482     std::vector< Matrix3D > expected1, expected2, h;
00483     std::vector< Matrix3D >::iterator h_iter;
00484     const unsigned N = vertices.size();
00485     expected1.resize( N * ( N + 1 ) / 2, Matrix3D( 0, 0, 0, 0, 0, 0, 0, 0, 0 ) );
00486     expected2 = expected1;
00487     Matrix3D outer;
00488     for( unsigned i = 0; i < handles.size(); ++i )
00489     {
00490         double v;
00491         m.evaluate_with_Hessian( pd, handles[i], v, tmpi, g, h, err );
00492         CPPUNIT_ASSERT( !err );
00493         h_iter = h.begin();
00494         for( unsigned r = 0; r < tmpi.size(); ++r )
00495         {
00496             unsigned R = index_of( vertices, tmpi[r] );
00497             CPPUNIT_ASSERT( R < N );
00498             for( unsigned c = r; c < tmpi.size(); ++c, ++h_iter )
00499             {
00500                 unsigned C = index_of( vertices, tmpi[c] );
00501                 CPPUNIT_ASSERT( C < N );
00502                 if( R <= C )
00503                 {
00504                     unsigned idx = N * R - R * ( R + 1 ) / 2 + C;
00505                     expected1[idx] += 1.0 / handles.size() * *h_iter;
00506                     expected2[idx] += 2.0 * v / handles.size() * *h_iter;
00507                     outer.outer_product( g[r], g[c] );
00508                     expected2[idx] += 2.0 / handles.size() * outer;
00509                 }
00510                 else
00511                 {
00512                     unsigned idx = N * C - C * ( C + 1 ) / 2 + R;
00513                     expected1[idx] += 1.0 / handles.size() * transpose( *h_iter );
00514                     expected2[idx] += 2.0 * v / handles.size() * transpose( *h_iter );
00515                     outer.outer_product( g[c], g[r] );
00516                     expected2[idx] += 2.0 / handles.size() * outer;
00517                 }
00518             }
00519         }
00520     }
00521 
00522     // compare Hessians
00523     unsigned H_idx = 0;
00524     for( unsigned R = 0; R < vertices.size(); ++R )
00525     {
00526         if( vertices[R] >= pd.num_free_vertices() ) continue;
00527         unsigned r = index_of( indices3, vertices[R] );
00528         CPPUNIT_ASSERT( r < indices3.size() );
00529         for( unsigned C = R; C < vertices.size(); ++C, ++H_idx )
00530         {
00531             if( vertices[C] >= pd.num_free_vertices() ) continue;
00532             unsigned c = index_of( indices3, vertices[C] );
00533             CPPUNIT_ASSERT( c < indices3.size() );
00534             if( r <= c )
00535             {
00536                 unsigned idx = indices3.size() * r - r * ( r + 1 ) / 2 + c;
00537                 CPPUNIT_ASSERT_MATRICES_EQUAL( expected1[H_idx], hess3[idx], 1e-5 );
00538                 CPPUNIT_ASSERT_MATRICES_EQUAL( expected2[H_idx], hess4[idx], 1e-5 );
00539             }
00540             else
00541             {
00542                 unsigned idx = indices3.size() * c - c * ( c + 1 ) / 2 + r;
00543                 CPPUNIT_ASSERT_MATRICES_EQUAL( transpose( expected1[H_idx] ), hess3[idx], 1e-5 );
00544                 CPPUNIT_ASSERT_MATRICES_EQUAL( transpose( expected2[H_idx] ), hess4[idx], 1e-5 );
00545             }
00546         }
00547     }
00548 }
00549 
00550 void PMeanPMetricTest::test_hessian_diagonal()
00551 {
00552     MsqError err;
00553     FauxMetric m;
00554     ElementPMeanP m1( 1.0, &m );
00555     ElementPMeanP m2( 2.0, &m );
00556 
00557     // we've already validated the Hessian results in the
00558     // previous test, so just check that the diagonal terms
00559     // match the terms for the full Hessian.
00560     std::vector< size_t > m1_indices_h, m1_indices_d, m2_indices_h, m2_indices_d;
00561     std::vector< Vector3D > m1_g_h, m1_g_d, m2_g_h, m2_g_d;
00562     std::vector< Matrix3D > m1_h_h, m2_h_h;
00563     std::vector< SymMatrix3D > m1_h_d, m2_h_d;
00564     double m1_v_h, m1_v_d, m2_v_h, m2_v_d;
00565     m1.evaluate_with_Hessian( pd, 0, m1_v_h, m1_indices_h, m1_g_h, m1_h_h, err );
00566     ASSERT_NO_ERROR( err );
00567     m2.evaluate_with_Hessian( pd, 0, m2_v_h, m2_indices_h, m2_g_h, m2_h_h, err );
00568     ASSERT_NO_ERROR( err );
00569     m1.evaluate_with_Hessian_diagonal( pd, 0, m1_v_d, m1_indices_d, m1_g_d, m1_h_d, err );
00570     ASSERT_NO_ERROR( err );
00571     m2.evaluate_with_Hessian_diagonal( pd, 0, m2_v_d, m2_indices_d, m2_g_d, m2_h_d, err );
00572     ASSERT_NO_ERROR( err );
00573 
00574     // compare values
00575     CPPUNIT_ASSERT_DOUBLES_EQUAL( m1_v_h, m1_v_d, 1e-6 );
00576     CPPUNIT_ASSERT_DOUBLES_EQUAL( m2_v_h, m2_v_d, 1e-6 );
00577 
00578     // Assume indices in same order because
00579     // it simplifiers later code in test
00580     CPPUNIT_ASSERT( m1_indices_h == m1_indices_d );
00581     CPPUNIT_ASSERT( m2_indices_h == m1_indices_d );
00582 
00583     // compare gradient values
00584     CPPUNIT_ASSERT_EQUAL( m1_indices_h.size(), m1_g_h.size() );
00585     CPPUNIT_ASSERT_EQUAL( m2_indices_h.size(), m2_g_h.size() );
00586     CPPUNIT_ASSERT_EQUAL( m1_indices_d.size(), m1_g_d.size() );
00587     CPPUNIT_ASSERT_EQUAL( m2_indices_d.size(), m2_g_d.size() );
00588     for( unsigned i = 0; i < m1_indices_h.size(); ++i )
00589     {
00590         CPPUNIT_ASSERT_VECTORS_EQUAL( m1_g_h[i], m1_g_d[i], 1e-6 );
00591         CPPUNIT_ASSERT_VECTORS_EQUAL( m2_g_h[i], m2_g_d[i], 1e-6 );
00592     }
00593 
00594     // compare hessian diagonal terms
00595     CPPUNIT_ASSERT_EQUAL( m1_indices_h.size() * ( m1_indices_h.size() + 1 ) / 2, m1_h_h.size() );
00596     CPPUNIT_ASSERT_EQUAL( m2_indices_h.size() * ( m2_indices_h.size() + 1 ) / 2, m2_h_h.size() );
00597     CPPUNIT_ASSERT_EQUAL( m1_indices_d.size(), m1_h_d.size() );
00598     CPPUNIT_ASSERT_EQUAL( m2_indices_d.size(), m2_h_d.size() );
00599     unsigned h = 0;
00600     for( unsigned r = 0; r < m1_indices_h.size(); ++r )
00601     {
00602         CPPUNIT_ASSERT_MATRICES_EQUAL( m1_h_h[h], m1_h_d[r], 1e-6 );
00603         CPPUNIT_ASSERT_MATRICES_EQUAL( m2_h_h[h], m2_h_d[r], 1e-6 );
00604         h += ( m1_indices_h.size() - r );
00605     }
00606 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines