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