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