MOAB: Mesh Oriented datABase  (version 5.3.1)
ElementTest.cpp
Go to the documentation of this file.
00001 #include "TestUtil.hpp"
00002 #include "ElemUtil.hpp"
00003 #include "moab/Core.hpp"
00004 #include "moab/Range.hpp"
00005 #include <iostream>
00006 
00007 #ifndef MESHDIR
00008 #error Specify MESHDIR to compile test
00009 #endif
00010 
00011 using namespace moab;
00012 
00013 void test_tet();
00014 void test_hex();
00015 void test_spectral_hex();
00016 void test_spectral_quad();
00017 void test_spherical_quad();
00018 void test_linear_tri();
00019 void test_spherical_tri();
00020 
00021 int main()
00022 {
00023     int rval = 0;
00024     rval += RUN_TEST( test_tet );
00025     rval += RUN_TEST( test_hex );
00026     rval += RUN_TEST( test_spectral_hex );
00027     rval += RUN_TEST( test_spectral_quad );
00028     rval += RUN_TEST( test_spherical_quad );
00029     rval += RUN_TEST( test_linear_tri );
00030     rval += RUN_TEST( test_spherical_tri );
00031     return rval;
00032 }
00033 
00034 void test_tet()
00035 {
00036     moab::Element::LinearTet tet;
00037 }  // test_tet()
00038 
00039 void test_hex()
00040 {
00041     double positions[] = { 236.80706050970281, -139.55422526228017, 193.27999999999997,
00042                            236.47511729348639, -141.33020962638582, 193.27999999999997,
00043                            237.8457938295229,  -142.57076074835663, 193.27999999999997,
00044                            239.12702305519684, -139.96608933577852, 193.27999999999997,
00045                            236.80841321361444, -139.55341321335499, 202.654,
00046                            236.47655014713746, -141.32980272396816, 202.654,
00047                            237.8477913707564,  -142.57047282187165, 202.654,
00048                            239.12865103844533, -139.96531051891105, 202.654 };
00049     CartVect x( 235.96518686964933, -142.43503000077749, 188.19999999999987 );
00050     std::vector< CartVect > vertices;
00051     for( int i = 0; i < 8; i++ )
00052         vertices.push_back( CartVect( positions + 3 * i ) );
00053 
00054     moab::Element::LinearHex hex( vertices );
00055     double tol( 0.0001 );
00056     if( hex.inside_box( x, tol ) )
00057     {
00058         CartVect nat_par = hex.ievaluate( x, 0.0001 );
00059         std::cout << nat_par << "\n";
00060     }
00061 
00062     double positions2[] = { 49.890500000000024, -20.376134375374882, 312.72000000000003,
00063                             52.015875000000044, -19.149048546996006, 312.72000000000003,
00064                             48.430375821458099, -18.548796774572125, 312.72000000000003,
00065                             47.717616239031223, -21.191360829777231, 312.72000000000003,
00066                             49.890500000000024, -20.376134375374882, 322.88,
00067                             52.015875000000044, -19.149048546996006, 322.88,
00068                             48.429930354643275, -18.52828610485021,  322.88,
00069                             47.720552036968819, -21.167591146685712, 322.88 };
00070 
00071     CartVect x2( 51.469000000000015, -20.145482942833631, 317.80000000000001 );
00072 
00073     vertices.clear();
00074     for( int i = 0; i < 8; i++ )
00075         vertices.push_back( CartVect( positions2 + 3 * i ) );
00076     moab::Element::LinearHex hex2( vertices );
00077     if( hex2.inside_box( x2, tol ) )
00078     {
00079         try
00080         {
00081             CartVect nat_par = hex.ievaluate( x, 0.0001 );
00082             std::cout << nat_par << "\n";
00083         }
00084         catch( Element::Map::EvaluationError )
00085         {
00086             // nothing
00087         }
00088     }
00089 
00090 }  // test_hex()
00091 
00092 void test_spectral_hex()
00093 {
00094     // first load a model that has spectral elements
00095     moab::Core* mb       = new moab::Core();
00096     std::string meshFile = STRINGIFY( MESHDIR ) "/spectral.h5m";
00097     moab::ErrorCode rval = mb->load_mesh( meshFile.c_str() );
00098     if( moab::MB_SUCCESS != rval )
00099     {
00100         std::cout << "Problems reading file " << meshFile << "." << std::endl;
00101         delete mb;
00102         return;
00103     }
00104 
00105     // get the ent set with SEM_DIMS tag
00106     moab::Range spectral_sets;
00107     moab::Tag sem_tag;
00108     rval = mb->tag_get_handle( "SEM_DIMS", 3, moab::MB_TYPE_INTEGER, sem_tag );
00109     if( moab::MB_SUCCESS != rval )
00110     {
00111         std::cout << "can't find tag, no spectral set\n";
00112         delete mb;
00113         return;
00114     }
00115     rval = mb->get_entities_by_type_and_tag( 0, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets );
00116     if( moab::MB_SUCCESS != rval || spectral_sets.empty() )
00117     {
00118         std::cout << "can't get sem set\n";
00119         delete mb;
00120         return;
00121     }
00122 
00123     moab::Range ents;
00124 
00125     int sem_dims[3];
00126     moab::EntityHandle firstSemSet = spectral_sets[0];
00127     rval                           = mb->tag_get_data( sem_tag, &firstSemSet, 1, (void*)sem_dims );
00128     if( moab::MB_SUCCESS != rval )
00129     {
00130         delete mb;
00131         return;
00132     }
00133 
00134     rval = mb->get_entities_by_dimension( firstSemSet, 3, ents );
00135     if( moab::MB_SUCCESS != rval )
00136     {
00137         delete mb;
00138         return;
00139     }
00140     std::cout << "Found " << ents.size() << " " << 3 << "-dimensional entities:" << std::endl;
00141 
00142     if( sem_dims[0] != sem_dims[1] || sem_dims[0] != sem_dims[2] )
00143     {
00144         std::cout << " dimensions are different. bail out\n";
00145         delete mb;
00146         return;
00147     }
00148 
00149     // get the SEM_X ...tags
00150     moab::Tag xm1Tag, ym1Tag, zm1Tag;
00151     int ntot = sem_dims[0] * sem_dims[1] * sem_dims[2];
00152     rval     = mb->tag_get_handle( "SEM_X", ntot, moab::MB_TYPE_DOUBLE, xm1Tag );
00153     if( moab::MB_SUCCESS != rval )
00154     {
00155         std::cout << "can't get xm1tag \n";
00156         delete mb;
00157         return;
00158     }
00159     rval = mb->tag_get_handle( "SEM_Y", ntot, moab::MB_TYPE_DOUBLE, ym1Tag );
00160     if( moab::MB_SUCCESS != rval )
00161     {
00162         std::cout << "can't get ym1tag \n";
00163         delete mb;
00164         return;
00165     }
00166     rval = mb->tag_get_handle( "SEM_Z", ntot, moab::MB_TYPE_DOUBLE, zm1Tag );
00167     if( moab::MB_SUCCESS != rval )
00168     {
00169         std::cout << "can't get zm1tag \n";
00170         delete mb;
00171         return;
00172     }
00173     moab::Tag velTag;
00174 
00175     rval = mb->tag_get_handle( "VX", ntot, moab::MB_TYPE_DOUBLE, velTag );
00176     if( moab::MB_SUCCESS != rval )
00177     {
00178         std::cout << "can't get veltag \n";
00179         delete mb;
00180         return;
00181     }
00182     moab::Element::SpectralHex specHex( sem_dims[0] );
00183 
00184     // compute the data for some elements
00185     for( moab::Range::iterator rit = ents.begin(); rit != ents.end(); ++rit )
00186     {
00187         // get the tag pointers to the internal storage for xm1, to not copy the values
00188         moab::EntityHandle eh = *rit;
00189         const double* xval;
00190         const double* yval;
00191         const double* zval;
00192         rval = mb->tag_get_by_ptr( xm1Tag, &eh, 1, (const void**)&xval );
00193         if( moab::MB_SUCCESS != rval )
00194         {
00195             std::cout << "can't get xm1 values \n";
00196             delete mb;
00197             return;
00198         }
00199         rval = mb->tag_get_by_ptr( ym1Tag, &eh, 1, (const void**)&yval );
00200         if( moab::MB_SUCCESS != rval )
00201         {
00202             std::cout << "can't get ym1 values \n";
00203             delete mb;
00204             return;
00205         }
00206         rval = mb->tag_get_by_ptr( zm1Tag, &eh, 1, (const void**)&zval );
00207         if( moab::MB_SUCCESS != rval )
00208         {
00209             std::cout << "can't get zm1 values \n";
00210             delete mb;
00211             return;
00212         }
00213         if( rit == ents.begin() )
00214         {
00215             std::cout << " xm1 for first element: \n";
00216             for( int i = 0; i < ntot; i++ )
00217                 std::cout << " " << xval[i];
00218             std::cout << "\n";
00219         }
00220         specHex.set_gl_points( (double*)xval, (double*)yval, (double*)zval );
00221         // first evaluate a point, then inverse it to see if we get the same thing
00222         moab::CartVect rst( 0.1, -0.1, 0.5 );
00223         moab::CartVect pos     = specHex.evaluate( rst );
00224         moab::CartVect inverse = specHex.ievaluate( pos, 0.0001 );
00225         std::cout << "difference" << rst - inverse << "\n";
00226         Matrix3 jac = specHex.jacobian( rst );
00227         std::cout << "jacobian: \n" << jac << " \n determinant: " << jac.determinant() << "\n";
00228         // evaluate vx at some point
00229         const double* vx;
00230         rval = mb->tag_get_by_ptr( velTag, &eh, 1, (const void**)&vx );
00231         if( moab::MB_SUCCESS != rval )
00232         {
00233             std::cout << "can't get vel values \n";
00234             delete mb;
00235             return;
00236         }
00237         double vel1 = specHex.evaluate_scalar_field( rst, vx );
00238         std::cout << "velocity: " << vel1 << "\n";
00239         // compute integral over vx:
00240         double integral = specHex.integrate_scalar_field( vx );
00241         std::cout << "integral over vx: " << integral << "\n";
00242     }
00243     std::cout << "success...\n";
00244 
00245     delete mb;
00246 }
00247 
00248 void test_spectral_quad()
00249 {
00250     // first load a model that has spectral elements
00251     moab::Core* mb = new moab::Core();
00252     // use the grid on Sphere from mbcslam folder
00253     std::string meshFile = STRINGIFY( MESHDIR ) "/mbcslam/eulerHomme.vtk";
00254     moab::ErrorCode rval = mb->load_mesh( meshFile.c_str() );
00255     if( moab::MB_SUCCESS != rval )
00256     {
00257         std::cout << "Problems reading file " << meshFile << "." << std::endl;
00258         delete mb;
00259         return;
00260     }
00261 
00262     // for each element, compute the gl points and project them on sphere
00263     // the radius is the same as the one from intersection test on sphere
00264     // double R = 6. * sqrt(3.) / 2; // input
00265 
00266     moab::Range ents;
00267 
00268     rval = mb->get_entities_by_type( 0, moab::MBQUAD, ents );  // get all quads
00269     if( moab::MB_SUCCESS != rval )
00270     {
00271         delete mb;
00272         return;
00273     }
00274     std::cout << "Found " << ents.size() << " " << 2 << "-dimensional entities:" << std::endl;
00275 
00276     //
00277     int NP = 4;  // test this....
00278     moab::Element::SpectralQuad specQuad( NP );
00279 
00280     // compute the gl points for some elements
00281     for( moab::Range::iterator rit = ents.begin(); rit != ents.end(); ++rit )
00282     {
00283 
00284         const moab::EntityHandle* conn4 = NULL;
00285         int num_nodes                   = 0;
00286         rval                            = mb->get_connectivity( *rit, conn4, num_nodes );
00287         if( moab::MB_SUCCESS != rval )
00288         {
00289             std::cout << "can't get connectivity for quad \n";
00290             delete mb;
00291             return;
00292         }
00293         assert( num_nodes == 4 );
00294 
00295         std::vector< CartVect > verts( 4 );
00296         rval = mb->get_coords( conn4, num_nodes, &( verts[0][0] ) );
00297         if( moab::MB_SUCCESS != rval )
00298         {
00299             std::cout << "can't get coords for quad \n";
00300             delete mb;
00301             return;
00302         }
00303 
00304         specQuad.set_vertices( verts );
00305         specQuad.compute_gl_positions();
00306         // do something with the gl positions, project them on a sphere, and create another mesh?
00307         if( rit == ents.begin() )
00308         {
00309             std::cout << " gl points for first element: \n";
00310             int size;
00311             double* xyz[3];
00312             specQuad.get_gl_points( xyz[0], xyz[1], xyz[2], size );
00313             for( int i = 0; i < size; i++ )
00314                 std::cout << xyz[0][i] << " " << xyz[1][i] << " " << xyz[2][i] << "\n";
00315         }
00316 
00317         // project them on a sphere, and create another mesh with it?
00318     }
00319     std::cout << "success...\n";
00320 
00321     delete mb;
00322 }
00323 void test_spherical_quad()
00324 {
00325     // example from one coupler test, run like this
00326     // ./mbcoupler_test -meshes sphere_16p.h5m mpas_p8.h5m -itag vertex_field -meth 4 -outfile
00327     // dd.h5m method 4 is spherical
00328     double positions[] = { -0.88882388032987436, -0.069951956448441419, 0.45287838714646161,   -0.88226455385461389,
00329                            -0.13973697758043971, 0.4495362433757738,    -0.84497006020160348,  -0.13383011007602069,
00330                            0.51779831884618843,  -0.85072691325794214,  -0.066953660115039074, 0.52132612293631853 };
00331     CartVect x( -0.85408569769999998, -0.12391301439999999, 0.50515659540000002 );
00332     std::vector< CartVect > vertices;
00333     for( int i = 0; i < 4; i++ )
00334         vertices.push_back( CartVect( positions + 3 * i ) );
00335 
00336     moab::Element::SphericalQuad squad( vertices );
00337     double tol( 0.0001 );
00338     if( squad.inside_box( x, tol ) )
00339     {
00340         CartVect nat_par = squad.ievaluate( x, 0.000001 );
00341         std::cout << nat_par << "\n";
00342     }
00343     std::cout << "success...\n";
00344 }
00345 
00346 void test_linear_tri()
00347 {
00348     double positions[] = { 0, 0, 0, 2, 0, 0, 0, 3, 0 };
00349     CartVect x( 1, 0.5, 0 );
00350     std::vector< CartVect > vertices;
00351     for( int i = 0; i < 3; i++ )
00352         vertices.push_back( CartVect( positions + 3 * i ) );
00353 
00354     moab::Element::LinearTri tri( vertices );
00355     double tol( 0.0001 );
00356     if( tri.inside_box( x, tol ) )
00357     {
00358         CartVect nat_par = tri.ievaluate( x );
00359         std::cout << x << " :" << nat_par << "\n";
00360     }
00361 
00362     x = CartVect( 0, 2, 0 );
00363     if( tri.inside_box( x, tol ) )
00364     {
00365         CartVect nat_par = tri.ievaluate( x );
00366         std::cout << x << " :" << nat_par << "\n";
00367     }
00368 
00369     x = CartVect( 1, 0, 0.5 );
00370     if( tri.inside_box( x, tol ) )
00371     {
00372         CartVect nat_par = tri.ievaluate( x );
00373         std::cout << x << " :" << nat_par << "\n";
00374     }
00375 
00376     double positions2[] = { -0.866026, -0.500001, 0., 0.866026, -0.500001, 0., 0.000000, 100.000000, 0. };
00377     x                   = CartVect( 0, 0, 0 );
00378     std::vector< CartVect > vertices2;
00379     for( int i = 0; i < 3; i++ )
00380         vertices2.push_back( CartVect( positions2 + 3 * i ) );
00381 
00382     moab::Element::LinearTri tri2( vertices2 );
00383 
00384     if( tri2.inside_box( x, tol ) )
00385     {
00386         CartVect nat_par = tri2.ievaluate( x );
00387         std::cout << x << " :" << nat_par << "\n";
00388     }
00389 
00390     std::cout << "vertices2 " << vertices2[0] << " " << vertices2[1] << " " << vertices2[2] << "\n";
00391 
00392     x = CartVect( -0.866026, -0.500001, 0. );
00393     std::cout << x << " :" << tri2.ievaluate( x ) << "\n";
00394 
00395     x = CartVect( +0.866026, -0.500001, 0. );
00396     std::cout << x << " :" << tri2.ievaluate( x ) << "\n";
00397     x = CartVect( 0.000000, 100.000000, 0. );
00398     std::cout << x << " :" << tri2.ievaluate( x ) << "\n";
00399 
00400     std::cout << "success...\n";
00401 }
00402 
00403 void test_spherical_tri()
00404 {
00405     // example from one coupler test, run like this
00406     // ./mbcoupler_test -meshes  tri_fl_8p.h5m mpas_p8.h5m -itag vertex_field -meth 4  -outfile
00407     // oo.h5m -eps 1.e-9 method 4 is spherical
00408     double positions[] = { -0.86339258282987197, -0.17004443185241255,  0.47501383044112816,
00409                            -0.80777478326268271, -0.15172299908552511,  0.5696314870803928,
00410                            -0.8655618847392077,  -0.061613422011313854, 0.49699739427361828 };
00411     CartVect x( -0.85408569769999998, -0.12391301439999999, 0.50515659540000002 );
00412     std::vector< CartVect > vertices;
00413     for( int i = 0; i < 3; i++ )
00414         vertices.push_back( CartVect( positions + 3 * i ) );
00415 
00416     moab::Element::SphericalTri sphtri( vertices );
00417     double tol( 0.000001 );
00418     if( sphtri.inside_box( x, tol ) )
00419     {
00420         CartVect nat_par = sphtri.ievaluate( x, 0.000001 );
00421         std::cout << nat_par << "\n";
00422     }
00423     std::cout << "success...\n";
00424 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines