ElementTest.cpp

Go to the documentation of this file.

#include "TestUtil.hpp"
#include "ElemUtil.hpp"
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include <iostream>
 
#ifndef MESHDIR
#error Specify MESHDIR to compile test
#endif
 
using namespace moab;
 
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 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_tet()
{
 moab::Element::LinearTet tet;
} // test_tet()
 
void test_hex()
{
 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_spectral_hex()
{
 // 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()
{
 // 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()
{
 // 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_linear_tri()
{
 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_spherical_tri()
{
 // 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";
}