MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 }