MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#include <iostream>
#include <cmath>
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
#include "moab/ProgOptions.hpp"
#include "MBTagConventions.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
#include "IntxUtilsCSLAM.hpp"
#include "TestUtil.hpp"
Go to the source code of this file.
Functions | |
ErrorCode | manufacture_lagrange_mesh_on_sphere (Interface *mb, EntityHandle euler_set, EntityHandle &lagr_set) |
int | main (int argc, char **argv) |
Variables | |
double | gtol = 0.0001 |
double | CubeSide = 6. |
double | t = 0.1 |
double | delta_t = 0.43 |
int main | ( | int | argc, |
char ** | argv | ||
) |
Definition at line 97 of file case1_test.cpp.
References moab::IntxAreaUtils::area_on_sphere(), moab::Interface::create_meshset(), CubeSide, delta_t, moab::IntxUtils::enforce_convexity(), ErrorCode, moab::Intx2Mesh::FindMaxEdges(), gtol, moab::Intx2Mesh::intersect_meshes(), moab::Interface::load_file(), manufacture_lagrange_mesh_on_sphere(), mb, MB_SUCCESS, MESHSET_SET, radius, moab::Intx2Mesh::set_error_tolerance(), moab::Intx2MeshOnSphere::set_radius_destination_mesh(), moab::Intx2MeshOnSphere::set_radius_source_mesh(), STRINGIFY, and moab::Interface::write_file().
{ const char* filename_mesh1 = STRINGIFY( MESHDIR ) "/mbcslam/eulerHomme.vtk"; if( argc > 1 ) { int index = 1; while( index < argc ) { if( !strcmp( argv[index], "-gtol" ) ) // this is for geometry tolerance { gtol = atof( argv[++index] ); } if( !strcmp( argv[index], "-cube" ) ) { CubeSide = atof( argv[++index] ); } if( !strcmp( argv[index], "-dt" ) ) { delta_t = atof( argv[++index] ); } if( !strcmp( argv[index], "-input" ) ) { filename_mesh1 = argv[++index]; } index++; } } std::cout << " case 1: use -gtol " << gtol << " -dt " << delta_t << " -cube " << CubeSide << " -input " << filename_mesh1 << "\n"; Core moab; Interface& mb = moab; EntityHandle euler_set; ErrorCode rval; rval = mb.create_meshset( MESHSET_SET, euler_set ); if( MB_SUCCESS != rval ) return 1; rval = mb.load_file( filename_mesh1, &euler_set ); if( MB_SUCCESS != rval ) return 1; // everybody will get a DP tag, including the non owned entities; so exchange tags is not // required for LOC (here) EntityHandle lagrange_set; rval = manufacture_lagrange_mesh_on_sphere( &mb, euler_set, lagrange_set ); if( MB_SUCCESS != rval ) return 1; rval = mb.write_file( "lagrIni.h5m", 0, 0, &lagrange_set, 1 ); if( MB_SUCCESS != rval ) std::cout << "can't write lagr set\n"; rval = moab::IntxUtils::enforce_convexity( &mb, lagrange_set ); if( MB_SUCCESS != rval ) return 1; rval = mb.write_file( "lagr.h5m", 0, 0, &lagrange_set, 1 ); if( MB_SUCCESS != rval ) std::cout << "can't write lagr set\n"; Intx2MeshOnSphere worker( &mb ); double radius = CubeSide / 2 * sqrt( 3. ); // input worker.set_radius_source_mesh( radius ); worker.set_radius_destination_mesh( radius ); // worker.SetEntityType(MBQUAD); worker.set_error_tolerance( gtol ); std::cout << "error tolerance epsilon_1=" << gtol << "\n"; EntityHandle outputSet; rval = mb.create_meshset( MESHSET_SET, outputSet ); if( MB_SUCCESS != rval ) return 1; rval = worker.FindMaxEdges( lagrange_set, euler_set ); if( MB_SUCCESS != rval ) return 1; rval = worker.intersect_meshes( lagrange_set, euler_set, outputSet ); if( MB_SUCCESS != rval ) return 1; // std::string opts_write(""); std::stringstream outf; outf << "intersect1" << ".h5m"; rval = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 ); if( MB_SUCCESS != rval ) std::cout << "can't write output\n"; 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 << " Arrival area: " << arrival_area << " intersection area:" << intx_area << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n"; if( MB_SUCCESS != rval ) return 1; return 0; }
ErrorCode manufacture_lagrange_mesh_on_sphere | ( | Interface * | mb, |
EntityHandle | euler_set, | ||
EntityHandle & | lagr_set | ||
) |
Definition at line 29 of file case1_test.cpp.
References moab::Interface::add_entities(), moab::Range::begin(), moab::Interface::create_element(), moab::Interface::create_meshset(), moab::Interface::create_vertex(), CubeSide, delta_t, IntxUtilsCSLAM::departure_point_case1(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_type(), MB_SUCCESS, MBQUAD, MESHSET_SET, radius, and t.
Referenced by main(), and test_intx_mpas().
{ /* * get all quads first, then vertices, then move them on the surface of the sphere * radius is in, it comes from MeshKit/python/examples/manufHomme.py : * length = 6. * each edge of the cube will be divided using this meshcount * meshcount = 11 * circumscribed sphere radius * radius = length * math.sqrt(3) /2 */ double radius = CubeSide / 2 * sqrt( 3. ); // our value depends on cube side Range quads; ErrorCode rval = mb->get_entities_by_type( euler_set, MBQUAD, quads ); if( MB_SUCCESS != rval ) return rval; Range connecVerts; rval = mb->get_connectivity( quads, connecVerts ); if( MB_SUCCESS != rval ) return rval; // create new set rval = mb->create_meshset( MESHSET_SET, lagr_set ); if( MB_SUCCESS != rval ) return rval; // get the coordinates of the old mesh, and move it around the sphere according to case 1 // now put the vertices in the right place.... // int vix=0; // vertex index in new array // first create departure points (vertices in the lagrange mesh) // then connect them in quads std::map< EntityHandle, EntityHandle > newNodes; for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit ) { EntityHandle oldV = *vit; CartVect posi; rval = mb->get_coords( &oldV, 1, &( posi[0] ) ); if( MB_SUCCESS != rval ) return rval; // Intx utils, case 1 CartVect newPos; IntxUtilsCSLAM::departure_point_case1( posi, t, delta_t, newPos ); newPos = radius * newPos; EntityHandle new_vert; rval = mb->create_vertex( &( newPos[0] ), new_vert ); if( MB_SUCCESS != rval ) return rval; newNodes[oldV] = new_vert; } for( Range::iterator it = quads.begin(); it != quads.end(); ++it ) { EntityHandle q = *it; int nnodes; const EntityHandle* conn4; rval = mb->get_connectivity( q, conn4, nnodes ); if( MB_SUCCESS != rval ) return rval; EntityHandle new_conn[4]; for( int i = 0; i < nnodes; i++ ) { EntityHandle v1 = conn4[i]; new_conn[i] = newNodes[v1]; } EntityHandle new_quad; rval = mb->create_element( MBQUAD, new_conn, 4, new_quad ); if( MB_SUCCESS != rval ) return rval; rval = mb->add_entities( lagr_set, &new_quad, 1 ); if( MB_SUCCESS != rval ) return rval; } return rval; }
double CubeSide = 6. |
Definition at line 26 of file case1_test.cpp.
Referenced by main(), and manufacture_lagrange_mesh_on_sphere().
double delta_t = 0.43 |
Definition at line 27 of file case1_test.cpp.
Referenced by compute_tracer_case1(), get_departure_grid(), main(), manufacture_lagrange_mesh_on_sphere(), and moab::SmoothCurve::u_from_position().
double gtol = 0.0001 |
Definition at line 25 of file case1_test.cpp.
double t = 0.1 |
Definition at line 27 of file case1_test.cpp.
Referenced by add_tag_counts(), MetisPartitioner::assemble_taggedsets_graph(), moab::ReadNASTRAN::assign_ids(), moab::HiReconstruction::average_vertex_tangent(), moab::box_from_axes(), check_identical_mesh(), check_sets_sizes(), moab::SequenceManager::clear(), moab::GeomUtil::closest_location_on_polygon(), moab::GeomUtil::closest_location_on_tri(), moab::NCWriteHelper::collect_variable_data(), compute_tracer_case1(), compute_velocity_case1(), count_owned_entities(), moab::ReadCGM::create_group_entsets(), moab::Coupler::create_tuples(), moab::BSPTreePoly::cut_polyhedron(), moab::ReadABAQUS::cyl2rect(), moab::Skinner::deinitialize(), departure_point_swirl_rot(), moab::HalfFacetRep::determine_border_vertices(), moab::CN::Dimension(), do_test(), moab::ParallelComm::estimate_ents_buffer_size(), moab::DenseTag::find_entities_with_value(), moab::VarLenDenseTag::find_entities_with_value(), moab::Skinner::find_skin_vertices_3D(), smoab::Interface::findEntitiesWithTag(), gather_tag_counts(), get_by_all_types_and_tag(), get_departure_grid(), moab::SequenceManager::get_entities(), moab::BitTag::get_entities_with_bits(), moab::Coupler::get_matching_entities(), moab::AEntityFactory::get_memory_use(), moab::BitTag::get_memory_use(), moab::DenseTag::get_memory_use(), moab::VarLenDenseTag::get_memory_use(), moab::SequenceManager::get_number_entities(), moab::ReadUtil::get_ordered_vertices(), moab::BitTag::get_tagged(), moab::get_tagged(), moab::DenseTag::get_tagged_entities(), moab::ReorderTool::handle_order_from_int_tag(), iMeshP_pushTags(), moab::NCHelperHOMME::init_mesh_vals(), moab::NCHelperMPAS::init_mesh_vals(), moab::NCHelperFV::init_mesh_vals(), moab::NCHelperGCRM::init_mesh_vals(), moab::NCHelperEuler::init_mesh_vals(), moab::Coupler::interpolate(), MBiMesh::is_ent_handle_tag(), MBiMesh::is_set_handle_tag(), itaps_tag_cast(), moab::RayIntersectSets::leaf(), moab::TempestOnlineMap::LinearRemapGLLtoGLL2_MOAB(), moab::TempestOnlineMap::LinearRemapGLLtoGLL2_Pointwise_MOAB(), moab::Core::load_file(), lob_bnd_2(), main(), make_basic_sequence(), manufacture_lagrange_mesh_on_sphere(), min_edge_length(), moab::Coupler::normalize_subset(), MBiMesh::note_ent_handle_tag(), MBiMesh::note_set_handle_tag(), MBiMesh::note_tag_destroyed(), moab::operator<<(), MetisPartitioner::partition_mesh(), ZoltanPartitioner::partition_mesh_and_geometry(), moab::ParallelMergeMesh::PerformRealSort(), permutation(), moab::permute_this(), moab::plane_cut_edge(), moab::point_perp(), moab::SmoothCurve::position_from_u(), moab::Core::print_database(), print_memory_stats(), print_partitioned_entities(), print_tag_counts(), moab::DGMSolver::qr_polyfit_safeguarded(), moab::TupleList::radix_index_sort(), moab::TupleList::radix_offsets(), moab::AdaptiveKDTree::ray_intersect_triangles(), moab::OrientedBoxTreeTool::ray_intersect_triangles(), moab::SmoothFace::ray_intersection_correct(), moab::GeomUtil::ray_tri_intersect(), moab::ReadHDF5::read_elems(), moab::ScdNCHelper::read_scd_variables_to_nonset(), moab::ScdNCHelper::read_scd_variables_to_nonset_allocate(), moab::NCHelperHOMME::read_ucd_variables_to_nonset(), moab::NCHelperMPAS::read_ucd_variables_to_nonset(), moab::NCHelperGCRM::read_ucd_variables_to_nonset(), moab::NCHelperHOMME::read_ucd_variables_to_nonset_allocate(), moab::NCHelperMPAS::read_ucd_variables_to_nonset_allocate(), moab::NCHelperGCRM::read_ucd_variables_to_nonset_allocate(), moab::NCHelper::read_variables_to_set(), moab::NCHelper::read_variables_to_set_allocate(), moab::SimplexTemplateRefiner::refine_3_simplex(), moab::BitTag::release_all_data(), moab::SequenceManager::release_tag_array(), moab::ReorderTool::reorder_entities(), ZoltanPartitioner::repartition(), moab::CN::resetPermutation(), moab::ParallelComm::resolve_shared_sets(), moab::rev_permute_this(), moab::CN::setPermutation(), moab::TagInfo::size_from_data_type(), skin_common(), moab::ReadABAQUS::sph2rect(), moab::OrientedBoxTreeTool::sphere_intersect_triangles(), moab::ParallelMergeMesh::SwapTuples(), tag_time(), test_create_tag(), test_create_var_len_tag(), test_ho_node_index(), test_ho_node_parent(), test_plucker_ray_tri_intersect(), test_ray_tri_intersect(), test_read_nothing_common(), test_type_names(), test_var_length_handle_tag(), test_write_elements(), test_write_read_many_tags(), moab::HomXform::three_pt_xform(), TIME(), TIME_DEL(), TIME_QRY(), moab::Core::valid_tag_handle(), moab::WriteVtk::write_bit_tag(), moab::WriteVtk::write_elems(), moab::NCWriteHOMME::write_nonset_variables(), moab::NCWriteGCRM::write_nonset_variables(), moab::NCWriteMPAS::write_nonset_variables(), moab::ScdNCWriteHelper::write_nonset_variables(), moab::WriteHDF5::write_qa(), moab::WriteVtk::write_tag(), and moab::TempestOnlineMap::WriteHDF5MapFile().