MOAB: Mesh Oriented datABase  (version 5.4.1)
QualityMetricTest.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     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     (2006) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 /*! \file QualityMetricTest.cpp
00028 
00029 Unit testing for the QualityMetric class
00030 \author Jason Kraftcheck
00031 */
00032 #include "Mesquite.hpp"
00033 #include "VertexQM.hpp"
00034 #include "ElementQM.hpp"
00035 #include "IdealElements.hpp"
00036 #include "UnitUtil.hpp"
00037 #include "TopologyInfo.hpp"
00038 #include "PatchData.hpp"
00039 #include "cppunit/extensions/HelperMacros.h"
00040 #include "TMPDerivs.hpp"
00041 
00042 #include <string>
00043 
00044 using namespace MBMesquite;
00045 
00046 class QualityMetricTest : public CppUnit::TestFixture
00047 {
00048     CPPUNIT_TEST_SUITE( QualityMetricTest );
00049     CPPUNIT_TEST( test_fixed_vertex_list );
00050     CPPUNIT_TEST( test_remove_fixed_gradients );
00051     CPPUNIT_TEST( test_remove_fixed_hessians );
00052     CPPUNIT_TEST( test_gradient_constant );
00053     CPPUNIT_TEST( test_gradient_linear );
00054     CPPUNIT_TEST( test_gradient_parabolic );
00055     CPPUNIT_TEST( test_gradient_tau );
00056     CPPUNIT_TEST( test_Hessian_constant );
00057     CPPUNIT_TEST( test_Hessian_linear );
00058     CPPUNIT_TEST( test_Hessian_parabolic );
00059     CPPUNIT_TEST( test_Hessian_tau );
00060     CPPUNIT_TEST( test_diagonal_constant );
00061     CPPUNIT_TEST( test_diagonal_linear );
00062     CPPUNIT_TEST( test_diagonal_parabolic );
00063     CPPUNIT_TEST( test_diagonal_tau );
00064     CPPUNIT_TEST_SUITE_END();
00065     PatchData tri_pd;
00066 
00067   public:
00068     void setUp();
00069 
00070     void test_fixed_vertex_list();
00071     void test_remove_fixed_gradients();
00072     void test_remove_fixed_hessians();
00073     void test_gradient_constant();
00074     void test_gradient_linear();
00075     void test_gradient_parabolic();
00076     void test_gradient_tau();
00077     void test_Hessian_constant();
00078     void test_Hessian_linear();
00079     void test_Hessian_parabolic();
00080     void test_Hessian_tau();
00081     void test_diagonal_constant();
00082     void test_diagonal_linear();
00083     void test_diagonal_parabolic();
00084     void test_diagonal_tau();
00085 
00086     void compare_indices( QualityMetric& qm,
00087                           PatchData& pd,
00088                           size_t sample,
00089                           double value,
00090                           const std::vector< size_t >& indices );
00091     void compare_gradient( QualityMetric& qm,
00092                            PatchData& pd,
00093                            size_t sample,
00094                            double value,
00095                            const std::vector< size_t >& indices,
00096                            const std::vector< Vector3D >& grad );
00097 };
00098 
00099 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( QualityMetricTest, "QualityMetricTest" );
00100 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( QualityMetricTest, "Unit" );
00101 
00102 // Create a single-triangle patch with one fixed vertex
00103 void QualityMetricTest::setUp()
00104 {
00105     MsqError err;
00106     const double vtx_coords[]   = { 1, 0, 0, 0, 1, 0, 3, 0, 0, 0, 0, 0, 2, 1, 0 };
00107     const size_t connectivity[] = { 0, 1, 3, 1, 2, 4 };
00108     const bool fixed[]          = { false, false, false, true, true };
00109     tri_pd.fill( 5, vtx_coords, 2, TRIANGLE, connectivity, fixed, err );
00110     CPPUNIT_ASSERT( !err );
00111 }
00112 
00113 // tolerance on numerical gradient/hessian values
00114 const double EPSILON = 5e-3;
00115 
00116 /**\brief Fake quality metric for testing numerial gradient
00117  *
00118  * Implement a vertex metric for which the "quality" at a given
00119  * vertex is the value of the x-coordinate of that vertex.  Thus
00120  * the gradient of the "quality" at every vertex should be {1,0,0}.
00121  */
00122 class LinearVertexMetric : public VertexQM
00123 {
00124   public:
00125     std::string get_name() const
00126     {
00127         return "Fake metric for testing numerical gradient";
00128     }
00129     int get_negate_flag() const
00130     {
00131         return 1;
00132     }
00133     bool evaluate( PatchData& pd, size_t vtx_idx, double& value, MsqError& )
00134     {
00135         value = pd.vertex_by_index( vtx_idx )[0];
00136         return true;
00137     }
00138     bool evaluate_with_indices( PatchData& pd,
00139                                 size_t vtx_idx,
00140                                 double& value,
00141                                 std::vector< size_t >& indices,
00142                                 MsqError& )
00143     {
00144         value = pd.vertex_by_index( vtx_idx )[0];
00145         indices.resize( 1 );
00146         indices[0] = vtx_idx;
00147         return true;
00148     }
00149 };
00150 
00151 /**\brief Fake quality metric for testing numerial gradient
00152  *
00153  * Implement an element metric for which the "quality" is always 1.0.
00154  * The resulting gradient and Hessian should always be zero.
00155  */
00156 class ConstantElementMetric : public ElementQM
00157 {
00158   public:
00159     std::string get_name() const
00160     {
00161         return "Fake metric for testing numerical gradient";
00162     }
00163     int get_negate_flag() const
00164     {
00165         return 1;
00166     }
00167     bool evaluate( PatchData& pd, size_t, double& value, MsqError& )
00168     {
00169         value = 1.0;
00170         return true;
00171     }
00172     bool evaluate_with_indices( PatchData& pd,
00173                                 size_t elem_idx,
00174                                 double& value,
00175                                 std::vector< size_t >& indices,
00176                                 MsqError& )
00177     {
00178         MsqMeshEntity& elem = pd.element_by_index( elem_idx );
00179         unsigned nv         = elem.node_count();
00180         const size_t* conn  = elem.get_vertex_index_array();
00181         indices.clear();
00182         for( unsigned i = 0; i < nv; ++i )
00183             if( conn[i] < pd.num_free_vertices() ) indices.push_back( conn[i] );
00184 
00185         value = 1.0;
00186         return true;
00187     }
00188 };
00189 
00190 /**\brief Fake quality metric for testing numerial gradient
00191  *
00192  * Implement a vertex metric for which the "quality" at a given
00193  * vertex is the square of the value of the y-coordinate of that vertex.  Thus
00194  * the Hessian of the "quality" at every vertex should be {2,0,0}.
00195  */
00196 class ParabolicVertexMetric : public VertexQM
00197 {
00198   public:
00199     std::string get_name() const
00200     {
00201         return "Fake metric for testing numerical gradient";
00202     }
00203     int get_negate_flag() const
00204     {
00205         return 1;
00206     }
00207     bool evaluate( PatchData& pd, size_t vtx_idx, double& value, MsqError& )
00208     {
00209         value = pd.vertex_by_index( vtx_idx )[1];
00210         value *= value;
00211         return true;
00212     }
00213     bool evaluate_with_indices( PatchData& pd,
00214                                 size_t vtx_idx,
00215                                 double& value,
00216                                 std::vector< size_t >& indices,
00217                                 MsqError& )
00218     {
00219         value = pd.vertex_by_index( vtx_idx )[1];
00220         value *= value;
00221         indices.resize( 1 );
00222         indices[0] = vtx_idx;
00223         return true;
00224     }
00225 };
00226 
00227 /**\brief More complex fake metric for testing finite difference
00228  *
00229  * Calculate /f$ (det(M)-1)^2 /f$ where the rows of M are the
00230  * coordinates of the three vertices in the triangle.
00231  */
00232 class TriTauMetric : public ElementQM
00233 {
00234   public:
00235     std::string get_name() const
00236     {
00237         return "Triangle Tau Metric";
00238     }
00239     int get_negate_flag() const
00240     {
00241         return 1;
00242     }
00243     MsqMatrix< 3, 3 > matrix( PatchData& pd, size_t elem_idx )
00244     {
00245         MsqMeshEntity& elem = pd.element_by_index( elem_idx );
00246         unsigned nv         = elem.node_count();
00247         const size_t* conn  = elem.get_vertex_index_array();
00248         CPPUNIT_ASSERT( nv == 3 );
00249 
00250         MsqMatrix< 3, 3 > M;
00251         M.set_row( 0, pd.vertex_by_index( conn[0] ).to_array() );
00252         M.set_row( 1, pd.vertex_by_index( conn[1] ).to_array() );
00253         M.set_row( 2, pd.vertex_by_index( conn[2] ).to_array() );
00254 
00255         return M;
00256     }
00257 
00258     bool evaluate( PatchData& pd, size_t element_idx, double& value, MsqError& )
00259     {
00260         MsqMatrix< 3, 3 > M = matrix( pd, element_idx );
00261         value               = det( M ) - 1;
00262         value *= value;
00263         return true;
00264     }
00265 
00266     bool evaluate_with_indices( PatchData& pd,
00267                                 size_t elem_idx,
00268                                 double& value,
00269                                 std::vector< size_t >& indices,
00270                                 MsqError& )
00271     {
00272         MsqMatrix< 3, 3 > M = matrix( pd, elem_idx );
00273         value               = det( M ) - 1;
00274         value *= value;
00275 
00276         indices.clear();
00277         unsigned nv        = pd.element_by_index( elem_idx ).node_count();
00278         const size_t* conn = pd.element_by_index( elem_idx ).get_vertex_index_array();
00279         for( unsigned i = 0; i < nv; ++i )
00280             if( pd.is_vertex_free( conn[i] ) ) indices.push_back( conn[i] );
00281 
00282         return true;
00283     }
00284 };
00285 
00286 void QualityMetricTest::test_fixed_vertex_list()
00287 {
00288     // define a patch of four quads such that
00289     // the number of fixed vertices in each quad
00290     // varies from 0 to 3
00291     /*   2------1-----8
00292      *   |      |     |\
00293      *   | (0)  | (3) | \
00294      *   |      |     |  \
00295      *   3------0-----7   \
00296      *   |      |     |\   \
00297      *   | (1)  | (2) | \   \
00298      *   |      |     |  \   \
00299      *   4------5-----6   \   \
00300      *           \____\____\___\__ fixed
00301      */
00302     const double coords[] = { 0, 0, 0, 0, 1, 0, -1, 1, 0, -1, 0, 0, -1, -1, 0, 0, -1, 0, 1, -1, 0, 1, 0, 0, 1, 1, 0 };
00303     const size_t conn[]   = { 0, 1, 2, 3, 3, 4, 5, 0, 6, 7, 0, 5, 0, 7, 8, 1 };
00304     const bool fixed[]    = { false, false, false, false, false, true, true, true, true };
00305 
00306     MsqPrintError err( std::cout );
00307     PatchData pd;
00308     pd.fill( 9, coords, 4, QUADRILATERAL, conn, fixed, err );
00309     ASSERT_NO_ERROR( err );
00310 
00311     uint32_t bits;
00312     std::vector< size_t > indices;
00313     std::vector< size_t >::iterator it;
00314     const size_t* verts;
00315     unsigned i;
00316 
00317     indices.clear();
00318     bits = QualityMetric::fixed_vertex_bitmap( pd, &pd.element_by_index( 0 ), indices );
00319     CPPUNIT_ASSERT_EQUAL( (size_t)4, indices.size() );
00320     CPPUNIT_ASSERT_EQUAL( (uint32_t)0, bits & 0xF );
00321     CPPUNIT_ASSERT( pd.num_free_vertices() > *std::max_element( indices.begin(), indices.end() ) );
00322     verts = pd.element_by_index( 0 ).get_vertex_index_array();
00323     for( i = 0; i < 4; ++i )
00324         CPPUNIT_ASSERT( std::find( verts, verts + 4, indices[i] ) != verts + 4 );
00325 
00326     indices.clear();
00327     bits = QualityMetric::fixed_vertex_bitmap( pd, &pd.element_by_index( 1 ), indices );
00328     CPPUNIT_ASSERT_EQUAL( (size_t)3, indices.size() );
00329     verts = pd.element_by_index( 1 ).get_vertex_index_array();
00330     for( i = 0; i < 4; ++i )
00331     {
00332         it = std::find( indices.begin(), indices.end(), verts[i] );
00333         if( verts[i] < pd.num_free_vertices() )
00334         {
00335             CPPUNIT_ASSERT( it != indices.end() );
00336             CPPUNIT_ASSERT_EQUAL( 0u, bits & ( 1 << i ) );
00337         }
00338         else
00339         {
00340             CPPUNIT_ASSERT( it == indices.end() );
00341             CPPUNIT_ASSERT_EQUAL( 1u, ( bits >> i ) & 1 );
00342         }
00343     }
00344 
00345     indices.clear();
00346     bits = QualityMetric::fixed_vertex_bitmap( pd, &pd.element_by_index( 2 ), indices );
00347     CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );
00348     verts = pd.element_by_index( 2 ).get_vertex_index_array();
00349     for( i = 0; i < 4; ++i )
00350     {
00351         it = std::find( indices.begin(), indices.end(), verts[i] );
00352         if( verts[i] < pd.num_free_vertices() )
00353         {
00354             CPPUNIT_ASSERT( it != indices.end() );
00355             CPPUNIT_ASSERT_EQUAL( 0u, bits & ( 1 << i ) );
00356         }
00357         else
00358         {
00359             CPPUNIT_ASSERT( it == indices.end() );
00360             CPPUNIT_ASSERT_EQUAL( 1u, ( bits >> i ) & 1 );
00361         }
00362     }
00363 
00364     indices.clear();
00365     bits = QualityMetric::fixed_vertex_bitmap( pd, &pd.element_by_index( 3 ), indices );
00366     CPPUNIT_ASSERT_EQUAL( (size_t)2, indices.size() );
00367     verts = pd.element_by_index( 3 ).get_vertex_index_array();
00368     for( i = 0; i < 4; ++i )
00369     {
00370         it = std::find( indices.begin(), indices.end(), verts[i] );
00371         if( verts[i] < pd.num_free_vertices() )
00372         {
00373             CPPUNIT_ASSERT( it != indices.end() );
00374             CPPUNIT_ASSERT_EQUAL( 0u, bits & ( 1 << i ) );
00375         }
00376         else
00377         {
00378             CPPUNIT_ASSERT( it == indices.end() );
00379             CPPUNIT_ASSERT_EQUAL( 1u, ( bits >> i ) & 1 );
00380         }
00381     }
00382 }
00383 
00384 void QualityMetricTest::test_remove_fixed_gradients()
00385 {
00386     // define a list of vectors
00387     std::vector< Vector3D > grads( 6 );
00388     grads[0] = Vector3D( 0, 0, 0 );
00389     grads[1] = Vector3D( 1, 1, 1 );
00390     grads[2] = Vector3D( 2, 2, 2 );
00391     grads[3] = Vector3D( 3, 3, 3 );
00392     grads[4] = Vector3D( 4, 4, 4 );
00393     grads[5] = Vector3D( 5, 5, 5 );
00394     // remove the first, third, and last
00395     uint32_t flags = 1u | 4u | 32u;
00396     // call function, choose element type w/ correct number of corners
00397     QualityMetric::remove_fixed_gradients( PRISM, flags, grads );
00398     // check results
00399     CPPUNIT_ASSERT_EQUAL( (size_t)3, grads.size() );
00400     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 1, 1, 1 ), grads[0], DBL_EPSILON );
00401     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 3, 3, 3 ), grads[1], DBL_EPSILON );
00402     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 4, 4, 4 ), grads[2], DBL_EPSILON );
00403 }
00404 
00405 void QualityMetricTest::test_remove_fixed_hessians()
00406 {
00407     // define Hessian matrix for a quadrilateral
00408     Matrix3D m[10] = { Matrix3D( 0.0 ), Matrix3D( 1.0 ), Matrix3D( 2.0 ), Matrix3D( 3.0 ), Matrix3D( 4.0 ),
00409                        Matrix3D( 5.0 ), Matrix3D( 6.0 ), Matrix3D( 7.0 ), Matrix3D( 8.0 ), Matrix3D( 9.0 ) };
00410     // convert to std::vector
00411     std::vector< Matrix3D > hess( 10 );
00412     std::copy( m, m + 10, hess.begin() );
00413     // mark fist and third vertices as fixed
00414     uint32_t flags = 1u | 4u;
00415     // call function to remove grads for fixed vertices
00416     QualityMetric::remove_fixed_hessians( QUADRILATERAL, flags, hess );
00417     // the submatrix with the first and third rows/columns removed should be
00418     // { 4, 6,
00419     //      9 }
00420     CPPUNIT_ASSERT_EQUAL( (size_t)3, hess.size() );
00421     CPPUNIT_ASSERT_MATRICES_EQUAL( Matrix3D( 4.0 ), hess[0], DBL_EPSILON );
00422     CPPUNIT_ASSERT_MATRICES_EQUAL( Matrix3D( 6.0 ), hess[1], DBL_EPSILON );
00423     CPPUNIT_ASSERT_MATRICES_EQUAL( Matrix3D( 9.0 ), hess[2], DBL_EPSILON );
00424 }
00425 
00426 void QualityMetricTest::test_gradient_constant()
00427 {
00428     MsqError err;
00429     std::vector< size_t > indices;
00430     std::vector< Vector3D > gradient;
00431     double value;
00432     size_t ELEMENT = 0;
00433 
00434     ConstantElementMetric constant;
00435     QualityMetric& qm = constant;
00436 
00437     bool valid = qm.evaluate_with_gradient( tri_pd, ELEMENT, value, indices, gradient, err );
00438     CPPUNIT_ASSERT( !err );
00439     CPPUNIT_ASSERT( valid );
00440     compare_indices( qm, tri_pd, ELEMENT, value, indices );
00441     CPPUNIT_ASSERT_EQUAL( indices.size(), gradient.size() );
00442     CPPUNIT_ASSERT_EQUAL( (size_t)2, indices.size() );  // two free vertices in triangle
00443 
00444     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, 0, 0 ), gradient[0], EPSILON );
00445     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, 0, 0 ), gradient[1], EPSILON );
00446 }
00447 
00448 void QualityMetricTest::test_gradient_linear()
00449 {
00450     MsqError err;
00451     std::vector< size_t > indices;
00452     std::vector< Vector3D > gradient;
00453     double value;
00454     size_t VERTEX = 0;
00455 
00456     LinearVertexMetric linear;
00457     QualityMetric& qm = linear;
00458 
00459     bool valid = qm.evaluate_with_gradient( tri_pd, VERTEX, value, indices, gradient, err );
00460     CPPUNIT_ASSERT( !err );
00461     CPPUNIT_ASSERT( valid );
00462     compare_indices( qm, tri_pd, VERTEX, value, indices );
00463     CPPUNIT_ASSERT_EQUAL( indices.size(), gradient.size() );
00464     CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );
00465     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 1, 0, 0 ), gradient[0], EPSILON );
00466 }
00467 
00468 void QualityMetricTest::test_gradient_parabolic()
00469 {
00470     const size_t VERTEX = 1;  // pick vertex with non-zero Y-coordinate
00471     MsqError err;
00472     std::vector< size_t > indices;
00473     std::vector< Vector3D > gradient;
00474     double value;
00475 
00476     ParabolicVertexMetric parab;
00477     QualityMetric& qm = parab;
00478 
00479     bool valid = qm.evaluate_with_gradient( tri_pd, VERTEX, value, indices, gradient, err );
00480     CPPUNIT_ASSERT( !err );
00481     CPPUNIT_ASSERT( valid );
00482     compare_indices( qm, tri_pd, VERTEX, value, indices );
00483     CPPUNIT_ASSERT_EQUAL( indices.size(), gradient.size() );
00484     CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );  // two free vertices in triangle
00485 
00486     const double expected_y = 2 * tri_pd.vertex_by_index( VERTEX )[1];
00487     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, expected_y, 0 ), gradient[0], EPSILON );
00488 }
00489 
00490 void QualityMetricTest::test_gradient_tau()
00491 {
00492     MsqError err;
00493     std::vector< size_t > indices;
00494     std::vector< Vector3D > gradient;
00495     double value;
00496     size_t ELEMENT = 1;
00497 
00498     TriTauMetric qm;
00499 
00500     bool valid = qm.evaluate_with_gradient( tri_pd, ELEMENT, value, indices, gradient, err );
00501     CPPUNIT_ASSERT( !err );
00502     CPPUNIT_ASSERT( valid );
00503     compare_indices( qm, tri_pd, ELEMENT, value, indices );
00504     CPPUNIT_ASSERT_EQUAL( indices.size(), gradient.size() );
00505     CPPUNIT_ASSERT_EQUAL( (size_t)2, indices.size() );  // two free vertices in triangle
00506 
00507     MsqMatrix< 3, 3 > M  = qm.matrix( tri_pd, ELEMENT );
00508     MsqMatrix< 3, 3 > dM = transpose( adj( M ) );
00509     dM *= 2 * ( det( M ) - 1 );
00510 
00511     if( indices[0] == 2 )
00512     {
00513         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices[1] );
00514         std::swap( indices[0], indices[1] );
00515         std::swap( gradient[0], gradient[1] );
00516     }
00517     else
00518     {
00519         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices[0] );
00520         CPPUNIT_ASSERT_EQUAL( (size_t)2, indices[1] );
00521     }
00522     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( dM.row( 0 ).data() ), gradient[0], EPSILON );
00523     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( dM.row( 1 ).data() ), gradient[1], EPSILON );
00524 }
00525 
00526 void QualityMetricTest::test_Hessian_constant()
00527 {
00528     MsqError err;
00529     std::vector< size_t > indices;
00530     std::vector< Vector3D > gradient;
00531     std::vector< Matrix3D > Hessian;
00532     double value;
00533     size_t ELEMENT = 0;
00534 
00535     ConstantElementMetric constant;
00536     QualityMetric& qm = constant;
00537 
00538     bool valid = qm.evaluate_with_Hessian( tri_pd, ELEMENT, value, indices, gradient, Hessian, err );
00539     CPPUNIT_ASSERT( !err );
00540     CPPUNIT_ASSERT( valid );
00541     compare_gradient( qm, tri_pd, ELEMENT, value, indices, gradient );
00542     CPPUNIT_ASSERT_EQUAL( indices.size() * ( indices.size() + 1 ) / 2, Hessian.size() );
00543     CPPUNIT_ASSERT_EQUAL( (size_t)2, indices.size() );
00544 
00545     CPPUNIT_ASSERT_MATRICES_EQUAL( Matrix3D( 0.0 ), Hessian[0], EPSILON );
00546     CPPUNIT_ASSERT_MATRICES_EQUAL( Matrix3D( 0.0 ), Hessian[1], EPSILON );
00547 }
00548 
00549 void QualityMetricTest::test_Hessian_linear()
00550 {
00551     MsqError err;
00552     std::vector< size_t > indices;
00553     std::vector< Vector3D > gradient;
00554     std::vector< Matrix3D > Hessian;
00555     double value;
00556     size_t VERTEX = 0;
00557 
00558     LinearVertexMetric linear;
00559     QualityMetric& qm = linear;
00560 
00561     bool valid = qm.evaluate_with_Hessian( tri_pd, VERTEX, value, indices, gradient, Hessian, err );
00562     CPPUNIT_ASSERT( !err );
00563     CPPUNIT_ASSERT( valid );
00564     compare_gradient( qm, tri_pd, VERTEX, value, indices, gradient );
00565     CPPUNIT_ASSERT_EQUAL( indices.size() * ( indices.size() + 1 ) / 2, Hessian.size() );
00566     CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );
00567     CPPUNIT_ASSERT_MATRICES_EQUAL( Matrix3D( 0.0 ), Hessian[0], EPSILON );
00568 }
00569 
00570 void QualityMetricTest::test_Hessian_parabolic()
00571 {
00572     MsqError err;
00573     std::vector< size_t > indices;
00574     std::vector< Vector3D > gradient;
00575     std::vector< Matrix3D > Hessian;
00576     double value;
00577     size_t VERTEX = 0;
00578 
00579     ParabolicVertexMetric parab;
00580     QualityMetric& qm = parab;
00581 
00582     bool valid = qm.evaluate_with_Hessian( tri_pd, VERTEX, value, indices, gradient, Hessian, err );
00583     CPPUNIT_ASSERT( !err );
00584     CPPUNIT_ASSERT( valid );
00585     compare_gradient( qm, tri_pd, VERTEX, value, indices, gradient );
00586     CPPUNIT_ASSERT_EQUAL( indices.size() * ( indices.size() + 1 ) / 2, Hessian.size() );
00587     CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );
00588 
00589     Matrix3D expected( 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0 );
00590     CPPUNIT_ASSERT_MATRICES_EQUAL( expected, Hessian[0], EPSILON );
00591 }
00592 
00593 void QualityMetricTest::test_Hessian_tau()
00594 {
00595     MsqError err;
00596     std::vector< size_t > indices;
00597     std::vector< Vector3D > gradient;
00598     std::vector< Matrix3D > Hessian;
00599     double value;
00600     size_t ELEMENT = 1;
00601 
00602     TriTauMetric qm;
00603 
00604     bool valid = qm.evaluate_with_Hessian( tri_pd, ELEMENT, value, indices, gradient, Hessian, err );
00605     CPPUNIT_ASSERT( !err );
00606     CPPUNIT_ASSERT( valid );
00607     compare_gradient( qm, tri_pd, ELEMENT, value, indices, gradient );
00608     CPPUNIT_ASSERT_EQUAL( indices.size() * ( indices.size() + 1 ) / 2, Hessian.size() );
00609     CPPUNIT_ASSERT_EQUAL( (size_t)2, indices.size() );  // two free vertices in triangle
00610 
00611     MsqMatrix< 3, 3 > M = qm.matrix( tri_pd, ELEMENT );
00612     MsqMatrix< 3, 3 > d2M[6];
00613     set_scaled_2nd_deriv_of_det( d2M, 2 * ( det( M ) - 1 ), M );
00614     pluseq_scaled_outer_product( d2M, 2, transpose( adj( M ) ) );
00615 
00616     if( indices[0] == 2 )
00617     {
00618         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices[1] );
00619         std::swap( indices[0], indices[1] );
00620         std::swap( gradient[0], gradient[1] );
00621         std::swap( Hessian[0], Hessian[2] );
00622         Hessian[1] = transpose( Hessian[1] );
00623     }
00624     else
00625     {
00626         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices[0] );
00627         CPPUNIT_ASSERT_EQUAL( (size_t)2, indices[1] );
00628     }
00629 
00630     CPPUNIT_ASSERT_MATRICES_EQUAL( Matrix3D( d2M[0].data() ), Hessian[0], EPSILON );
00631     CPPUNIT_ASSERT_MATRICES_EQUAL( Matrix3D( d2M[1].data() ), Hessian[1], EPSILON );
00632     CPPUNIT_ASSERT_MATRICES_EQUAL( Matrix3D( d2M[3].data() ), Hessian[2], EPSILON );
00633 }
00634 
00635 void QualityMetricTest::test_diagonal_constant()
00636 {
00637     MsqError err;
00638     std::vector< size_t > indices;
00639     std::vector< Vector3D > gradient;
00640     std::vector< SymMatrix3D > Hessian;
00641     double value;
00642     size_t ELEMENT = 0;
00643 
00644     ConstantElementMetric constant;
00645     QualityMetric& qm = constant;
00646 
00647     bool valid = qm.evaluate_with_Hessian_diagonal( tri_pd, ELEMENT, value, indices, gradient, Hessian, err );
00648     CPPUNIT_ASSERT( !err );
00649     CPPUNIT_ASSERT( valid );
00650     compare_gradient( qm, tri_pd, ELEMENT, value, indices, gradient );
00651     CPPUNIT_ASSERT_EQUAL( indices.size(), Hessian.size() );
00652     CPPUNIT_ASSERT_EQUAL( (size_t)2, indices.size() );
00653 
00654     CPPUNIT_ASSERT_MATRICES_EQUAL( SymMatrix3D( 0.0 ), Hessian[0], EPSILON );
00655     CPPUNIT_ASSERT_MATRICES_EQUAL( SymMatrix3D( 0.0 ), Hessian[1], EPSILON );
00656 }
00657 
00658 void QualityMetricTest::test_diagonal_linear()
00659 {
00660     MsqError err;
00661     std::vector< size_t > indices;
00662     std::vector< Vector3D > gradient;
00663     std::vector< SymMatrix3D > Hessian;
00664     double value;
00665     size_t VERTEX = 0;
00666 
00667     LinearVertexMetric linear;
00668     QualityMetric& qm = linear;
00669 
00670     bool valid = qm.evaluate_with_Hessian_diagonal( tri_pd, VERTEX, value, indices, gradient, Hessian, err );
00671     CPPUNIT_ASSERT( !err );
00672     CPPUNIT_ASSERT( valid );
00673     compare_gradient( qm, tri_pd, VERTEX, value, indices, gradient );
00674     CPPUNIT_ASSERT_EQUAL( indices.size(), Hessian.size() );
00675     CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );
00676     CPPUNIT_ASSERT_MATRICES_EQUAL( SymMatrix3D( 0.0 ), Hessian[0], EPSILON );
00677 }
00678 
00679 void QualityMetricTest::test_diagonal_parabolic()
00680 {
00681     MsqError err;
00682     std::vector< size_t > indices;
00683     std::vector< Vector3D > gradient;
00684     std::vector< SymMatrix3D > Hessian;
00685     double value;
00686     size_t VERTEX = 0;
00687 
00688     ParabolicVertexMetric parab;
00689     QualityMetric& qm = parab;
00690 
00691     bool valid = qm.evaluate_with_Hessian_diagonal( tri_pd, VERTEX, value, indices, gradient, Hessian, err );
00692     CPPUNIT_ASSERT( !err );
00693     CPPUNIT_ASSERT( valid );
00694     compare_gradient( qm, tri_pd, VERTEX, value, indices, gradient );
00695     CPPUNIT_ASSERT_EQUAL( indices.size(), Hessian.size() );
00696     CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );
00697 
00698     SymMatrix3D expected( 0.0, 0.0, 0.0, 2.0, 0.0, 0.0 );
00699     CPPUNIT_ASSERT_MATRICES_EQUAL( expected, Hessian[0], EPSILON );
00700 }
00701 
00702 void QualityMetricTest::test_diagonal_tau()
00703 {
00704     MsqError err;
00705     std::vector< size_t > indices;
00706     std::vector< Vector3D > gradient;
00707     std::vector< SymMatrix3D > Hessian;
00708     double value;
00709     size_t ELEMENT = 1;
00710 
00711     TriTauMetric qm;
00712 
00713     bool valid = qm.evaluate_with_Hessian_diagonal( tri_pd, ELEMENT, value, indices, gradient, Hessian, err );
00714     CPPUNIT_ASSERT( !err );
00715     CPPUNIT_ASSERT( valid );
00716     compare_gradient( qm, tri_pd, ELEMENT, value, indices, gradient );
00717     CPPUNIT_ASSERT_EQUAL( indices.size(), Hessian.size() );
00718     CPPUNIT_ASSERT_EQUAL( (size_t)2, indices.size() );  // two free vertices in triangle
00719 
00720     MsqMatrix< 3, 3 > M = qm.matrix( tri_pd, ELEMENT );
00721     MsqMatrix< 3, 3 > d2M[6];
00722     set_scaled_2nd_deriv_of_det( d2M, 2 * ( det( M ) - 1 ), M );
00723     pluseq_scaled_outer_product( d2M, 2, transpose( adj( M ) ) );
00724 
00725     if( indices[0] == 2 )
00726     {
00727         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices[1] );
00728         std::swap( indices[0], indices[1] );
00729         std::swap( gradient[0], gradient[1] );
00730         std::swap( Hessian[0], Hessian[1] );
00731     }
00732     else
00733     {
00734         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices[0] );
00735         CPPUNIT_ASSERT_EQUAL( (size_t)2, indices[1] );
00736     }
00737 
00738     SymMatrix3D exp0( d2M[0]( 0, 0 ), d2M[0]( 0, 1 ), d2M[0]( 0, 2 ), d2M[0]( 1, 1 ), d2M[0]( 1, 2 ), d2M[0]( 2, 2 ) );
00739     SymMatrix3D exp1( d2M[3]( 0, 0 ), d2M[3]( 0, 1 ), d2M[3]( 0, 2 ), d2M[3]( 1, 1 ), d2M[3]( 1, 2 ), d2M[3]( 2, 2 ) );
00740 
00741     CPPUNIT_ASSERT_MATRICES_EQUAL( exp0, Hessian[0], EPSILON );
00742     CPPUNIT_ASSERT_MATRICES_EQUAL( exp1, Hessian[1], EPSILON );
00743 }
00744 
00745 void QualityMetricTest::compare_indices( QualityMetric& qm,
00746                                          PatchData& pd,
00747                                          size_t sample,
00748                                          double value,
00749                                          const std::vector< size_t >& indices )
00750 {
00751     double value2;
00752     std::vector< size_t > indices2;
00753     MsqError err;
00754     bool valid = qm.evaluate_with_indices( pd, sample, value2, indices2, err );
00755     ASSERT_NO_ERROR( err );
00756     CPPUNIT_ASSERT( valid );
00757     CPPUNIT_ASSERT_DOUBLES_EQUAL( value2, value, EPSILON );
00758     std::vector< size_t > indices1( indices );
00759     std::sort( indices1.begin(), indices1.end() );
00760     std::sort( indices2.begin(), indices2.end() );
00761     ASSERT_STD_VECTORS_EQUAL( indices2, indices1 );
00762 }
00763 
00764 void QualityMetricTest::compare_gradient( QualityMetric& qm,
00765                                           PatchData& pd,
00766                                           size_t sample,
00767                                           double value,
00768                                           const std::vector< size_t >& indices,
00769                                           const std::vector< Vector3D >& grad )
00770 {
00771     double value2;
00772     std::vector< size_t > indices2;
00773     std::vector< Vector3D > grad2;
00774     MsqError err;
00775     bool valid = qm.evaluate_with_gradient( pd, sample, value2, indices2, grad2, err );
00776     ASSERT_NO_ERROR( err );
00777     CPPUNIT_ASSERT( valid );
00778     CPPUNIT_ASSERT_DOUBLES_EQUAL( value2, value, EPSILON );
00779     CPPUNIT_ASSERT_EQUAL( indices2.size(), indices.size() );
00780     CPPUNIT_ASSERT_EQUAL( indices.size(), grad.size() );
00781 
00782     for( size_t i = 0; i < indices.size(); ++i )
00783     {
00784         size_t j = std::find( indices2.begin(), indices2.end(), indices[i] ) - indices2.begin();
00785         CPPUNIT_ASSERT( j < indices2.size() );  // not found->indices don't match
00786         CPPUNIT_ASSERT_VECTORS_EQUAL( grad2[j], grad[i], EPSILON );
00787     }
00788 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines