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