MOAB: Mesh Oriented datABase  (version 5.2.1)
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     };
00098     const int tri_offset       = 0;
00099     const int quad_offset      = tri_offset + 3;
00100     const int tet_offset       = quad_offset + 4;
00101     const int pyr_offset       = tet_offset + 4;
00102     const int wedge_offset     = pyr_offset + 5;
00103     const int hex_offset       = wedge_offset + 6;
00104     const int num_vtx          = hex_offset + 8;
00105     const int fixed[num_vtx]   = { 0 };
00106     EntityTopology types[]     = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON };
00107     const unsigned long conn[] = {
00108         0,  1,  2,                      // triangle
00109         3,  4,  5,  6,                  // quadrilateral
00110         7,  8,  9,  10,                 // tetrahedron
00111         11, 12, 13, 14, 15,             // pyramid
00112         16, 17, 18, 19, 20, 21,         // prism,
00113         22, 23, 24, 25, 26, 27, 28, 29  // hex
00114     };
00115     ArrayMesh mesh( 3, num_vtx, coords, fixed, 6, types, conn );
00116 
00117     // calculate expected value
00118     const int edges[][2] = { // triangle
00119                              { 0, 1 },
00120                              { 1, 2 },
00121                              { 2, 0 },
00122                              // quad
00123                              { 3, 4 },
00124                              { 4, 5 },
00125                              { 5, 6 },
00126                              { 6, 3 },
00127                              // tet
00128                              { 7, 8 },
00129                              { 8, 9 },
00130                              { 9, 7 },
00131                              { 7, 10 },
00132                              { 8, 10 },
00133                              { 9, 10 },
00134                              // pyramid
00135                              { 11, 12 },
00136                              { 12, 13 },
00137                              { 13, 14 },
00138                              { 14, 11 },
00139                              { 11, 15 },
00140                              { 12, 15 },
00141                              { 13, 15 },
00142                              { 14, 15 },
00143                              // prism
00144                              { 16, 17 },
00145                              { 17, 18 },
00146                              { 18, 16 },
00147                              { 16, 19 },
00148                              { 17, 20 },
00149                              { 18, 21 },
00150                              { 19, 20 },
00151                              { 20, 21 },
00152                              { 21, 19 },
00153                              // hex
00154                              { 22, 23 },
00155                              { 23, 24 },
00156                              { 24, 25 },
00157                              { 25, 22 },
00158                              { 22, 26 },
00159                              { 23, 27 },
00160                              { 24, 28 },
00161                              { 25, 29 },
00162                              { 26, 27 },
00163                              { 27, 28 },
00164                              { 28, 29 },
00165                              { 29, 26 }
00166     };
00167     const int num_edges = sizeof( edges ) / sizeof( edges[0] );
00168 
00169     double exp_min = HUGE_VAL, exp_max = -HUGE_VAL;
00170     double exp_avg = 0, exp_rms = 0;
00171     for( int i = 0; i < num_edges; ++i )
00172     {
00173         int v1 = edges[i][0];
00174         int v2 = edges[i][1];
00175         Vector3D c0( coords + 3 * v1 );
00176         Vector3D c1( coords + 3 * v2 );
00177         Vector3D diff = c1 - c0;
00178         double len    = sqrt( diff % diff );
00179         if( len < exp_min ) exp_min = len;
00180         if( len > exp_max ) exp_max = len;
00181         exp_avg += len;
00182         exp_rms += diff % diff;
00183     }
00184 
00185     exp_avg /= num_edges;
00186     exp_rms /= num_edges;
00187     double exp_std_dev = sqrt( exp_rms - exp_avg * exp_avg );
00188     exp_rms            = sqrt( exp_rms );
00189 
00190     SimpleStats act_results;
00191     MsqError err;
00192     MeshUtil tool( &mesh );
00193     tool.edge_length_distribution( act_results, err );
00194     ASSERT_NO_ERROR( err );
00195     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_min, act_results.minimum(), EPSILON );
00196     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_max, act_results.maximum(), EPSILON );
00197     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_avg, act_results.average(), EPSILON );
00198     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_rms, act_results.rms(), EPSILON );
00199     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_std_dev, act_results.standard_deviation(), EPSILON );
00200 }
00201 
00202 void MeshUtilTest::test_edge_length_distribution_unique()
00203 {
00204     // define two-tet mesh where tets share a face
00205     double coords[]            = { 0, -5, 0, 0, 5, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1 };
00206     const unsigned long conn[] = { 4, 3, 2, 0, 2, 3, 4, 1 };
00207     int fixed[5]               = { 0 };
00208     ArrayMesh mesh( 3, 5, coords, fixed, 2, TETRAHEDRON, conn );
00209 
00210     const int edges[][2] = { { 0, 2 }, { 0, 3 }, { 0, 4 }, { 1, 2 }, { 1, 3 }, { 1, 4 }, { 2, 3 }, { 3, 4 }, { 4, 2 } };
00211     const int num_edges  = sizeof( edges ) / sizeof( edges[0] );
00212 
00213     double exp_min = HUGE_VAL, exp_max = -HUGE_VAL;
00214     double exp_avg = 0, exp_rms = 0;
00215     for( int i = 0; i < num_edges; ++i )
00216     {
00217         int v1 = edges[i][0];
00218         int v2 = edges[i][1];
00219         Vector3D c0( coords + 3 * v1 );
00220         Vector3D c1( coords + 3 * v2 );
00221         Vector3D diff = c1 - c0;
00222         double len    = sqrt( diff % diff );
00223         if( len < exp_min ) exp_min = len;
00224         if( len > exp_max ) exp_max = len;
00225         exp_avg += len;
00226         exp_rms += diff % diff;
00227     }
00228 
00229     exp_avg /= num_edges;
00230     exp_rms /= num_edges;
00231     double exp_std_dev = sqrt( exp_rms - exp_avg * exp_avg );
00232     exp_rms            = sqrt( exp_rms );
00233 
00234     SimpleStats act_results;
00235     MsqError err;
00236     MeshUtil tool( &mesh );
00237     tool.edge_length_distribution( act_results, err );
00238     ASSERT_NO_ERROR( err );
00239     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_min, act_results.minimum(), EPSILON );
00240     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_max, act_results.maximum(), EPSILON );
00241     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_avg, act_results.average(), EPSILON );
00242     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_rms, act_results.rms(), EPSILON );
00243     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_std_dev, act_results.standard_deviation(), EPSILON );
00244 }
00245 
00246 void MeshUtilTest::test_edge_length_distribution_empty()
00247 {
00248     ArrayMesh empty( 3, 0, 0, 0, 0, 0, 0 );
00249     MeshUtil tool( &empty );
00250     SimpleStats act_results;
00251     MsqError err;
00252     tool.edge_length_distribution( act_results, err );
00253     CPPUNIT_ASSERT_EQUAL( MsqError::INVALID_MESH, err.error_code() );
00254 }
00255 
00256 void MeshUtilTest::test_lambda_distribution()
00257 {
00258     // simple regular unit quad mesh one deformed corner
00259     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 };
00260     const unsigned long conn[] = { 3, 4, 1, 0, 4, 5, 2, 1, 6, 7, 4, 3, 7, 8, 5, 4 };
00261     const int fixed[9]         = { 0 };
00262 
00263     // first mesh doesn't contain deformed element
00264     SimpleStats results;
00265     MsqError err;
00266     ArrayMesh mesh3( 3, 8, coords, fixed, 3, QUADRILATERAL, conn );
00267     MeshUtil tool3( &mesh3 );
00268     tool3.lambda_distribution( results, err );
00269 
00270     // lambda should be one everywhere
00271     ASSERT_NO_ERROR( err );
00272     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, results.minimum(), 1e-12 );
00273     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, results.maximum(), 1e-12 );
00274     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, results.average(), 1e-12 );
00275     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, results.rms(), 1e-12 );
00276     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, results.variance(), 1e-12 );
00277 
00278     // now try with deformed element
00279     ArrayMesh mesh4( 3, 9, coords, fixed, 4, QUADRILATERAL, conn );
00280     MeshUtil tool4( &mesh4 );
00281     results.clear();
00282     tool4.lambda_distribution( results, err );
00283 
00284     // expect 13 values of 1 and a few other values
00285     Vector3D v4( coords + 3 * 4 );
00286     Vector3D v5( coords + 3 * 5 );
00287     Vector3D v7( coords + 3 * 7 );
00288     Vector3D v8( coords + 3 * 8 );
00289     double val14   = sqrt( ( ( v4 - v5 ) * ( v8 - v5 ) ).length() );
00290     double val15   = sqrt( ( ( v5 - v8 ) * ( v7 - v8 ) ).length() );
00291     double val16   = sqrt( ( ( v4 - v7 ) * ( v8 - v7 ) ).length() );
00292     double exp_min = std::min( std::min( val14, val15 ), std::min( val16, 1.0 ) );
00293     double exp_max = std::max( std::max( val14, val15 ), val16 );
00294     double exp_avg = ( 13 + val14 + val15 + val16 ) / 16.0;
00295     double exp_rms = sqrt( ( 13 + val14 * val14 + val15 * val16 + val16 * val16 ) / 16.0 );
00296 
00297     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_min, results.minimum(), EPSILON );
00298     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_max, results.maximum(), EPSILON );
00299     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_avg, results.average(), EPSILON );
00300     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_rms, results.rms(), EPSILON );
00301 }
00302 
00303 void MeshUtilTest::test_lambda_distribution_mixed()
00304 {
00305     // create a mesh that contains one of each element type,
00306     // completely disconnected from any other element
00307     double coords[] = { // triangle : 0 - 2
00308                         0, 0, 0, 2, 0, 0, 1, 1, 0,
00309                         // quad : 3 - 6
00310                         3, 0, 0, 3, 1, 0, 3, 2, 2, 3, 0, 1,
00311                         // tet : 7 - 10
00312                         0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1,
00313                         // pyramid : 11 - 15
00314                         0, 0, 0, 0, -2, 0, -2, -2, 0, -1, 0, 0, 0, 0, -3,
00315                         // wedge : 16 - 21
00316                         0, 0, -2, 2, 0, -2, 1, 1, -2, 0, 0, -1.5, 2, 0, -1.5, 1, 1, -1, 5,
00317                         // hex : 22 - 29
00318                         4, 0, 0, 4, 1, 0, 4, 2, 2, 4, 0, 1, 5, 0, 0, 5, 1, 0, 5, 2, 2, 5, 0, 1
00319     };
00320     const int tri_offset       = 0;
00321     const int quad_offset      = tri_offset + 3;
00322     const int tet_offset       = quad_offset + 4;
00323     const int pyr_offset       = tet_offset + 4;
00324     const int wedge_offset     = pyr_offset + 5;
00325     const int hex_offset       = wedge_offset + 6;
00326     const int num_vtx          = hex_offset + 8;
00327     const int fixed[num_vtx]   = { 0 };
00328     EntityTopology types[]     = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON };
00329     const unsigned long conn[] = {
00330         0,  1,  2,                      // triangle
00331         3,  4,  5,  6,                  // quadrilateral
00332         7,  8,  9,  10,                 // tetrahedron
00333         11, 12, 13, 14, 15,             // pyramid
00334         16, 17, 18, 19, 20, 21,         // prism,
00335         22, 23, 24, 25, 26, 27, 28, 29  // hex
00336     };
00337     ArrayMesh mesh( 3, num_vtx, coords, fixed, 6, types, conn );
00338 
00339     // just test that this doesn't fail for any of the element types
00340     SimpleStats act_results;
00341     MsqError err;
00342     MeshUtil tool( &mesh );
00343     tool.edge_length_distribution( act_results, err );
00344     ASSERT_NO_ERROR( err );
00345 }
00346 
00347 void MeshUtilTest::test_lambda_distribution_empty()
00348 {
00349     ArrayMesh empty( 3, 0, 0, 0, 0, 0, 0 );
00350     MeshUtil tool( &empty );
00351     SimpleStats act_results;
00352     MsqError err;
00353     tool.lambda_distribution( act_results, err );
00354     CPPUNIT_ASSERT_EQUAL( MsqError::INVALID_MESH, err.error_code() );
00355 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines