Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
#include "TestUtil.hpp"
#include "ElemUtil.hpp"
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include <iostream>
Go to the source code of this file.
Functions | |
void | test_tet () |
void | test_hex () |
void | test_spectral_hex () |
void | test_spectral_quad () |
void | test_spherical_quad () |
void | test_linear_tri () |
void | test_spherical_tri () |
int | main () |
int main | ( | ) |
Definition at line 21 of file ElementTest.cpp.
References RUN_TEST, test_hex(), test_linear_tri(), test_spectral_hex(), test_spectral_quad(), test_spherical_quad(), test_spherical_tri(), and test_tet().
{ int rval = 0; rval += RUN_TEST( test_tet ); rval += RUN_TEST( test_hex ); rval += RUN_TEST( test_spectral_hex ); rval += RUN_TEST( test_spectral_quad ); rval += RUN_TEST( test_spherical_quad ); rval += RUN_TEST( test_linear_tri ); rval += RUN_TEST( test_spherical_tri ); return rval; }
void test_hex | ( | ) |
Definition at line 39 of file ElementTest.cpp.
References moab::Element::Map::ievaluate(), and moab::Element::Map::inside_box().
Referenced by main().
{ double positions[] = { 236.80706050970281, -139.55422526228017, 193.27999999999997, 236.47511729348639, -141.33020962638582, 193.27999999999997, 237.8457938295229, -142.57076074835663, 193.27999999999997, 239.12702305519684, -139.96608933577852, 193.27999999999997, 236.80841321361444, -139.55341321335499, 202.654, 236.47655014713746, -141.32980272396816, 202.654, 237.8477913707564, -142.57047282187165, 202.654, 239.12865103844533, -139.96531051891105, 202.654 }; CartVect x( 235.96518686964933, -142.43503000077749, 188.19999999999987 ); std::vector< CartVect > vertices; for( int i = 0; i < 8; i++ ) vertices.push_back( CartVect( positions + 3 * i ) ); moab::Element::LinearHex hex( vertices ); double tol( 0.0001 ); if( hex.inside_box( x, tol ) ) { CartVect nat_par = hex.ievaluate( x, 0.0001 ); std::cout << nat_par << "\n"; } double positions2[] = { 49.890500000000024, -20.376134375374882, 312.72000000000003, 52.015875000000044, -19.149048546996006, 312.72000000000003, 48.430375821458099, -18.548796774572125, 312.72000000000003, 47.717616239031223, -21.191360829777231, 312.72000000000003, 49.890500000000024, -20.376134375374882, 322.88, 52.015875000000044, -19.149048546996006, 322.88, 48.429930354643275, -18.52828610485021, 322.88, 47.720552036968819, -21.167591146685712, 322.88 }; CartVect x2( 51.469000000000015, -20.145482942833631, 317.80000000000001 ); vertices.clear(); for( int i = 0; i < 8; i++ ) vertices.push_back( CartVect( positions2 + 3 * i ) ); moab::Element::LinearHex hex2( vertices ); if( hex2.inside_box( x2, tol ) ) { try { CartVect nat_par = hex.ievaluate( x, 0.0001 ); std::cout << nat_par << "\n"; } catch( Element::Map::EvaluationError ) { // nothing } } } // test_hex()
void test_linear_tri | ( | ) |
Definition at line 346 of file ElementTest.cpp.
References moab::Element::LinearTri::ievaluate(), and moab::Element::Map::inside_box().
Referenced by main().
{ double positions[] = { 0, 0, 0, 2, 0, 0, 0, 3, 0 }; CartVect x( 1, 0.5, 0 ); std::vector< CartVect > vertices; for( int i = 0; i < 3; i++ ) vertices.push_back( CartVect( positions + 3 * i ) ); moab::Element::LinearTri tri( vertices ); double tol( 0.0001 ); if( tri.inside_box( x, tol ) ) { CartVect nat_par = tri.ievaluate( x ); std::cout << x << " :" << nat_par << "\n"; } x = CartVect( 0, 2, 0 ); if( tri.inside_box( x, tol ) ) { CartVect nat_par = tri.ievaluate( x ); std::cout << x << " :" << nat_par << "\n"; } x = CartVect( 1, 0, 0.5 ); if( tri.inside_box( x, tol ) ) { CartVect nat_par = tri.ievaluate( x ); std::cout << x << " :" << nat_par << "\n"; } double positions2[] = { -0.866026, -0.500001, 0., 0.866026, -0.500001, 0., 0.000000, 100.000000, 0. }; x = CartVect( 0, 0, 0 ); std::vector< CartVect > vertices2; for( int i = 0; i < 3; i++ ) vertices2.push_back( CartVect( positions2 + 3 * i ) ); moab::Element::LinearTri tri2( vertices2 ); if( tri2.inside_box( x, tol ) ) { CartVect nat_par = tri2.ievaluate( x ); std::cout << x << " :" << nat_par << "\n"; } std::cout << "vertices2 " << vertices2[0] << " " << vertices2[1] << " " << vertices2[2] << "\n"; x = CartVect( -0.866026, -0.500001, 0. ); std::cout << x << " :" << tri2.ievaluate( x ) << "\n"; x = CartVect( +0.866026, -0.500001, 0. ); std::cout << x << " :" << tri2.ievaluate( x ) << "\n"; x = CartVect( 0.000000, 100.000000, 0. ); std::cout << x << " :" << tri2.ievaluate( x ) << "\n"; std::cout << "success...\n"; }
void test_spectral_hex | ( | ) |
Definition at line 92 of file ElementTest.cpp.
References moab::Range::begin(), moab::Matrix3::determinant(), moab::Range::empty(), moab::Range::end(), ErrorCode, moab::Element::SpectralHex::evaluate(), moab::Element::SpectralHex::evaluate_scalar_field(), moab::Core::get_entities_by_dimension(), moab::Core::get_entities_by_type_and_tag(), moab::Element::SpectralHex::ievaluate(), moab::Element::SpectralHex::integrate_scalar_field(), inverse(), moab::Element::SpectralHex::jacobian(), moab::Core::load_mesh(), mb, MB_SUCCESS, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, MBENTITYSET, moab::Element::SpectralHex::set_gl_points(), moab::Range::size(), STRINGIFY, moab::Core::tag_get_by_ptr(), moab::Core::tag_get_data(), and moab::Core::tag_get_handle().
Referenced by main().
{ // first load a model that has spectral elements moab::Core* mb = new moab::Core(); std::string meshFile = STRINGIFY( MESHDIR ) "/spectral.h5m"; moab::ErrorCode rval = mb->load_mesh( meshFile.c_str() ); if( moab::MB_SUCCESS != rval ) { std::cout << "Problems reading file " << meshFile << "." << std::endl; delete mb; return; } // get the ent set with SEM_DIMS tag moab::Range spectral_sets; moab::Tag sem_tag; rval = mb->tag_get_handle( "SEM_DIMS", 3, moab::MB_TYPE_INTEGER, sem_tag ); if( moab::MB_SUCCESS != rval ) { std::cout << "can't find tag, no spectral set\n"; delete mb; return; } rval = mb->get_entities_by_type_and_tag( 0, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets ); if( moab::MB_SUCCESS != rval || spectral_sets.empty() ) { std::cout << "can't get sem set\n"; delete mb; return; } moab::Range ents; int sem_dims[3]; moab::EntityHandle firstSemSet = spectral_sets[0]; rval = mb->tag_get_data( sem_tag, &firstSemSet, 1, (void*)sem_dims ); if( moab::MB_SUCCESS != rval ) { delete mb; return; } rval = mb->get_entities_by_dimension( firstSemSet, 3, ents ); if( moab::MB_SUCCESS != rval ) { delete mb; return; } std::cout << "Found " << ents.size() << " " << 3 << "-dimensional entities:" << std::endl; if( sem_dims[0] != sem_dims[1] || sem_dims[0] != sem_dims[2] ) { std::cout << " dimensions are different. bail out\n"; delete mb; return; } // get the SEM_X ...tags moab::Tag xm1Tag, ym1Tag, zm1Tag; int ntot = sem_dims[0] * sem_dims[1] * sem_dims[2]; rval = mb->tag_get_handle( "SEM_X", ntot, moab::MB_TYPE_DOUBLE, xm1Tag ); if( moab::MB_SUCCESS != rval ) { std::cout << "can't get xm1tag \n"; delete mb; return; } rval = mb->tag_get_handle( "SEM_Y", ntot, moab::MB_TYPE_DOUBLE, ym1Tag ); if( moab::MB_SUCCESS != rval ) { std::cout << "can't get ym1tag \n"; delete mb; return; } rval = mb->tag_get_handle( "SEM_Z", ntot, moab::MB_TYPE_DOUBLE, zm1Tag ); if( moab::MB_SUCCESS != rval ) { std::cout << "can't get zm1tag \n"; delete mb; return; } moab::Tag velTag; rval = mb->tag_get_handle( "VX", ntot, moab::MB_TYPE_DOUBLE, velTag ); if( moab::MB_SUCCESS != rval ) { std::cout << "can't get veltag \n"; delete mb; return; } moab::Element::SpectralHex specHex( sem_dims[0] ); // compute the data for some elements for( moab::Range::iterator rit = ents.begin(); rit != ents.end(); ++rit ) { // get the tag pointers to the internal storage for xm1, to not copy the values moab::EntityHandle eh = *rit; const double* xval; const double* yval; const double* zval; rval = mb->tag_get_by_ptr( xm1Tag, &eh, 1, (const void**)&xval ); if( moab::MB_SUCCESS != rval ) { std::cout << "can't get xm1 values \n"; delete mb; return; } rval = mb->tag_get_by_ptr( ym1Tag, &eh, 1, (const void**)&yval ); if( moab::MB_SUCCESS != rval ) { std::cout << "can't get ym1 values \n"; delete mb; return; } rval = mb->tag_get_by_ptr( zm1Tag, &eh, 1, (const void**)&zval ); if( moab::MB_SUCCESS != rval ) { std::cout << "can't get zm1 values \n"; delete mb; return; } if( rit == ents.begin() ) { std::cout << " xm1 for first element: \n"; for( int i = 0; i < ntot; i++ ) std::cout << " " << xval[i]; std::cout << "\n"; } specHex.set_gl_points( (double*)xval, (double*)yval, (double*)zval ); // first evaluate a point, then inverse it to see if we get the same thing moab::CartVect rst( 0.1, -0.1, 0.5 ); moab::CartVect pos = specHex.evaluate( rst ); moab::CartVect inverse = specHex.ievaluate( pos, 0.0001 ); std::cout << "difference" << rst - inverse << "\n"; Matrix3 jac = specHex.jacobian( rst ); std::cout << "jacobian: \n" << jac << " \n determinant: " << jac.determinant() << "\n"; // evaluate vx at some point const double* vx; rval = mb->tag_get_by_ptr( velTag, &eh, 1, (const void**)&vx ); if( moab::MB_SUCCESS != rval ) { std::cout << "can't get vel values \n"; delete mb; return; } double vel1 = specHex.evaluate_scalar_field( rst, vx ); std::cout << "velocity: " << vel1 << "\n"; // compute integral over vx: double integral = specHex.integrate_scalar_field( vx ); std::cout << "integral over vx: " << integral << "\n"; } std::cout << "success...\n"; delete mb; }
void test_spectral_quad | ( | ) |
Definition at line 248 of file ElementTest.cpp.
References moab::Range::begin(), moab::Element::SpectralQuad::compute_gl_positions(), moab::Range::end(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), moab::Core::get_entities_by_type(), moab::Element::SpectralQuad::get_gl_points(), moab::Core::load_mesh(), mb, MB_SUCCESS, MBQUAD, moab::Element::Map::set_vertices(), size, moab::Range::size(), and STRINGIFY.
Referenced by main().
{ // first load a model that has spectral elements moab::Core* mb = new moab::Core(); // use the grid on Sphere from mbcslam folder std::string meshFile = STRINGIFY( MESHDIR ) "/mbcslam/eulerHomme.vtk"; moab::ErrorCode rval = mb->load_mesh( meshFile.c_str() ); if( moab::MB_SUCCESS != rval ) { std::cout << "Problems reading file " << meshFile << "." << std::endl; delete mb; return; } // for each element, compute the gl points and project them on sphere // the radius is the same as the one from intersection test on sphere // double R = 6. * sqrt(3.) / 2; // input moab::Range ents; rval = mb->get_entities_by_type( 0, moab::MBQUAD, ents ); // get all quads if( moab::MB_SUCCESS != rval ) { delete mb; return; } std::cout << "Found " << ents.size() << " " << 2 << "-dimensional entities:" << std::endl; // int NP = 4; // test this.... moab::Element::SpectralQuad specQuad( NP ); // compute the gl points for some elements for( moab::Range::iterator rit = ents.begin(); rit != ents.end(); ++rit ) { const moab::EntityHandle* conn4 = NULL; int num_nodes = 0; rval = mb->get_connectivity( *rit, conn4, num_nodes ); if( moab::MB_SUCCESS != rval ) { std::cout << "can't get connectivity for quad \n"; delete mb; return; } assert( num_nodes == 4 ); std::vector< CartVect > verts( 4 ); rval = mb->get_coords( conn4, num_nodes, &( verts[0][0] ) ); if( moab::MB_SUCCESS != rval ) { std::cout << "can't get coords for quad \n"; delete mb; return; } specQuad.set_vertices( verts ); specQuad.compute_gl_positions(); // do something with the gl positions, project them on a sphere, and create another mesh? if( rit == ents.begin() ) { std::cout << " gl points for first element: \n"; int size; double* xyz[3]; specQuad.get_gl_points( xyz[0], xyz[1], xyz[2], size ); for( int i = 0; i < size; i++ ) std::cout << xyz[0][i] << " " << xyz[1][i] << " " << xyz[2][i] << "\n"; } // project them on a sphere, and create another mesh with it? } std::cout << "success...\n"; delete mb; }
void test_spherical_quad | ( | ) |
Definition at line 323 of file ElementTest.cpp.
References moab::Element::SphericalQuad::ievaluate(), and moab::Element::SphericalQuad::inside_box().
Referenced by main().
{ // example from one coupler test, run like this // ./mbcoupler_test -meshes sphere_16p.h5m mpas_p8.h5m -itag vertex_field -meth 4 -outfile // dd.h5m method 4 is spherical double positions[] = { -0.88882388032987436, -0.069951956448441419, 0.45287838714646161, -0.88226455385461389, -0.13973697758043971, 0.4495362433757738, -0.84497006020160348, -0.13383011007602069, 0.51779831884618843, -0.85072691325794214, -0.066953660115039074, 0.52132612293631853 }; CartVect x( -0.85408569769999998, -0.12391301439999999, 0.50515659540000002 ); std::vector< CartVect > vertices; for( int i = 0; i < 4; i++ ) vertices.push_back( CartVect( positions + 3 * i ) ); moab::Element::SphericalQuad squad( vertices ); double tol( 0.0001 ); if( squad.inside_box( x, tol ) ) { CartVect nat_par = squad.ievaluate( x, 0.000001 ); std::cout << nat_par << "\n"; } std::cout << "success...\n"; }
void test_spherical_tri | ( | ) |
Definition at line 403 of file ElementTest.cpp.
References moab::Element::SphericalTri::ievaluate(), and moab::Element::SphericalTri::inside_box().
Referenced by main().
{ // example from one coupler test, run like this // ./mbcoupler_test -meshes tri_fl_8p.h5m mpas_p8.h5m -itag vertex_field -meth 4 -outfile // oo.h5m -eps 1.e-9 method 4 is spherical double positions[] = { -0.86339258282987197, -0.17004443185241255, 0.47501383044112816, -0.80777478326268271, -0.15172299908552511, 0.5696314870803928, -0.8655618847392077, -0.061613422011313854, 0.49699739427361828 }; CartVect x( -0.85408569769999998, -0.12391301439999999, 0.50515659540000002 ); std::vector< CartVect > vertices; for( int i = 0; i < 3; i++ ) vertices.push_back( CartVect( positions + 3 * i ) ); moab::Element::SphericalTri sphtri( vertices ); double tol( 0.000001 ); if( sphtri.inside_box( x, tol ) ) { CartVect nat_par = sphtri.ievaluate( x, 0.000001 ); std::cout << nat_par << "\n"; } std::cout << "success...\n"; }
void test_tet | ( | ) |
Definition at line 34 of file ElementTest.cpp.
Referenced by main().
{ moab::Element::LinearTet tet; } // test_tet()