MOAB: Mesh Oriented datABase  (version 5.2.1)
par_intx_sph.cpp
Go to the documentation of this file.
00001 /*
00002  * par_intx_sph.cpp
00003  *  test to trigger intersection on a sphere in parallel
00004  *  it will start from an eulerian mesh
00005  *  the mesh is read in parallel; lagrangian mesh is manufactured on the fly (part of the test), in
00006  *    a different set.
00007  *
00008  *  intersections will be performed on the euler mesh
00009  *  Created on: Nov 14, 2012
00010  */
00011 
00012 #include <iostream>
00013 #include <sstream>
00014 #include <time.h>
00015 #include <stdlib.h>
00016 #include <stdio.h>
00017 #include <string.h>
00018 #include "moab/Core.hpp"
00019 #include "moab/Interface.hpp"
00020 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
00021 #include <math.h>
00022 #include "TestUtil.hpp"
00023 #include "moab/ParallelComm.hpp"
00024 #include "moab/ProgOptions.hpp"
00025 #include "MBParallelConventions.h"
00026 #include "moab/ReadUtilIface.hpp"
00027 #include "MBTagConventions.hpp"
00028 
00029 #include "moab/IntxMesh/IntxUtils.hpp"
00030 
00031 // for M_PI
00032 #include <math.h>
00033 
00034 using namespace moab;
00035 // some input data
00036 double EPS1 = 0.2;                               // this is for box error
00037 std::string input_mesh_file( "Homme_2pt.h5m" );  // input file, partitioned correctly
00038 std::string mpas_file( "mpas_p8.h5m" );
00039 double CubeSide = 6.;  // the above file starts with cube side 6; radius depends on cube side
00040 double radius;
00041 void test_intx_in_parallel_elem_based();
00042 void test_intx_mpas();
00043 
00044 int main( int argc, char** argv )
00045 {
00046     MPI_Init( &argc, &argv );
00047     EPS1       = 0.2;
00048     int result = 0;
00049 
00050     if( argc > 1 )
00051     {
00052         int index = 1;
00053         while( index < argc )
00054         {
00055             if( !strcmp( argv[index], "-eps" ) )  // this is for box error
00056             { EPS1 = atof( argv[++index] ); }
00057             if( !strcmp( argv[index], "-input" ) ) { input_mesh_file = argv[++index]; }
00058             if( !strcmp( argv[index], "-cube" ) ) { CubeSide = atof( argv[++index] ); }
00059             index++;
00060         }
00061     }
00062 
00063     radius = CubeSide / 2 * sqrt( 3. );
00064     result += RUN_TEST( test_intx_in_parallel_elem_based );
00065     radius = 1.;
00066     result += RUN_TEST( test_intx_mpas );
00067 
00068     MPI_Finalize();
00069     return result;
00070 }
00071 // will save the LOC tag on the euler nodes
00072 ErrorCode manufacture_lagrange_mesh_on_sphere( Interface* mb, EntityHandle euler_set )
00073 {
00074     /*
00075      * get all quads first, then vertices, then move them on the surface of the sphere
00076      *  radius is in, it comes from MeshKit/python/examples/manufHomme.py :
00077      *  length = 6.
00078      *  each edge of the cube will be divided using this meshcount
00079      *  meshcount = 11
00080      *   circumscribed sphere radius
00081      *   radius = length * math.sqrt(3) /2
00082      */
00083     // radius = CubeSide/2*sqrt(3.);// our value depends on cube side
00084     Range quads;
00085     ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, quads );CHECK_ERR( rval );
00086 
00087     Range connecVerts;
00088     rval = mb->get_connectivity( quads, connecVerts );CHECK_ERR( rval );
00089 
00090     // the LOC tag, should be provided by the user?
00091     Tag tagh = 0;
00092     std::string tag_name( "DP" );
00093     rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval );
00094     void* data;  // pointer to the LOC in memory, for each vertex
00095     int count;
00096 
00097     rval = mb->tag_iterate( tagh, connecVerts.begin(), connecVerts.end(), count, data );CHECK_ERR( rval );
00098     // here we are checking contiguity
00099     assert( count == (int)connecVerts.size() );
00100     double* ptr_DP = (double*)data;
00101     // get the coordinates of the old mesh, and move it around the sphere in the same way as in the
00102     // python script
00103 
00104     // now put the vertices in the right place....
00105     // int vix=0; // vertex index in new array
00106     double t = 0.1, T = 5;  // check the script
00107     double time = 0.05;
00108     double rot  = M_PI / 10;
00109     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00110     {
00111         EntityHandle oldV = *vit;
00112         CartVect posi;
00113         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
00114         // do some mumbo jumbo, as in python script
00115         moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
00116         double lat1                            = sphCoord.lat - 2 * M_PI * t / T;  // 0.1/5
00117         double uu = 3 * radius / T * pow( sin( lat1 ), 2 ) * sin( 2 * sphCoord.lon ) * cos( M_PI * t / T );
00118         uu += 2 * M_PI * cos( sphCoord.lon ) / T;
00119         double vv = 3 * radius / T * ( sin( 2 * lat1 ) ) * cos( sphCoord.lon ) * cos( M_PI * t / T );
00120         double vx = -uu * sin( sphCoord.lon ) - vv * sin( sphCoord.lat ) * cos( sphCoord.lon );
00121         double vy = -uu * cos( sphCoord.lon ) - vv * sin( sphCoord.lat ) * sin( sphCoord.lon );
00122         double vz = vv * cos( sphCoord.lat );
00123         posi      = posi + time * CartVect( vx, vy, vz );
00124         double x2 = posi[0] * cos( rot ) - posi[1] * sin( rot );
00125         double y2 = posi[0] * sin( rot ) + posi[1] * cos( rot );
00126         CartVect newPos( x2, y2, posi[2] );
00127         double len1 = newPos.length();
00128         newPos      = radius * newPos / len1;
00129 
00130         ptr_DP[0] = newPos[0];
00131         ptr_DP[1] = newPos[1];
00132         ptr_DP[2] = newPos[2];
00133         ptr_DP += 3;  // increment to the next node
00134     }
00135 
00136     return rval;
00137 }
00138 
00139 void test_intx_in_parallel_elem_based()
00140 {
00141     std::string opts = std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
00142                        std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" );
00143     Core moab;
00144     Interface& mb = moab;
00145     EntityHandle euler_set;
00146     ErrorCode rval;
00147     rval = mb.create_meshset( MESHSET_SET, euler_set );CHECK_ERR( rval );
00148     std::string example( TestDir + "/" + input_mesh_file );
00149 
00150     rval = mb.load_file( example.c_str(), &euler_set, opts.c_str() );CHECK_ERR( rval );
00151 
00152     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00153     CHECK( NULL != pcomm );
00154 
00155     rval = pcomm->check_all_shared_handles();CHECK_ERR( rval );
00156 
00157     // everybody will get a DP tag, including the non owned entities; so exchange tags is not
00158     // required for LOC (here)
00159     rval = manufacture_lagrange_mesh_on_sphere( &mb, euler_set );CHECK_ERR( rval );
00160 
00161     int rank = pcomm->proc_config().proc_rank();
00162 
00163     std::stringstream ste;
00164     ste << "initial" << rank << ".vtk";
00165     mb.write_file( ste.str().c_str(), 0, 0, &euler_set, 1 );
00166 
00167     Intx2MeshOnSphere worker( &mb );
00168 
00169     // double radius= CubeSide/2 * sqrt(3.) ; // input
00170     worker.set_radius_source_mesh( radius );
00171     worker.set_radius_destination_mesh( radius );
00172     worker.set_box_error( EPS1 );  //
00173     // worker.SetEntityType(MBQUAD);
00174 
00175     worker.set_error_tolerance( radius * 1.e-8 );
00176     std::cout << "error tolerance epsilon_1=" << radius * 1.e-8 << "\n";
00177     //  worker.locate_departure_points(euler_set);
00178     // set pcomm
00179     worker.set_parallel_comm( pcomm );
00180 
00181     rval = worker.FindMaxEdges( euler_set, euler_set );  // use euler set for both, it is just finding max_edges_1 and 2
00182     // we need to make sure the covering set is bigger than the euler mesh
00183     EntityHandle covering_lagr_set;
00184     rval = mb.create_meshset( MESHSET_SET, covering_lagr_set );CHECK_ERR( rval );
00185 
00186     rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );CHECK_ERR( rval );
00187 
00188     std::stringstream ss;
00189     ss << "partial" << rank << ".vtk";
00190     mb.write_file( ss.str().c_str(), 0, 0, &covering_lagr_set, 1 );
00191     EntityHandle outputSet;
00192     rval = mb.create_meshset( MESHSET_SET, outputSet );CHECK_ERR( rval );
00193     rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );CHECK_ERR( rval );
00194 
00195     IntxAreaUtils areaAdaptor;
00196     // std::string opts_write("PARALLEL=WRITE_PART");
00197     // rval = mb.write_file("manuf.h5m", 0, opts_write.c_str(), &outputSet, 1);
00198     // std::string opts_write("");
00199     std::stringstream outf;
00200     outf << "intersect" << rank << ".h5m";
00201     rval                = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );
00202     double intx_area    = areaAdaptor.area_on_sphere( &mb, outputSet, radius );
00203     double arrival_area = areaAdaptor.area_on_sphere( &mb, euler_set, radius );
00204     std::cout << "On rank : " << rank << " arrival area: " << arrival_area << "  intersection area:" << intx_area
00205               << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";CHECK_ERR( rval );
00206 }
00207 
00208 void test_intx_mpas()
00209 {
00210     std::string opts = std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
00211                        std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" );
00212     Core moab;
00213     Interface& mb = moab;
00214     EntityHandle euler_set;
00215     ErrorCode rval;
00216     rval = mb.create_meshset( MESHSET_SET, euler_set );CHECK_ERR( rval );
00217     std::string example( TestDir + "/" + mpas_file );
00218 
00219     rval = mb.load_file( example.c_str(), &euler_set, opts.c_str() );
00220 
00221     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00222 
00223     rval = pcomm->check_all_shared_handles();CHECK_ERR( rval );
00224 
00225     // everybody will get a DP tag, including the non owned entities; so exchange tags is not
00226     // required for LOC (here)
00227     rval = manufacture_lagrange_mesh_on_sphere( &mb, euler_set );CHECK_ERR( rval );
00228 
00229     int rank = pcomm->proc_config().proc_rank();
00230 
00231     std::stringstream ste;
00232     ste << "initial" << rank << ".vtk";
00233     mb.write_file( ste.str().c_str(), 0, 0, &euler_set, 1 );
00234 
00235     Intx2MeshOnSphere worker( &mb );
00236 
00237     // double radius= CubeSide/2 * sqrt(3.) ; // input
00238     worker.set_radius_source_mesh( radius );
00239     worker.set_radius_destination_mesh( radius );
00240     worker.set_box_error( EPS1 );  //
00241     // worker.SetEntityType(MBQUAD);
00242 
00243     worker.set_error_tolerance( radius * 1.e-8 );
00244     std::cout << "error tolerance epsilon_1=" << radius * 1.e-8 << "\n";
00245     //  worker.locate_departure_points(euler_set);
00246 
00247     rval = worker.FindMaxEdges( euler_set, euler_set );
00248 
00249     // we need to make sure the covering set is bigger than the euler mesh
00250     EntityHandle covering_lagr_set;
00251     rval = mb.create_meshset( MESHSET_SET, covering_lagr_set );CHECK_ERR( rval );
00252 
00253     // set pcomm
00254     worker.set_parallel_comm( pcomm );
00255 
00256     rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );CHECK_ERR( rval );
00257 
00258     std::stringstream ss;
00259     ss << "partial" << rank << ".vtk";
00260     mb.write_file( ss.str().c_str(), 0, 0, &covering_lagr_set, 1 );
00261     rval = moab::IntxUtils::enforce_convexity( &mb, covering_lagr_set, rank );CHECK_ERR( rval );
00262     std::stringstream ss2;
00263     ss2 << "partialConvex" << rank << ".vtk";
00264     mb.write_file( ss2.str().c_str(), 0, 0, &covering_lagr_set, 1 );
00265     EntityHandle outputSet;
00266     rval = mb.create_meshset( MESHSET_SET, outputSet );CHECK_ERR( rval );
00267     rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );CHECK_ERR( rval );
00268 
00269     // std::string opts_write("PARALLEL=WRITE_PART");
00270     // rval = mb.write_file("manuf.h5m", 0, opts_write.c_str(), &outputSet, 1);
00271     // std::string opts_write("");
00272     std::stringstream outf;
00273     outf << "intersect" << rank << ".h5m";
00274     rval = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );
00275 
00276     IntxAreaUtils areaAdaptor;
00277     double intx_area    = areaAdaptor.area_on_sphere( &mb, outputSet, radius );
00278     double arrival_area = areaAdaptor.area_on_sphere( &mb, euler_set, radius );
00279     std::cout << "On rank : " << rank << " arrival area: " << arrival_area << "  intersection area:" << intx_area
00280               << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";CHECK_ERR( rval );
00281 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines