MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /** 00002 * \file elem_eval_test.cpp 00003 * 00004 * \brief test ElemEvaluator and the various element types in MOAB 00005 * 00006 */ 00007 #include "moab/Core.hpp" 00008 #include "moab/ReadUtilIface.hpp" 00009 #include "moab/ElemEvaluator.hpp" 00010 #include "moab/LocalDiscretization/LinearTri.hpp" 00011 #include "moab/LocalDiscretization/LinearQuad.hpp" 00012 #include "moab/LocalDiscretization/LinearHex.hpp" 00013 #include "moab/LocalDiscretization/LinearTet.hpp" 00014 #include "moab/LocalDiscretization/QuadraticHex.hpp" 00015 #include "moab/CartVect.hpp" 00016 #include "TestUtil.hpp" 00017 00018 using namespace moab; 00019 00020 void test_linear_tri(); 00021 void test_linear_quad(); 00022 void test_linear_hex(); 00023 void test_linear_tet(); 00024 void test_quadratic_hex(); 00025 void test_normal_linear_tri(); 00026 void test_normal_linear_quad(); 00027 void test_normal_linear_tet(); 00028 void test_normal_linear_hex(); 00029 ErrorCode create_mesh( Core& mb, EntityType type ); 00030 00031 CartVect hex_verts[] = { 00032 // corners 00033 CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( 1, 1, -1 ), CartVect( -1, 1, -1 ), CartVect( -1, -1, 1 ), 00034 CartVect( 1, -1, 1 ), CartVect( 1, 1, 1 ), CartVect( -1, 1, 1 ), 00035 // mid-edge (bottom, middle, top) 00036 CartVect( 0, -1, -1 ), CartVect( 1, 0, -1 ), CartVect( 0, 1, -1 ), CartVect( -1, 0, -1 ), CartVect( -1, -1, 0 ), 00037 CartVect( 1, -1, 0 ), CartVect( 1, 1, 0 ), CartVect( -1, 1, 0 ), CartVect( 0, -1, 1 ), CartVect( 1, 0, 1 ), 00038 CartVect( 0, 1, 1 ), CartVect( -1, 0, 1 ), 00039 // mid-face (middle, bottom, top) 00040 CartVect( 0, -1, 0 ), CartVect( 1, 0, 0 ), CartVect( 0, 1, 0 ), CartVect( -1, 0, 0 ), CartVect( 0, 0, -1 ), 00041 CartVect( 0, 0, 1 ), 00042 // mid-element 00043 CartVect( 0, 0, 0 ) 00044 }; 00045 00046 const double EPS1 = 1.0e-6; 00047 00048 void test_eval( ElemEvaluator& ee, bool test_integrate ) 00049 { 00050 00051 CartVect params, posn, params2; 00052 int is_inside; 00053 Matrix3 jacob; 00054 ErrorCode rval; 00055 int ent_dim = ee.get_moab()->dimension_from_handle( ee.get_ent_handle() ); 00056 00057 for( params[0] = -1; params[0] <= 1; params[0] += 0.2 ) 00058 { 00059 for( params[1] = -1; params[1] <= 1; params[1] += 0.2 ) 00060 { 00061 for( params[2] = -1; params[2] <= 1; params[2] += 0.2 ) 00062 { 00063 00064 // forward/reverse evaluation should get back to the same point, within tol 00065 if( !ee.inside( params.array(), EPS1 ) ) continue; 00066 00067 rval = ee.eval( params.array(), posn.array() );CHECK_ERR( rval ); 00068 rval = ee.reverse_eval( posn.array(), EPS1, EPS1, params2.array(), &is_inside );CHECK_ERR( rval ); 00069 CHECK_REAL_EQUAL( 0.0, params[0] - params2[0], 3 * EPS1 ); 00070 if( ent_dim > 1 ) CHECK_REAL_EQUAL( 0.0, params[1] - params2[1], 3 * EPS1 ); 00071 if( ent_dim > 2 ) CHECK_REAL_EQUAL( 0.0, params[2] - params2[2], 3 * EPS1 ); 00072 00073 // jacobian should be >= 0 00074 rval = ee.jacobian( params.array(), jacob.array() );CHECK_ERR( rval ); 00075 CHECK( jacob.determinant() >= 0.0 ); 00076 } 00077 } 00078 } 00079 00080 if( !test_integrate ) return; 00081 00082 // tag equal to coordinates should integrate to avg position, test that 00083 Tag tag; 00084 // make a temporary tag and set it on vertices to the vertex positions 00085 rval = ee.get_moab()->tag_get_handle( NULL, 3, MB_TYPE_DOUBLE, tag, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval ); 00086 rval = ee.get_moab()->tag_set_data( tag, ee.get_vert_handles(), ee.get_num_verts(), ee.get_vert_pos() );CHECK_ERR( rval ); 00087 // set that temporary tag on the evaluator so that's what gets integrated 00088 rval = ee.set_tag_handle( tag, 0 );CHECK_ERR( rval ); 00089 00090 CartVect integral, avg; 00091 rval = ee.integrate( integral.array() );CHECK_ERR( rval ); 00092 00093 // now integrate a const 1-valued function, using direct call to the integrate function 00094 std::vector< double > one( ee.get_num_verts(), 1.0 ); 00095 double measure; 00096 EvalSet es; 00097 EntityHandle eh = ee.get_ent_handle(); 00098 rval = EvalSet::get_eval_set( ee.get_moab(), eh, es );CHECK_ERR( rval ); 00099 rval = ( *es.integrateFcn )( &one[0], ee.get_vert_pos(), ee.get_num_verts(), 00100 ee.get_moab()->dimension_from_handle( eh ), 1, ee.get_work_space(), &measure );CHECK_ERR( rval ); 00101 if( measure ) integral /= measure; 00102 00103 // check against avg of entity's vertices' positions 00104 rval = ee.get_moab()->get_coords( &eh, 1, avg.array() );CHECK_ERR( rval ); 00105 CHECK_REAL_EQUAL( 0.0, ( avg - integral ).length(), EPS1 ); 00106 00107 // delete the temporary tag 00108 rval = ee.get_moab()->tag_delete( tag );CHECK_ERR( rval ); 00109 // set the ee's tag back to coords 00110 rval = ee.set_tag_handle( 0, 0 );CHECK_ERR( rval ); 00111 } 00112 00113 void test_evals( ElemEvaluator& ee, bool test_integrate, EntityHandle* ents, int num_ents, Tag onetag, 00114 double total_vol ) 00115 { 00116 for( int i = 0; i < num_ents; i++ ) 00117 { 00118 ee.set_ent_handle( ents[i] ); 00119 test_eval( ee, false ); 00120 } 00121 00122 if( !test_integrate ) return; 00123 if( !num_ents || !ents ) CHECK_ERR( MB_FAILURE ); 00124 00125 ErrorCode rval = ee.set_tag_handle( onetag, 0 );CHECK_ERR( rval ); 00126 00127 double tot_vol = 0.0; 00128 for( int i = 0; i < num_ents; i++ ) 00129 { 00130 double tmp_vol; 00131 ee.set_ent_handle( ents[i] ); 00132 rval = ee.integrate( &tmp_vol );CHECK_ERR( rval ); 00133 tot_vol += tmp_vol; 00134 } 00135 CHECK_REAL_EQUAL( total_vol, tot_vol, EPS1 ); 00136 } 00137 00138 int main() 00139 { 00140 int failures = 0; 00141 00142 failures += RUN_TEST( test_linear_tri ); // currently failing linear tri, bad formulation, working on it... 00143 failures += RUN_TEST( test_linear_quad ); 00144 failures += RUN_TEST( test_linear_hex ); 00145 failures += RUN_TEST( test_quadratic_hex ); 00146 failures += RUN_TEST( test_linear_tet ); 00147 failures += RUN_TEST( test_normal_linear_tri ); 00148 failures += RUN_TEST( test_normal_linear_quad ); 00149 failures += RUN_TEST( test_normal_linear_tet ); 00150 failures += RUN_TEST( test_normal_linear_hex ); 00151 00152 return failures; 00153 } 00154 00155 void test_linear_tri() 00156 { 00157 Core mb; 00158 Range verts; 00159 double tri_verts[] = { -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0 }; 00160 00161 ErrorCode rval = mb.create_vertices( tri_verts, 3, verts );CHECK_ERR( rval ); 00162 EntityHandle tri; 00163 std::vector< EntityHandle > connect; 00164 std::copy( verts.begin(), verts.end(), std::back_inserter( connect ) ); 00165 rval = mb.create_element( MBTRI, &connect[0], 3, tri );CHECK_ERR( rval ); 00166 00167 ElemEvaluator ee( &mb, tri, 0 ); 00168 ee.set_tag_handle( 0, 0 ); 00169 ee.set_eval_set( MBTRI, LinearTri::eval_set() ); 00170 00171 test_eval( ee, true ); 00172 } 00173 00174 void test_linear_quad() 00175 { 00176 Core mb; 00177 Range verts; 00178 ErrorCode rval = mb.create_vertices( (double*)hex_verts, 4, verts );CHECK_ERR( rval ); 00179 EntityHandle quad; 00180 std::vector< EntityHandle > connect; 00181 std::copy( verts.begin(), verts.end(), std::back_inserter( connect ) ); 00182 rval = mb.create_element( MBQUAD, &connect[0], 4, quad );CHECK_ERR( rval ); 00183 00184 ElemEvaluator ee( &mb, quad, 0 ); 00185 ee.set_tag_handle( 0, 0 ); 00186 ee.set_eval_set( MBQUAD, LinearQuad::eval_set() ); 00187 00188 test_eval( ee, true ); 00189 } 00190 00191 void test_linear_hex() 00192 { 00193 Core mb; 00194 Range verts; 00195 ErrorCode rval = mb.create_vertices( (double*)hex_verts, 8, verts );CHECK_ERR( rval ); 00196 EntityHandle hex; 00197 std::vector< EntityHandle > connect; 00198 std::copy( verts.begin(), verts.end(), std::back_inserter( connect ) ); 00199 rval = mb.create_element( MBHEX, &connect[0], 8, hex );CHECK_ERR( rval ); 00200 00201 ElemEvaluator ee( &mb, hex, 0 ); 00202 ee.set_tag_handle( 0, 0 ); 00203 ee.set_eval_set( MBHEX, LinearHex::eval_set() ); 00204 00205 test_eval( ee, true ); 00206 } 00207 00208 void test_quadratic_hex() 00209 { 00210 Core mb; 00211 Range verts; 00212 ErrorCode rval = mb.create_vertices( (double*)hex_verts, 27, verts );CHECK_ERR( rval ); 00213 EntityHandle hex; 00214 std::vector< EntityHandle > connect; 00215 std::copy( verts.begin(), verts.end(), std::back_inserter( connect ) ); 00216 rval = mb.create_element( MBHEX, &connect[0], 27, hex );CHECK_ERR( rval ); 00217 00218 ElemEvaluator ee( &mb, hex, 0 ); 00219 ee.set_tag_handle( 0, 0 ); 00220 ee.set_eval_set( MBHEX, QuadraticHex::eval_set() ); 00221 00222 test_eval( ee, false ); 00223 } 00224 00225 void test_linear_tet() 00226 { 00227 Core mb; 00228 Range verts; 00229 ErrorCode rval = mb.create_vertices( (double*)hex_verts[0].array(), 8, verts );CHECK_ERR( rval ); 00230 EntityHandle starth = 1, *conn; 00231 int conn_inds[] = { 1, 6, 4, 5, 1, 4, 6, 3, 0, 1, 3, 4, 1, 2, 3, 6, 3, 4, 6, 7 }; 00232 ReadUtilIface* ru; 00233 rval = mb.query_interface( ru );CHECK_ERR( rval ); 00234 rval = ru->get_element_connect( 5, 4, MBTET, 1, starth, conn );CHECK_ERR( rval ); 00235 for( unsigned int i = 0; i < 20; i++ ) 00236 conn[i] = verts[conn_inds[i]]; 00237 EntityHandle tets[5]; 00238 for( unsigned int i = 0; i < 5; i++ ) 00239 tets[i] = starth + i; 00240 ElemEvaluator ee( &mb, 0, 0 ); 00241 ee.set_tag_handle( 0, 0 ); 00242 ee.set_eval_set( MBTET, LinearTet::eval_set() ); 00243 00244 // make a tag whose value is one on all vertices 00245 Tag tag; 00246 rval = mb.tag_get_handle( NULL, 1, MB_TYPE_DOUBLE, tag, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval ); 00247 std::vector< double > vals( verts.size(), 1.0 ); 00248 rval = mb.tag_set_data( tag, verts, &vals[0] );CHECK_ERR( rval ); 00249 00250 test_evals( ee, true, tets, 5, tag, 8.0 ); 00251 } 00252 00253 void test_normal_linear_tri() 00254 { 00255 ErrorCode error; 00256 Core mb; 00257 00258 error = create_mesh( mb, MBTRI );CHECK_ERR( error ); 00259 Range faces; 00260 error = mb.get_entities_by_dimension( 0, 2, faces );CHECK_ERR( error ); 00261 00262 ElemEvaluator ee( &mb, 0, 0 ); 00263 ee.set_eval_set( MBTRI, LinearTri::eval_set() ); 00264 00265 double nrms[8][3]; 00266 00267 for( int i = 0; i < 4; i++ ) 00268 { 00269 ee.set_ent_handle( faces[i] ); 00270 ee.get_normal( 1, 0, nrms[2 * i] ); 00271 ee.get_normal( 1, 2, nrms[2 * i + 1] ); 00272 } 00273 00274 for( int i = 0; i < 4; i++ ) 00275 { 00276 if( i == 3 ) 00277 { 00278 double val = nrms[7][0] * nrms[0][0] + nrms[7][1] * nrms[0][1] + nrms[7][2] * nrms[0][2]; 00279 CHECK_EQUAL( val, -1 ); 00280 } 00281 else 00282 { 00283 double val = nrms[2 * i + 1][0] * nrms[2 * ( i + 1 )][0] + nrms[2 * i + 1][1] * nrms[2 * ( i + 1 )][1] + 00284 nrms[2 * i + 1][2] * nrms[2 * ( i + 1 )][2]; 00285 CHECK_EQUAL( val, -1 ); 00286 } 00287 } 00288 } 00289 00290 void test_normal_linear_quad() 00291 { 00292 ErrorCode error; 00293 Core mb; 00294 00295 error = create_mesh( mb, MBQUAD );CHECK_ERR( error ); 00296 Range faces; 00297 error = mb.get_entities_by_dimension( 0, 2, faces );CHECK_ERR( error ); 00298 00299 ElemEvaluator ee( &mb, 0, 0 ); 00300 ee.set_eval_set( MBQUAD, LinearQuad::eval_set() ); 00301 00302 double nrms[8][3]; 00303 00304 for( int i = 0; i < 4; i++ ) 00305 { 00306 ee.set_ent_handle( faces[i] ); 00307 ee.get_normal( 1, 0, nrms[2 * i] ); 00308 ee.get_normal( 1, 3, nrms[2 * i + 1] ); 00309 } 00310 00311 for( int i = 0; i < 4; i++ ) 00312 { 00313 if( i == 3 ) 00314 { 00315 double val = nrms[7][0] * nrms[0][0] + nrms[7][1] * nrms[0][1] + nrms[7][2] * nrms[0][2]; 00316 CHECK_EQUAL( val, -1 ); 00317 } 00318 else 00319 { 00320 double val = nrms[2 * i + 1][0] * nrms[2 * ( i + 1 )][0] + nrms[2 * i + 1][1] * nrms[2 * ( i + 1 )][1] + 00321 nrms[2 * i + 1][2] * nrms[2 * ( i + 1 )][2]; 00322 CHECK_EQUAL( val, -1 ); 00323 } 00324 } 00325 } 00326 00327 void test_normal_linear_tet() 00328 { 00329 ErrorCode error; 00330 Core mb; 00331 00332 error = create_mesh( mb, MBTET );CHECK_ERR( error ); 00333 Range cells; 00334 error = mb.get_entities_by_dimension( 0, 3, cells );CHECK_ERR( error ); 00335 00336 ElemEvaluator ee( &mb, 0, 0 ); 00337 ee.set_eval_set( MBTET, LinearTet::eval_set() ); 00338 00339 double nrms[8][3]; 00340 00341 for( int i = 0; i < 4; i++ ) 00342 { 00343 ee.set_ent_handle( cells[i] ); 00344 ee.get_normal( 2, 0, nrms[2 * i] ); 00345 ee.get_normal( 2, 2, nrms[2 * i + 1] ); 00346 } 00347 00348 for( int i = 0; i < 4; i++ ) 00349 { 00350 if( i == 3 ) 00351 { 00352 double val = nrms[7][0] * nrms[0][0] + nrms[7][1] * nrms[0][1] + nrms[7][2] * nrms[0][2]; 00353 CHECK_EQUAL( val, -1 ); 00354 } 00355 else 00356 { 00357 double val = nrms[2 * i + 1][0] * nrms[2 * ( i + 1 )][0] + nrms[2 * i + 1][1] * nrms[2 * ( i + 1 )][1] + 00358 nrms[2 * i + 1][2] * nrms[2 * ( i + 1 )][2]; 00359 CHECK_EQUAL( val, -1 ); 00360 } 00361 } 00362 } 00363 00364 void test_normal_linear_hex() 00365 { 00366 ErrorCode error; 00367 Core mb; 00368 00369 error = create_mesh( mb, MBHEX );CHECK_ERR( error ); 00370 Range cells; 00371 error = mb.get_entities_by_dimension( 0, 3, cells );CHECK_ERR( error ); 00372 00373 ElemEvaluator ee( &mb, 0, 0 ); 00374 ee.set_eval_set( MBHEX, LinearHex::eval_set() ); 00375 00376 double nrms[8][3]; 00377 00378 for( int i = 0; i < 4; i++ ) 00379 { 00380 ee.set_ent_handle( cells[i] ); 00381 ee.get_normal( 2, 0, nrms[2 * i] ); 00382 ee.get_normal( 2, 3, nrms[2 * i + 1] ); 00383 } 00384 00385 for( int i = 0; i < 4; i++ ) 00386 { 00387 if( i == 3 ) 00388 { 00389 double val = nrms[7][0] * nrms[0][0] + nrms[7][1] * nrms[0][1] + nrms[7][2] * nrms[0][2]; 00390 CHECK_EQUAL( val, -1 ); 00391 } 00392 else 00393 { 00394 double val = nrms[2 * i + 1][0] * nrms[2 * ( i + 1 )][0] + nrms[2 * i + 1][1] * nrms[2 * ( i + 1 )][1] + 00395 nrms[2 * i + 1][2] * nrms[2 * ( i + 1 )][2]; 00396 CHECK_EQUAL( val, -1 ); 00397 } 00398 } 00399 } 00400 00401 ErrorCode create_mesh( Core& mb, EntityType type ) 00402 { 00403 ErrorCode error; 00404 if( type == MBTRI ) 00405 { 00406 const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0 }; 00407 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00408 00409 const int conn[] = { 0, 1, 2, 0, 2, 3, 0, 3, 4, 0, 4, 1 }; 00410 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3; 00411 00412 EntityHandle verts[num_vtx], faces[num_elems]; 00413 for( size_t i = 0; i < num_vtx; ++i ) 00414 { 00415 error = mb.create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00416 } 00417 00418 for( size_t i = 0; i < num_elems; ++i ) 00419 { 00420 EntityHandle c[3]; 00421 for( int j = 0; j < 3; j++ ) 00422 c[j] = verts[conn[3 * i + j]]; 00423 00424 error = mb.create_element( MBTRI, c, 3, faces[i] );CHECK_ERR( error ); 00425 } 00426 } 00427 else if( type == MBQUAD ) 00428 { 00429 const double coords[] = { 00430 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, -1, 1, 0, -1, 0, 0, -1, -1, 0, 0, -1, 0, 1, -1, 0 00431 }; 00432 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00433 00434 const int conn[] = { 0, 1, 2, 3, 0, 3, 4, 5, 0, 5, 6, 7, 0, 7, 8, 1 }; 00435 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4; 00436 00437 EntityHandle verts[num_vtx], faces[num_elems]; 00438 for( size_t i = 0; i < num_vtx; ++i ) 00439 { 00440 error = mb.create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00441 } 00442 00443 for( size_t i = 0; i < num_elems; ++i ) 00444 { 00445 EntityHandle c[4]; 00446 for( int j = 0; j < 4; j++ ) 00447 c[j] = verts[conn[4 * i + j]]; 00448 00449 error = mb.create_element( MBQUAD, c, 4, faces[i] );CHECK_ERR( error ); 00450 } 00451 } 00452 else if( type == MBTET ) 00453 { 00454 const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 1 }; 00455 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00456 00457 const int conn[] = { 0, 1, 2, 5, 0, 2, 3, 5, 0, 3, 4, 5, 0, 4, 1, 5 }; 00458 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4; 00459 00460 EntityHandle verts[num_vtx], cells[num_elems]; 00461 for( size_t i = 0; i < num_vtx; ++i ) 00462 { 00463 error = mb.create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00464 } 00465 00466 for( size_t i = 0; i < num_elems; ++i ) 00467 { 00468 EntityHandle c[4]; 00469 for( int j = 0; j < 4; j++ ) 00470 c[j] = verts[conn[4 * i + j]]; 00471 00472 error = mb.create_element( MBTET, c, 4, cells[i] );CHECK_ERR( error ); 00473 } 00474 } 00475 else if( type == MBHEX ) 00476 { 00477 const double coords[] = { 00478 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, -1, 1, 0, -1, 0, 0, -1, -1, 0, 0, -1, 0, 1, -1, 0, 00479 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, -1, 1, 1, -1, 0, 1, -1, -1, 1, 0, -1, 1, 1, -1, 1 00480 }; 00481 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00482 00483 const int conn[] = { 0, 1, 2, 3, 9, 10, 11, 12, 0, 3, 4, 5, 9, 12, 13, 14, 00484 0, 5, 6, 7, 9, 14, 15, 16, 0, 7, 8, 1, 9, 16, 17, 10 }; 00485 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 8; 00486 00487 EntityHandle verts[num_vtx], cells[num_elems]; 00488 for( size_t i = 0; i < num_vtx; ++i ) 00489 { 00490 error = mb.create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00491 } 00492 00493 for( size_t i = 0; i < num_elems; ++i ) 00494 { 00495 EntityHandle c[8]; 00496 for( int j = 0; j < 8; j++ ) 00497 c[j] = verts[conn[8 * i + j]]; 00498 00499 error = mb.create_element( MBHEX, c, 8, cells[i] );CHECK_ERR( error ); 00500 } 00501 } 00502 00503 return MB_SUCCESS; 00504 }