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