Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
intx_mpas.cpp
Go to the documentation of this file.
00001 /*
00002  * mpas file test
00003  *
00004  *  Created on: Feb 12, 2013
00005  */
00006 
00007 // copy from case1 test
00008 #include "moab/MOABConfig.h"
00009 
00010 #ifdef MOAB_HAVE_ZOLTAN
00011 
00012 #include <iostream>
00013 #include <sstream>
00014 #include <ctime>
00015 #include <cstdlib>
00016 #include <cstdio>
00017 #include <cstring>
00018 #include <cmath>  // for M_PI
00019 #include <ctime>
00020 #include "moab/Core.hpp"
00021 #include "moab/Interface.hpp"
00022 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
00023 #include "moab/ProgOptions.hpp"
00024 #include "MBTagConventions.hpp"
00025 #include "moab/ParallelComm.hpp"
00026 #include "moab/IntxMesh/IntxUtils.hpp"
00027 #include "IntxUtilsCSLAM.hpp"
00028 #include "TestUtil.hpp"
00029 
00030 using namespace moab;
00031 // some input data
00032 double gtol = 1.e-9;  // this is for geometry tolerance
00033 
00034 double radius = 1.;  // in m:  6371220.
00035 
00036 double t = 0.1, delta_t = 0.05;  // check the script
00037 bool Verbose = false;
00038 double rot   = M_PI / 10;
00039 
00040 ErrorCode manufacture_lagrange_mesh_on_sphere( Interface* mb, EntityHandle euler_set )
00041 {
00042     ErrorCode rval = MB_SUCCESS;
00043 
00044     /*
00045      * get all plys first, then vertices, then move them on the surface of the sphere
00046      *  radius is 1., most of the time
00047      *
00048      */
00049     Range polygons;
00050     rval = mb->get_entities_by_dimension( euler_set, 2, polygons );CHECK_ERR( rval );
00051 
00052     Range connecVerts;
00053     rval = mb->get_connectivity( polygons, connecVerts );CHECK_ERR( rval );
00054 
00055     Tag tagh = 0;
00056     std::string tag_name( "DP" );
00057     rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval );
00058     void* data;  // pointer to the LOC in memory, for each vertex
00059     int count;
00060 
00061     rval = mb->tag_iterate( tagh, connecVerts.begin(), connecVerts.end(), count, data );CHECK_ERR( rval );
00062     // here we are checking contiguity
00063     assert( count == (int)connecVerts.size() );
00064     double* ptr_DP = (double*)data;
00065     // get the coordinates of the old mesh, and move it around the sphere in the same way as in the
00066     // python script
00067 
00068     // now put the vertices in the right place....
00069     // int vix=0; // vertex index in new array
00070     double T = 5;  // check the script
00071 
00072     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00073     {
00074         EntityHandle oldV = *vit;
00075         CartVect posi;
00076         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
00077         // do some mumbo jumbo, as in python script
00078         IntxUtils::SphereCoords sphCoord = IntxUtils::cart_to_spherical( posi );
00079         double lat1                      = sphCoord.lat - 2 * M_PI * t / T;  // 0.1/5
00080         double uu = 3 * radius / T * pow( sin( lat1 ), 2 ) * sin( 2 * sphCoord.lon ) * cos( M_PI * t / T );
00081         uu += 2 * radius * M_PI * cos( sphCoord.lon ) / T;
00082         double vv = 3 * radius / T * ( sin( 2 * lat1 ) ) * cos( sphCoord.lon ) * cos( M_PI * t / T );
00083         double vx = -uu * sin( sphCoord.lon ) - vv * sin( sphCoord.lat ) * cos( sphCoord.lon );
00084         double vy = -uu * cos( sphCoord.lon ) - vv * sin( sphCoord.lat ) * sin( sphCoord.lon );
00085         double vz = vv * cos( sphCoord.lat );
00086         posi      = posi + delta_t * CartVect( vx, vy, vz );
00087         double x2 = posi[0] * cos( rot ) - posi[1] * sin( rot );
00088         double y2 = posi[0] * sin( rot ) + posi[1] * cos( rot );
00089         CartVect newPos( x2, y2, posi[2] );
00090         double len1 = newPos.length();
00091         newPos      = radius * newPos / len1;
00092 
00093         ptr_DP[0] = newPos[0];
00094         ptr_DP[1] = newPos[1];
00095         ptr_DP[2] = newPos[2];
00096         ptr_DP += 3;  // increment to the next node
00097     }
00098     return rval;
00099 }
00100 
00101 int main( int argc, char** argv )
00102 {
00103 
00104     MPI_Init( &argc, &argv );
00105 
00106     std::string extra_read_opts;
00107     std::string fileN          = TestDir + "unittest/io/mpasx1.642.t.2.nc";
00108     const char* filename_mesh1 = fileN.c_str();
00109     bool flux_form             = false;
00110     if( argc > 1 )
00111     {
00112         int index = 1;
00113         while( index < argc )
00114         {
00115             if( !strcmp( argv[index], "-gtol" ) )  // this is for geometry tolerance
00116             {
00117                 gtol = atof( argv[++index] );
00118             }
00119             if( !strcmp( argv[index], "-dt" ) )
00120             {
00121                 delta_t = atof( argv[++index] );
00122             }
00123             if( !strcmp( argv[index], "-input" ) )
00124             {
00125                 filename_mesh1 = argv[++index];
00126             }
00127             if( !strcmp( argv[index], "-R" ) )
00128             {
00129                 radius = atof( argv[++index] );
00130             }
00131             if( !strcmp( argv[index], "-O" ) )
00132             {
00133                 extra_read_opts = std::string( argv[++index] );
00134             }
00135             if( !strcmp( argv[index], "-FF" ) )
00136             {
00137                 flux_form = true;
00138             }
00139             if( !strcmp( argv[index], "-v" ) )
00140             {
00141                 Verbose = true;
00142             }
00143             if( !strcmp( argv[index], "-t" ) )
00144             {
00145                 t = atof( argv[++index] );
00146             }
00147             if( !strcmp( argv[index], "-t" ) )
00148             {
00149                 t = atof( argv[++index] );
00150             }
00151             if( !strcmp( argv[index], "-rot" ) )
00152             {
00153                 rot = atof( argv[++index] );
00154                 rot = M_PI / rot;  // so rot 50 means rotate with M_PI/50 radians
00155             }
00156 
00157             index++;
00158         }
00159     }
00160     // start copy
00161     std::string opts = std::string( "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN" ) +
00162                        std::string( ";PARALLEL_RESOLVE_SHARED_ENTS;VARIABLE=;NO_EDGES;" ) + extra_read_opts;
00163     Core moab;
00164     Interface& mb = moab;
00165     EntityHandle euler_set;
00166     ErrorCode rval;
00167     rval = mb.create_meshset( MESHSET_SET, euler_set );CHECK_ERR( rval );
00168 
00169     clock_t tt = clock();
00170 
00171     rval = mb.load_file( filename_mesh1, &euler_set, opts.c_str() );CHECK_ERR( rval );
00172 
00173     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );CHECK_ERR( rval );
00174 
00175     /*rval = pcomm->check_all_shared_handles();CHECK_ERR(rval);*/
00176     // end copy
00177     int rank  = pcomm->proc_config().proc_rank();
00178     int procs = pcomm->proc_config().proc_size();
00179 
00180     if( 0 == rank )
00181         std::cout << " case 1: use -gtol " << gtol << " -dt " << delta_t << " -R " << radius << " -input "
00182                   << filename_mesh1 << " -t " << t << " -rot " << rot << "\n";
00183 
00184     if( 0 == rank )
00185     {
00186         std::cout << "load mesh from " << filename_mesh1 << "\n  on " << procs << " processors in "
00187                   << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << std::endl;
00188         tt = clock();
00189     }
00190 
00191     rval = manufacture_lagrange_mesh_on_sphere( &mb, euler_set );CHECK_ERR( rval );
00192 
00193     // create a set with quads corresponding to each initial edge spanned with the displacement
00194     // field
00195     if( flux_form )
00196     {
00197         rval = IntxUtilsCSLAM::create_span_quads( &mb, euler_set, rank );CHECK_ERR( rval );
00198     }
00199 
00200     EntityHandle covering_lagr_set;
00201     rval = mb.create_meshset( MESHSET_SET, covering_lagr_set );CHECK_ERR( rval );
00202 
00203     Intx2MeshOnSphere worker( &mb );
00204     // double radius = 1.; // input
00205     worker.set_radius_source_mesh( radius );
00206     worker.set_radius_destination_mesh( radius );
00207     worker.set_parallel_comm( pcomm );
00208     if( 0 == rank )
00209     {
00210         std::cout << "manufacture departure mesh " << filename_mesh1 << "\n  on " << procs << " processors in "
00211                   << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << std::endl;
00212         tt = clock();
00213     }
00214     rval = worker.FindMaxEdges( euler_set, euler_set );CHECK_ERR( rval );
00215     worker.set_error_tolerance( gtol );
00216 
00217     rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );CHECK_ERR( rval );
00218 
00219     if( 0 == rank )
00220     {
00221         std::cout << "communicate covering mesh on " << procs << " processors in "
00222                   << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << std::endl;
00223         tt = clock();
00224     }
00225 
00226     if( Verbose )
00227     {
00228         std::stringstream lagrIni;
00229         lagrIni << "lagr0" << rank << ".h5m";
00230         rval = mb.write_file( lagrIni.str().c_str(), 0, 0, &covering_lagr_set, 1 );
00231     }
00232 
00233     rval = IntxUtils::enforce_convexity( &mb, covering_lagr_set, rank );CHECK_ERR( rval );
00234     if( Verbose )
00235     {
00236         std::stringstream ste;
00237         ste << "euler0" << rank << ".h5m";
00238         rval = mb.write_file( ste.str().c_str(), 0, 0, &euler_set, 1 );
00239     }
00240 
00241     if( MB_SUCCESS != rval ) std::cout << "can't write lagr set\n";
00242 
00243     EntityHandle outputSet;
00244     rval = mb.create_meshset( MESHSET_SET, outputSet );CHECK_ERR( rval );
00245     rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );CHECK_ERR( rval );
00246 
00247     if( 0 == rank )
00248     {
00249         std::cout << "intersect meshes in " << procs << " processors in " << ( clock() - tt ) / (double)CLOCKS_PER_SEC
00250                   << " seconds" << std::endl;
00251         tt = clock();
00252     }
00253     if( Verbose && rank <= 4 )
00254     {
00255         std::string opts_write( "" );
00256         std::stringstream outf;
00257         outf << "intersect0" << rank << ".h5m";
00258         rval = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );
00259         if( MB_SUCCESS != rval ) std::cout << "can't write output\n";
00260     }
00261 
00262     if( rank <= 4 )
00263     {
00264         moab::IntxAreaUtils areaAdaptor;
00265         double intx_area    = areaAdaptor.area_on_sphere( &mb, outputSet, radius, rank );
00266         double arrival_area = areaAdaptor.area_on_sphere( &mb, euler_set, radius, rank );
00267         std::cout << "On proc " << rank << "  arrival area: " << arrival_area << "  intersection area:" << intx_area
00268                   << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";
00269     }
00270 
00271     MPI_Finalize();
00272     if( MB_SUCCESS != rval ) return 1;
00273 
00274     return 0;
00275 }
00276 
00277 #else
00278 
00279 int main( int /*argc*/, char** /*argv*/ )
00280 {
00281     return 0;
00282 }
00283 
00284 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines