MOAB: Mesh Oriented datABase  (version 5.4.0)
MeshUtilTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2010 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     (2010) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file MeshUtilTest.cpp
00028  *  \brief Test MeshUtil class
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "ArrayMesh.hpp"
00033 #include "MeshUtil.hpp"
00034 #include "UnitUtil.hpp"
00035 #include "SimpleStats.hpp"
00036 
00037 using namespace MBMesquite;
00038 
00039 class MeshUtilTest : public CppUnit::TestFixture
00040 {
00041   private:
00042     CPPUNIT_TEST_SUITE( MeshUtilTest );
00043     CPPUNIT_TEST( test_edge_length_distribution_types );
00044     CPPUNIT_TEST( test_edge_length_distribution_unique );
00045     CPPUNIT_TEST( test_edge_length_distribution_empty );
00046     CPPUNIT_TEST( test_lambda_distribution );
00047     CPPUNIT_TEST( test_lambda_distribution_mixed );
00048     CPPUNIT_TEST( test_lambda_distribution_empty );
00049     CPPUNIT_TEST_SUITE_END();
00050 
00051   public:
00052     // test that edge_length_distribution function works
00053     // for all supported element types
00054     void test_edge_length_distribution_types();
00055 
00056     // test that edge_length_distribution function works
00057     // counts each implicit edge only once
00058     void test_edge_length_distribution_unique();
00059 
00060     // test that edge_length_distribution function
00061     // handles the case of an empty mesh
00062     void test_edge_length_distribution_empty();
00063 
00064     // test numerical results of lambda_distribution
00065     void test_lambda_distribution();
00066 
00067     // test that lambda_distribution function works
00068     // for all supported element types
00069     void test_lambda_distribution_mixed();
00070 
00071     // test that lambda_distribution function
00072     // handles the case of an empty mesh
00073     void test_lambda_distribution_empty();
00074 };
00075 
00076 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( MeshUtilTest, "MeshUtilTest" );
00077 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( MeshUtilTest, "Unit" );
00078 
00079 const double EPSILON = 1e-8;
00080 
00081 void MeshUtilTest::test_edge_length_distribution_types()
00082 {
00083     // create a mesh that contains one of each element type,
00084     // completely disconnected from any other element
00085     double coords[]            = { // triangle : 0 - 2
00086                         0, 0, 0, 2, 0, 0, 1, 1, 0,
00087                         // quad : 3 - 6
00088                         3, 0, 0, 3, 1, 0, 3, 2, 2, 3, 0, 1,
00089                         // tet : 7 - 10
00090                         0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1,
00091                         // pyramid : 11 - 15
00092                         0, 0, 0, 0, -2, 0, -2, -2, 0, -1, 0, 0, 0, 0, -3,
00093                         // wedge : 16 - 21
00094                         0, 0, -2, 2, 0, -2, 1, 1, -2, 0, 0, -1.5, 2, 0, -1.5, 1, 1, -1, 5,
00095                         // hex : 22 - 29
00096                         4, 0, 0, 4, 1, 0, 4, 2, 2, 4, 0, 1, 5, 0, 0, 5, 1, 0, 5, 2, 2, 5, 0, 1 };
00097     const int tri_offset       = 0;
00098     const int quad_offset      = tri_offset + 3;
00099     const int tet_offset       = quad_offset + 4;
00100     const int pyr_offset       = tet_offset + 4;
00101     const int wedge_offset     = pyr_offset + 5;
00102     const int hex_offset       = wedge_offset + 6;
00103     const int num_vtx          = hex_offset + 8;
00104     const int fixed[num_vtx]   = { 0 };
00105     EntityTopology types[]     = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON };
00106     const unsigned long conn[] = {
00107         0,  1,  2,                      // triangle
00108         3,  4,  5,  6,                  // quadrilateral
00109         7,  8,  9,  10,                 // tetrahedron
00110         11, 12, 13, 14, 15,             // pyramid
00111         16, 17, 18, 19, 20, 21,         // prism,
00112         22, 23, 24, 25, 26, 27, 28, 29  // hex
00113     };
00114     ArrayMesh mesh( 3, num_vtx, coords, fixed, 6, types, conn );
00115 
00116     // calculate expected value
00117     const int edges[][2] = { // triangle
00118                              { 0, 1 },
00119                              { 1, 2 },
00120                              { 2, 0 },
00121                              // quad
00122                              { 3, 4 },
00123                              { 4, 5 },
00124                              { 5, 6 },
00125                              { 6, 3 },
00126                              // tet
00127                              { 7, 8 },
00128                              { 8, 9 },
00129                              { 9, 7 },
00130                              { 7, 10 },
00131                              { 8, 10 },
00132                              { 9, 10 },
00133                              // pyramid
00134                              { 11, 12 },
00135                              { 12, 13 },
00136                              { 13, 14 },
00137                              { 14, 11 },
00138                              { 11, 15 },
00139                              { 12, 15 },
00140                              { 13, 15 },
00141                              { 14, 15 },
00142                              // prism
00143                              { 16, 17 },
00144                              { 17, 18 },
00145                              { 18, 16 },
00146                              { 16, 19 },
00147                              { 17, 20 },
00148                              { 18, 21 },
00149                              { 19, 20 },
00150                              { 20, 21 },
00151                              { 21, 19 },
00152                              // hex
00153                              { 22, 23 },
00154                              { 23, 24 },
00155                              { 24, 25 },
00156                              { 25, 22 },
00157                              { 22, 26 },
00158                              { 23, 27 },
00159                              { 24, 28 },
00160                              { 25, 29 },
00161                              { 26, 27 },
00162                              { 27, 28 },
00163                              { 28, 29 },
00164                              { 29, 26 } };
00165     const int num_edges  = sizeof( edges ) / sizeof( edges[0] );
00166 
00167     double exp_min = HUGE_VAL, exp_max = -HUGE_VAL;
00168     double exp_avg = 0, exp_rms = 0;
00169     for( int i = 0; i < num_edges; ++i )
00170     {
00171         int v1 = edges[i][0];
00172         int v2 = edges[i][1];
00173         Vector3D c0( coords + 3 * v1 );
00174         Vector3D c1( coords + 3 * v2 );
00175         Vector3D diff = c1 - c0;
00176         double len    = sqrt( diff % diff );
00177         if( len < exp_min ) exp_min = len;
00178         if( len > exp_max ) exp_max = len;
00179         exp_avg += len;
00180         exp_rms += diff % diff;
00181     }
00182 
00183     exp_avg /= num_edges;
00184     exp_rms /= num_edges;
00185     double exp_std_dev = sqrt( exp_rms - exp_avg * exp_avg );
00186     exp_rms            = sqrt( exp_rms );
00187 
00188     SimpleStats act_results;
00189     MsqError err;
00190     MeshUtil tool( &mesh );
00191     tool.edge_length_distribution( act_results, err );
00192     ASSERT_NO_ERROR( err );
00193     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_min, act_results.minimum(), EPSILON );
00194     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_max, act_results.maximum(), EPSILON );
00195     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_avg, act_results.average(), EPSILON );
00196     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_rms, act_results.rms(), EPSILON );
00197     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_std_dev, act_results.standard_deviation(), EPSILON );
00198 }
00199 
00200 void MeshUtilTest::test_edge_length_distribution_unique()
00201 {
00202     // define two-tet mesh where tets share a face
00203     double coords[]            = { 0, -5, 0, 0, 5, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1 };
00204     const unsigned long conn[] = { 4, 3, 2, 0, 2, 3, 4, 1 };
00205     int fixed[5]               = { 0 };
00206     ArrayMesh mesh( 3, 5, coords, fixed, 2, TETRAHEDRON, conn );
00207 
00208     const int edges[][2] = { { 0, 2 }, { 0, 3 }, { 0, 4 }, { 1, 2 }, { 1, 3 }, { 1, 4 }, { 2, 3 }, { 3, 4 }, { 4, 2 } };
00209     const int num_edges  = sizeof( edges ) / sizeof( edges[0] );
00210 
00211     double exp_min = HUGE_VAL, exp_max = -HUGE_VAL;
00212     double exp_avg = 0, exp_rms = 0;
00213     for( int i = 0; i < num_edges; ++i )
00214     {
00215         int v1 = edges[i][0];
00216         int v2 = edges[i][1];
00217         Vector3D c0( coords + 3 * v1 );
00218         Vector3D c1( coords + 3 * v2 );
00219         Vector3D diff = c1 - c0;
00220         double len    = sqrt( diff % diff );
00221         if( len < exp_min ) exp_min = len;
00222         if( len > exp_max ) exp_max = len;
00223         exp_avg += len;
00224         exp_rms += diff % diff;
00225     }
00226 
00227     exp_avg /= num_edges;
00228     exp_rms /= num_edges;
00229     double exp_std_dev = sqrt( exp_rms - exp_avg * exp_avg );
00230     exp_rms            = sqrt( exp_rms );
00231 
00232     SimpleStats act_results;
00233     MsqError err;
00234     MeshUtil tool( &mesh );
00235     tool.edge_length_distribution( act_results, err );
00236     ASSERT_NO_ERROR( err );
00237     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_min, act_results.minimum(), EPSILON );
00238     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_max, act_results.maximum(), EPSILON );
00239     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_avg, act_results.average(), EPSILON );
00240     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_rms, act_results.rms(), EPSILON );
00241     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_std_dev, act_results.standard_deviation(), EPSILON );
00242 }
00243 
00244 void MeshUtilTest::test_edge_length_distribution_empty()
00245 {
00246     ArrayMesh empty( 3, 0, 0, 0, 0, 0, 0 );
00247     MeshUtil tool( &empty );
00248     SimpleStats act_results;
00249     MsqError err;
00250     tool.edge_length_distribution( act_results, err );
00251     CPPUNIT_ASSERT_EQUAL( MsqError::INVALID_MESH, err.error_code() );
00252 }
00253 
00254 void MeshUtilTest::test_lambda_distribution()
00255 {
00256     // simple regular unit quad mesh one deformed corner
00257     double coords[]            = { 0, 2, 0, 1, 2, 0, 2, 2, 0, 0, 1, 0, 1, 1, 0, 2, 1, 0, 0, 0, 0, 1, 0, 0, 3, 0, 0 };
00258     const unsigned long conn[] = { 3, 4, 1, 0, 4, 5, 2, 1, 6, 7, 4, 3, 7, 8, 5, 4 };
00259     const int fixed[9]         = { 0 };
00260 
00261     // first mesh doesn't contain deformed element
00262     SimpleStats results;
00263     MsqError err;
00264     ArrayMesh mesh3( 3, 8, coords, fixed, 3, QUADRILATERAL, conn );
00265     MeshUtil tool3( &mesh3 );
00266     tool3.lambda_distribution( results, err );
00267 
00268     // lambda should be one everywhere
00269     ASSERT_NO_ERROR( err );
00270     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, results.minimum(), 1e-12 );
00271     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, results.maximum(), 1e-12 );
00272     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, results.average(), 1e-12 );
00273     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, results.rms(), 1e-12 );
00274     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, results.variance(), 1e-12 );
00275 
00276     // now try with deformed element
00277     ArrayMesh mesh4( 3, 9, coords, fixed, 4, QUADRILATERAL, conn );
00278     MeshUtil tool4( &mesh4 );
00279     results.clear();
00280     tool4.lambda_distribution( results, err );
00281 
00282     // expect 13 values of 1 and a few other values
00283     Vector3D v4( coords + 3 * 4 );
00284     Vector3D v5( coords + 3 * 5 );
00285     Vector3D v7( coords + 3 * 7 );
00286     Vector3D v8( coords + 3 * 8 );
00287     double val14   = sqrt( ( ( v4 - v5 ) * ( v8 - v5 ) ).length() );
00288     double val15   = sqrt( ( ( v5 - v8 ) * ( v7 - v8 ) ).length() );
00289     double val16   = sqrt( ( ( v4 - v7 ) * ( v8 - v7 ) ).length() );
00290     double exp_min = std::min( std::min( val14, val15 ), std::min( val16, 1.0 ) );
00291     double exp_max = std::max( std::max( val14, val15 ), val16 );
00292     double exp_avg = ( 13 + val14 + val15 + val16 ) / 16.0;
00293     double exp_rms = sqrt( ( 13 + val14 * val14 + val15 * val16 + val16 * val16 ) / 16.0 );
00294 
00295     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_min, results.minimum(), EPSILON );
00296     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_max, results.maximum(), EPSILON );
00297     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_avg, results.average(), EPSILON );
00298     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_rms, results.rms(), EPSILON );
00299 }
00300 
00301 void MeshUtilTest::test_lambda_distribution_mixed()
00302 {
00303     // create a mesh that contains one of each element type,
00304     // completely disconnected from any other element
00305     double coords[]            = { // triangle : 0 - 2
00306                         0, 0, 0, 2, 0, 0, 1, 1, 0,
00307                         // quad : 3 - 6
00308                         3, 0, 0, 3, 1, 0, 3, 2, 2, 3, 0, 1,
00309                         // tet : 7 - 10
00310                         0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1,
00311                         // pyramid : 11 - 15
00312                         0, 0, 0, 0, -2, 0, -2, -2, 0, -1, 0, 0, 0, 0, -3,
00313                         // wedge : 16 - 21
00314                         0, 0, -2, 2, 0, -2, 1, 1, -2, 0, 0, -1.5, 2, 0, -1.5, 1, 1, -1, 5,
00315                         // hex : 22 - 29
00316                         4, 0, 0, 4, 1, 0, 4, 2, 2, 4, 0, 1, 5, 0, 0, 5, 1, 0, 5, 2, 2, 5, 0, 1 };
00317     const int tri_offset       = 0;
00318     const int quad_offset      = tri_offset + 3;
00319     const int tet_offset       = quad_offset + 4;
00320     const int pyr_offset       = tet_offset + 4;
00321     const int wedge_offset     = pyr_offset + 5;
00322     const int hex_offset       = wedge_offset + 6;
00323     const int num_vtx          = hex_offset + 8;
00324     const int fixed[num_vtx]   = { 0 };
00325     EntityTopology types[]     = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON };
00326     const unsigned long conn[] = {
00327         0,  1,  2,                      // triangle
00328         3,  4,  5,  6,                  // quadrilateral
00329         7,  8,  9,  10,                 // tetrahedron
00330         11, 12, 13, 14, 15,             // pyramid
00331         16, 17, 18, 19, 20, 21,         // prism,
00332         22, 23, 24, 25, 26, 27, 28, 29  // hex
00333     };
00334     ArrayMesh mesh( 3, num_vtx, coords, fixed, 6, types, conn );
00335 
00336     // just test that this doesn't fail for any of the element types
00337     SimpleStats act_results;
00338     MsqError err;
00339     MeshUtil tool( &mesh );
00340     tool.edge_length_distribution( act_results, err );
00341     ASSERT_NO_ERROR( err );
00342 }
00343 
00344 void MeshUtilTest::test_lambda_distribution_empty()
00345 {
00346     ArrayMesh empty( 3, 0, 0, 0, 0, 0, 0 );
00347     MeshUtil tool( &empty );
00348     SimpleStats act_results;
00349     MsqError err;
00350     tool.lambda_distribution( act_results, err );
00351     CPPUNIT_ASSERT_EQUAL( MsqError::INVALID_MESH, err.error_code() );
00352 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines