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