MOAB: Mesh Oriented datABase  (version 5.2.1)
elem_eval_test.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines