|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#include <iostream>#include <sstream>#include <ctime>#include <cstdlib>#include <cstdio>#include <cstring>#include "moab/Core.hpp"#include "moab/Interface.hpp"#include "moab/IntxMesh/Intx2MeshOnSphere.hpp"#include <cmath>#include "TestUtil.hpp"#include "moab/ParallelComm.hpp"#include "moab/ProgOptions.hpp"#include "MBParallelConventions.h"#include "moab/ReadUtilIface.hpp"#include "MBTagConventions.hpp"#include "moab/IntxMesh/IntxUtils.hpp"#include "IntxUtilsCSLAM.hpp"
Include dependency graph for cslam_par_test.cpp:Go to the source code of this file.
Functions | |
| std::string | input_mesh_file ("VELO00_16p.h5m") |
| void | test_intx_in_parallel_elem_based () |
| int | main (int argc, char **argv) |
| ErrorCode | compute_lagrange_mesh_on_sphere (Interface *mb, EntityHandle euler_set) |
Variables | |
| double | EPS1 = 0.2 |
| double | Radius = 1.0 |
| double | deltaT = 1.e-6 |
| ErrorCode compute_lagrange_mesh_on_sphere | ( | Interface * | mb, |
| EntityHandle | euler_set | ||
| ) |
Definition at line 86 of file cslam_par_test.cpp.
References moab::Range::begin(), deltaT, moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_type(), moab::CartVect::length(), MB_CHK_ERR, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MBQUAD, Radius, moab::Range::size(), moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), and moab::Interface::tag_iterate().
Referenced by test_intx_in_parallel_elem_based().
{
/*
* get all quads first, then vertices, then move them on the surface of the sphere
* radius is 1, usually
* pos (t-dt) = pos(t) -Velo(t)*dt; this will be lagrange mesh, on each processor
*/
Range quads;
ErrorCode rval = mb->get_entities_by_type( euler_set, MBQUAD, quads );MB_CHK_ERR( rval );
Range connecVerts;
rval = mb->get_connectivity( quads, connecVerts );MB_CHK_ERR( rval );
// the LOC tag, should be provided by the user?
Tag tagh = 0;
std::string tag_name( "DP" );
rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
void* data; // pointer to the DP in memory, for each vertex
int count;
rval = mb->tag_iterate( tagh, connecVerts.begin(), connecVerts.end(), count, data );MB_CHK_ERR( rval );
// here we are checking contiguity
assert( count == (int)connecVerts.size() );
double* ptr_DP = (double*)data;
// get the coordinates of the old mesh, and move it around using velocity tag
Tag tagv = 0;
std::string velo_tag_name( "VELO" );
rval = mb->tag_get_handle( velo_tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagv, MB_TAG_DENSE );MB_CHK_ERR( rval );
/*void *datavelo; // pointer to the VELO in memory, for each vertex
rval = mb->tag_iterate(tagv, connecVerts.begin(), connecVerts.end(), count, datavelo);MB_CHK_ERR(rval);*/
// here we are checking contiguity
assert( count == (int)connecVerts.size() );
// now put the vertices in the right place....
// int vix=0; // vertex index in new array
for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
{
EntityHandle oldV = *vit;
CartVect posi;
rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
CartVect velo;
rval = mb->tag_get_data( tagv, &oldV, 1, (void*)&( velo[0] ) );MB_CHK_ERR( rval );
// do some mumbo jumbo, as in python script
CartVect newPos = posi - deltaT * velo;
double len1 = newPos.length();
newPos = Radius * newPos / len1;
ptr_DP[0] = newPos[0];
ptr_DP[1] = newPos[1];
ptr_DP[2] = newPos[2];
ptr_DP += 3; // increment to the next node
}
return rval;
}
| std::string input_mesh_file | ( | "VELO00_16p.h5m" | ) |
Referenced by main(), and test_intx_in_parallel_elem_based().
| int main | ( | int | argc, |
| char ** | argv | ||
| ) |
Definition at line 47 of file cslam_par_test.cpp.
References deltaT, EPS1, input_mesh_file(), Radius, RUN_TEST, and test_intx_in_parallel_elem_based().
{
MPI_Init( &argc, &argv );
EPS1 = 0.000002;
int result = 0;
if( argc > 1 )
{
int index = 1;
while( index < argc )
{
if( !strcmp( argv[index], "-eps" ) ) // this is for box error
{
EPS1 = atof( argv[++index] );
}
if( !strcmp( argv[index], "-input" ) )
{
input_mesh_file = argv[++index];
}
if( !strcmp( argv[index], "-radius" ) )
{
Radius = atof( argv[++index] );
}
if( !strcmp( argv[index], "-deltaT" ) )
{
deltaT = atof( argv[++index] );
}
index++;
}
}
std::cout << " run: -input " << input_mesh_file << " -eps " << EPS1 << " -radius " << Radius << " -deltaT "
<< deltaT << "\n";
result += RUN_TEST( test_intx_in_parallel_elem_based );
MPI_Finalize();
return result;
}
| void test_intx_in_parallel_elem_based | ( | ) |
Definition at line 145 of file cslam_par_test.cpp.
References moab::IntxAreaUtils::area_on_sphere(), moab::ParallelComm::check_all_shared_handles(), compute_lagrange_mesh_on_sphere(), moab::Intx2Mesh::create_departure_mesh_2nd_alg(), moab::Interface::create_meshset(), EPS1, ErrorCode, example, moab::Intx2Mesh::FindMaxEdges(), moab::ParallelComm::get_pcomm(), input_mesh_file(), moab::Intx2Mesh::intersect_meshes(), moab::Interface::load_file(), mb, MB_CHK_ERR_RET, MESHSET_SET, moab::ParallelComm::proc_config(), moab::ProcConfig::proc_rank(), Radius, rank, moab::Intx2Mesh::set_box_error(), moab::Intx2Mesh::set_error_tolerance(), moab::Intx2MeshOnSphere::set_radius_destination_mesh(), moab::Intx2MeshOnSphere::set_radius_source_mesh(), and moab::Interface::write_file().
Referenced by main().
{
std::string opts = std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" );
Core moab;
Interface& mb = moab;
EntityHandle euler_set;
ErrorCode rval;
rval = mb.create_meshset( MESHSET_SET, euler_set );MB_CHK_ERR_RET( rval );
std::string example( TestDir + "unittest/" + input_mesh_file );
rval = mb.load_file( example.c_str(), &euler_set, opts.c_str() );
ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );MB_CHK_ERR_RET( rval );
rval = pcomm->check_all_shared_handles();MB_CHK_ERR_RET( rval );
// everybody will get a DP tag, including the non owned entities; so exchange tags is not
// required for LOC (here)
rval = compute_lagrange_mesh_on_sphere( &mb, euler_set );MB_CHK_ERR_RET( rval );
int rank = pcomm->proc_config().proc_rank();
std::stringstream ste;
ste << "initial" << rank << ".vtk";
mb.write_file( ste.str().c_str(), 0, 0, &euler_set, 1 );
Intx2MeshOnSphere worker( &mb );
worker.set_radius_source_mesh( Radius );
worker.set_radius_destination_mesh( Radius );
worker.set_box_error( EPS1 ); //
// worker.SetEntityType(MBQUAD);
worker.set_error_tolerance( Radius * 1.e-8 );
std::cout << "error tolerance epsilon_1=" << Radius * 1.e-8 << "\n";
// worker.locate_departure_points(euler_set);
rval = worker.FindMaxEdges( euler_set, euler_set ); // departure will be the same max_edges
// we need to make sure the covering set is bigger than the euler mesh
EntityHandle covering_lagr_set;
rval = mb.create_meshset( MESHSET_SET, covering_lagr_set );MB_CHK_ERR_RET( rval );
rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );MB_CHK_ERR_RET( rval );
std::stringstream ss;
ss << "partial" << rank << ".vtk";
mb.write_file( ss.str().c_str(), 0, 0, &covering_lagr_set, 1 );
EntityHandle outputSet;
rval = mb.create_meshset( MESHSET_SET, outputSet );MB_CHK_ERR_RET( rval );
rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );MB_CHK_ERR_RET( rval );
// std::string opts_write("PARALLEL=WRITE_PART");
// rval = mb.write_file("manuf.h5m", 0, opts_write.c_str(), &outputSet, 1);
// std::string opts_write("");
std::stringstream outf;
outf << "intersect" << rank << ".h5m";
rval = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );MB_CHK_ERR_RET( rval );
moab::IntxAreaUtils sphAreaUtils;
double intx_area = sphAreaUtils.area_on_sphere( &mb, outputSet, Radius );
double arrival_area = sphAreaUtils.area_on_sphere( &mb, euler_set, Radius );
std::cout << "On rank : " << rank << " arrival area: " << arrival_area << " intersection area:" << intx_area
<< " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";
}
| double deltaT = 1.e-6 |
Definition at line 44 of file cslam_par_test.cpp.
Referenced by compute_lagrange_mesh_on_sphere(), and main().
| double EPS1 = 0.2 |
Definition at line 41 of file cslam_par_test.cpp.
Referenced by main(), test_eval(), test_evals(), test_hex_nat_coords(), test_intx_in_parallel_elem_based(), and test_intx_mpas().
| double Radius = 1.0 |
Definition at line 43 of file cslam_par_test.cpp.
Referenced by moab::IntxAreaUtils::area_spherical_triangle_lHuiller(), compute_lagrange_mesh_on_sphere(), main(), moab::IntxAreaUtils::spherical_angle(), and test_intx_in_parallel_elem_based().
1.7.6.1