MOAB: Mesh Oriented datABase  (version 5.4.0)
cslam_par_test.cpp
Go to the documentation of this file.
00001 /*
00002  * cslam_par_test.cpp
00003  *  test to trigger intersection on a sphere in parallel
00004  *  it will start from an eulerian mesh (mesh + velocity from Mark Taylor)
00005  *   file: VELO00.h5m; mesh + velo at 850 milibar, in an unstructured mesh refined from
00006  *   a cube sphere grid
00007  *  the mesh is read in parallel (euler mesh);
00008  *
00009  *   lagrangian mesh is obtained using
00010  *   pos (t-dt) = pos(t) -Velo(t)*dt
00011  *   then, project on the sphere with a given radius
00012  *
00013  *  Created on: Apr 22, 2013
00014  */
00015 
00016 #include <iostream>
00017 #include <sstream>
00018 #include <ctime>
00019 #include <cstdlib>
00020 #include <cstdio>
00021 #include <cstring>
00022 #include "moab/Core.hpp"
00023 #include "moab/Interface.hpp"
00024 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
00025 #include <cmath>
00026 #include "TestUtil.hpp"
00027 #include "moab/ParallelComm.hpp"
00028 #include "moab/ProgOptions.hpp"
00029 #include "MBParallelConventions.h"
00030 #include "moab/ReadUtilIface.hpp"
00031 #include "MBTagConventions.hpp"
00032 
00033 #include "moab/IntxMesh/IntxUtils.hpp"
00034 #include "IntxUtilsCSLAM.hpp"
00035 
00036 // for M_PI
00037 #include <cmath>
00038 
00039 using namespace moab;
00040 // some input data
00041 double EPS1 = 0.2;                                // this is for box error
00042 std::string input_mesh_file( "VELO00_16p.h5m" );  // input file, partitioned correctly
00043 double Radius = 1.0;                              // change to radius
00044 double deltaT = 1.e-6;
00045 void test_intx_in_parallel_elem_based();
00046 
00047 int main( int argc, char** argv )
00048 {
00049     MPI_Init( &argc, &argv );
00050     EPS1       = 0.000002;
00051     int result = 0;
00052 
00053     if( argc > 1 )
00054     {
00055         int index = 1;
00056         while( index < argc )
00057         {
00058             if( !strcmp( argv[index], "-eps" ) )  // this is for box error
00059             {
00060                 EPS1 = atof( argv[++index] );
00061             }
00062             if( !strcmp( argv[index], "-input" ) )
00063             {
00064                 input_mesh_file = argv[++index];
00065             }
00066             if( !strcmp( argv[index], "-radius" ) )
00067             {
00068                 Radius = atof( argv[++index] );
00069             }
00070             if( !strcmp( argv[index], "-deltaT" ) )
00071             {
00072                 deltaT = atof( argv[++index] );
00073             }
00074             index++;
00075         }
00076     }
00077     std::cout << " run: -input " << input_mesh_file << "  -eps " << EPS1 << " -radius " << Radius << " -deltaT "
00078               << deltaT << "\n";
00079 
00080     result += RUN_TEST( test_intx_in_parallel_elem_based );
00081 
00082     MPI_Finalize();
00083     return result;
00084 }
00085 // will save the LOC tag on the euler nodes
00086 ErrorCode compute_lagrange_mesh_on_sphere( Interface* mb, EntityHandle euler_set )
00087 {
00088     /*
00089      * get all quads first, then vertices, then move them on the surface of the sphere
00090      *  radius is 1, usually
00091      *  pos (t-dt) = pos(t) -Velo(t)*dt; this will be lagrange mesh, on each processor
00092      */
00093     Range quads;
00094     ErrorCode rval = mb->get_entities_by_type( euler_set, MBQUAD, quads );MB_CHK_ERR( rval );
00095 
00096     Range connecVerts;
00097     rval = mb->get_connectivity( quads, connecVerts );MB_CHK_ERR( rval );
00098 
00099     // the LOC tag, should be provided by the user?
00100     Tag tagh = 0;
00101     std::string tag_name( "DP" );
00102     rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00103     void* data;  // pointer to the DP in memory, for each vertex
00104     int count;
00105 
00106     rval = mb->tag_iterate( tagh, connecVerts.begin(), connecVerts.end(), count, data );MB_CHK_ERR( rval );
00107     // here we are checking contiguity
00108     assert( count == (int)connecVerts.size() );
00109     double* ptr_DP = (double*)data;
00110     // get the coordinates of the old mesh, and move it around using velocity tag
00111 
00112     Tag tagv = 0;
00113     std::string velo_tag_name( "VELO" );
00114     rval = mb->tag_get_handle( velo_tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagv, MB_TAG_DENSE );MB_CHK_ERR( rval );
00115 
00116     /*void *datavelo; // pointer to the VELO in memory, for each vertex
00117 
00118     rval = mb->tag_iterate(tagv, connecVerts.begin(), connecVerts.end(), count, datavelo);MB_CHK_ERR(rval);*/
00119     // here we are checking contiguity
00120     assert( count == (int)connecVerts.size() );
00121     // now put the vertices in the right place....
00122     // int vix=0; // vertex index in new array
00123 
00124     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00125     {
00126         EntityHandle oldV = *vit;
00127         CartVect posi;
00128         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
00129         CartVect velo;
00130         rval = mb->tag_get_data( tagv, &oldV, 1, (void*)&( velo[0] ) );MB_CHK_ERR( rval );
00131         // do some mumbo jumbo, as in python script
00132         CartVect newPos = posi - deltaT * velo;
00133         double len1     = newPos.length();
00134         newPos          = Radius * newPos / len1;
00135 
00136         ptr_DP[0] = newPos[0];
00137         ptr_DP[1] = newPos[1];
00138         ptr_DP[2] = newPos[2];
00139         ptr_DP += 3;  // increment to the next node
00140     }
00141 
00142     return rval;
00143 }
00144 
00145 void test_intx_in_parallel_elem_based()
00146 {
00147     std::string opts = std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
00148                        std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" );
00149     Core moab;
00150     Interface& mb = moab;
00151     EntityHandle euler_set;
00152     ErrorCode rval;
00153     rval = mb.create_meshset( MESHSET_SET, euler_set );MB_CHK_ERR_RET( rval );
00154     std::string example( TestDir + "unittest/" + input_mesh_file );
00155 
00156     rval = mb.load_file( example.c_str(), &euler_set, opts.c_str() );
00157 
00158     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );MB_CHK_ERR_RET( rval );
00159 
00160     rval = pcomm->check_all_shared_handles();MB_CHK_ERR_RET( rval );
00161 
00162     // everybody will get a DP tag, including the non owned entities; so exchange tags is not
00163     // required for LOC (here)
00164     rval = compute_lagrange_mesh_on_sphere( &mb, euler_set );MB_CHK_ERR_RET( rval );
00165 
00166     int rank = pcomm->proc_config().proc_rank();
00167 
00168     std::stringstream ste;
00169     ste << "initial" << rank << ".vtk";
00170     mb.write_file( ste.str().c_str(), 0, 0, &euler_set, 1 );
00171 
00172     Intx2MeshOnSphere worker( &mb );
00173 
00174     worker.set_radius_source_mesh( Radius );
00175     worker.set_radius_destination_mesh( Radius );
00176     worker.set_box_error( EPS1 );  //
00177     // worker.SetEntityType(MBQUAD);
00178 
00179     worker.set_error_tolerance( Radius * 1.e-8 );
00180     std::cout << "error tolerance epsilon_1=" << Radius * 1.e-8 << "\n";
00181     //  worker.locate_departure_points(euler_set);
00182 
00183     rval = worker.FindMaxEdges( euler_set, euler_set );  // departure will be the same max_edges
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 );MB_CHK_ERR_RET( rval );
00187 
00188     rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );MB_CHK_ERR_RET( 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 );MB_CHK_ERR_RET( rval );
00195     rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );MB_CHK_ERR_RET( rval );
00196 
00197     // std::string opts_write("PARALLEL=WRITE_PART");
00198     // rval = mb.write_file("manuf.h5m", 0, opts_write.c_str(), &outputSet, 1);
00199     // std::string opts_write("");
00200     std::stringstream outf;
00201     outf << "intersect" << rank << ".h5m";
00202     rval = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );MB_CHK_ERR_RET( rval );
00203 
00204     moab::IntxAreaUtils sphAreaUtils;
00205     double intx_area    = sphAreaUtils.area_on_sphere( &mb, outputSet, Radius );
00206     double arrival_area = sphAreaUtils.area_on_sphere( &mb, euler_set, Radius );
00207     std::cout << "On rank : " << rank << " arrival area: " << arrival_area << "  intersection area:" << intx_area
00208               << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";
00209 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines