MOAB: Mesh Oriented datABase  (version 5.3.0)
mbtempest.cpp
Go to the documentation of this file.
00001 /*
00002  * Usage: MOAB-Tempest tool
00003  *
00004  * Generate a Cubed-Sphere mesh: ./mbtempest -t 0 -res 25 -f cubed_sphere_mesh.exo
00005  * Generate a RLL mesh: ./mbtempest -t 1 -res 25 -f rll_mesh.exo
00006  * Generate a Icosahedral-Sphere mesh: ./mbtempest -t 2 -res 25 <-dual> -f icosa_mesh.exo
00007  *
00008  * Now you can compute the intersections between the meshes too!
00009  *
00010  * Generate the overlap mesh: ./mbtempest -t 3 -l cubed_sphere_mesh.exo -l rll_mesh.exo -f
00011  * overlap_mesh.exo
00012  *
00013  */
00014 
00015 #include <iostream>
00016 #include <cstdlib>
00017 #include <vector>
00018 #include <string>
00019 #include <sstream>
00020 #include <cassert>
00021 
00022 #include "moab/Core.hpp"
00023 #include "moab/IntxMesh/IntxUtils.hpp"
00024 #include "moab/Remapping/TempestRemapper.hpp"
00025 #include "moab/Remapping/TempestOnlineMap.hpp"
00026 #include "moab/ProgOptions.hpp"
00027 #include "moab/CpuTimer.hpp"
00028 #include "DebugOutput.hpp"
00029 
00030 //#ifndef MOAB_HAVE_MPI
00031 //    #error mbtempest tool requires MPI configuration
00032 //#endif
00033 
00034 #ifdef MOAB_HAVE_MPI
00035 // MPI includes
00036 #include "moab_mpi.h"
00037 #include "moab/ParallelComm.hpp"
00038 #include "MBParallelConventions.h"
00039 #endif
00040 
00041 struct ToolContext
00042 {
00043     moab::Interface* mbcore;
00044 #ifdef MOAB_HAVE_MPI
00045     moab::ParallelComm* pcomm;
00046 #endif
00047     const int proc_id, n_procs;
00048     moab::DebugOutput outputFormatter;
00049     int blockSize;
00050     std::vector< std::string > inFilenames;
00051     std::vector< Mesh* > meshes;
00052     std::vector< moab::EntityHandle > meshsets;
00053     std::vector< int > disc_orders;
00054     std::vector< std::string > disc_methods;
00055     std::vector< std::string > doftag_names;
00056     std::string outFilename;
00057     std::string intxFilename;
00058     moab::TempestRemapper::TempestMeshType meshType;
00059     bool computeDual;
00060     bool computeWeights;
00061     bool verifyWeights;
00062     bool enforceConvexity;
00063     int ensureMonotonicity;
00064     bool fNoConservation;
00065     bool fVolumetric;
00066     bool rrmGrids;
00067     bool kdtreeSearch;
00068     bool fNoBubble, fInputConcave, fOutputConcave, fCheck;
00069 
00070 #ifdef MOAB_HAVE_MPI
00071     ToolContext( moab::Interface* icore, moab::ParallelComm* p_pcomm )
00072         : mbcore( icore ), pcomm( p_pcomm ), proc_id( pcomm->rank() ), n_procs( pcomm->size() ),
00073           outputFormatter( std::cout, pcomm->rank(), 0 ),
00074 #else
00075     ToolContext( moab::Interface* icore )
00076         : mbcore( icore ), proc_id( 0 ), n_procs( 1 ), outputFormatter( std::cout, 0, 0 ),
00077 #endif
00078           blockSize( 5 ), outFilename( "output.exo" ), intxFilename( "" ), meshType( moab::TempestRemapper::DEFAULT ),
00079           computeDual( false ), computeWeights( false ), verifyWeights( false ), enforceConvexity( false ),
00080           ensureMonotonicity( 0 ), fNoConservation( false ), fVolumetric( false ), rrmGrids( false ),
00081           kdtreeSearch( true ), fNoBubble( true ), fInputConcave( false ), fOutputConcave( false ),
00082           fCheck( n_procs > 1 ? false : true )
00083     {
00084         inFilenames.resize( 2 );
00085         doftag_names.resize( 2 );
00086         timer = new moab::CpuTimer();
00087 
00088         outputFormatter.set_prefix( "[MBTempest]: " );
00089     }
00090 
00091     ~ToolContext()
00092     {
00093         // for (unsigned i=0; i < meshes.size(); ++i) delete meshes[i];
00094         meshes.clear();
00095         inFilenames.clear();
00096         disc_orders.clear();
00097         disc_methods.clear();
00098         doftag_names.clear();
00099         outFilename.clear();
00100         intxFilename.clear();
00101         meshsets.clear();
00102         delete timer;
00103     }
00104 
00105     void timer_push( std::string operation )
00106     {
00107         timer_ops = timer->time_since_birth();
00108         opName    = operation;
00109     }
00110 
00111     void timer_pop()
00112     {
00113         double locElapsed = timer->time_since_birth() - timer_ops, avgElapsed = 0, maxElapsed = 0;
00114 #ifdef MOAB_HAVE_MPI
00115         MPI_Reduce( &locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm() );
00116         MPI_Reduce( &locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, pcomm->comm() );
00117 #else
00118         maxElapsed = locElapsed;
00119         avgElapsed = locElapsed;
00120 #endif
00121         if( !proc_id )
00122         {
00123             avgElapsed /= n_procs;
00124             std::cout << "[LOG] Time taken to " << opName.c_str() << ": max = " << maxElapsed
00125                       << ", avg = " << avgElapsed << "\n";
00126         }
00127         // std::cout << "\n[LOG" << proc_id << "] Time taken to " << opName << " = " <<
00128         // timer->time_since_birth() - timer_ops << std::endl;
00129         opName.clear();
00130     }
00131 
00132     void ParseCLOptions( int argc, char* argv[] )
00133     {
00134         ProgOptions opts;
00135         int imeshType                  = 0;
00136         std::string expectedFName      = "output.exo";
00137         std::string expectedMethod     = "fv";
00138         std::string expectedDofTagName = "GLOBAL_ID";
00139         int expectedOrder              = 1;
00140         bool fBubble                   = false;
00141 
00142         if( !proc_id )
00143         {
00144             std::cout << "Command line options provided to mbtempest:\n  ";
00145             for( int iarg = 0; iarg < argc; ++iarg )
00146                 std::cout << argv[iarg] << " ";
00147             std::cout << std::endl << std::endl;
00148         }
00149 
00150         opts.addOpt< int >( "type,t",
00151                             "Type of mesh (default=CS; Choose from [CS=0, RLL=1, ICO=2, OVERLAP_FILES=3, "
00152                             "OVERLAP_MEMORY=4, OVERLAP_MOAB=5])",
00153                             &imeshType );
00154         opts.addOpt< int >( "res,r", "Resolution of the mesh (default=5)", &blockSize );
00155         opts.addOpt< void >( "dual,d", "Output the dual of the mesh (relevant only for ICO mesh type)", &computeDual );
00156         opts.addOpt< std::string >( "file,f", "Output computed mesh or remapping weights to specified filename",
00157                                     &outFilename );
00158         opts.addOpt< std::string >(
00159             "load,l", "Input mesh filenames for source and target meshes. (relevant only when computing weights)",
00160             &expectedFName );
00161         opts.addOpt< void >( "noconserve,c",
00162                              "Do not apply conservation to the resultant weights (relevant only "
00163                              "when computing weights)",
00164                              &fNoConservation );
00165         opts.addOpt< void >( "volumetric,v",
00166                              "Apply a volumetric projection to compute the weights (relevant only "
00167                              "when computing weights)",
00168                              &fVolumetric );
00169         opts.addOpt< void >( "advfront,a",
00170                              "Use the advancing front intersection instead of the Kd-tree based algorithm "
00171                              "to compute mesh intersections. (relevant only when computing weights)" );
00172         opts.addOpt< int >( "monotonic,n", "Ensure monotonicity in the weight generation", &ensureMonotonicity );
00173         opts.addOpt< int >( "order,o", "Discretization orders for the source and target solution fields",
00174                             &expectedOrder );
00175         opts.addOpt< std::string >( "method,m", "Discretization method for the source and target solution fields",
00176                                     &expectedMethod );
00177         opts.addOpt< std::string >( "global_id,g",
00178                                     "Tag name that contains the global DoF IDs for source and target solution fields",
00179                                     &expectedDofTagName );
00180 
00181         opts.addOpt< void >( "weights,w",
00182                              "Compute and output the weights using the overlap mesh (generally "
00183                              "relevant only for OVERLAP mesh)",
00184                              &computeWeights );
00185         opts.addOpt< std::string >( "intx,i", "Output TempestRemap intersection mesh filename", &intxFilename );
00186 
00187         opts.addOpt< void >( "rrmgrids",
00188                              "At least one of the meshes is a regionally refined grid (relevant to "
00189                              "accelerate intersection computation)",
00190                              &rrmGrids );
00191         opts.addOpt< void >( "checkmap", "Check the generated map for conservation and consistency", &fCheck );
00192         opts.addOpt< void >( "verify",
00193                              "Verify the accuracy of the maps by projecting analytical functions "
00194                              "from source to target "
00195                              "grid by applying the maps",
00196                              &verifyWeights );
00197         opts.addOpt< void >( "enforce_convexity", "check convexity of input meshes to compute mesh intersections",
00198                              &enforceConvexity );
00199         opts.addOpt< void >( "bubble", "use bubble on interior of spectral element nodes", &fBubble );
00200 
00201         opts.parseCommandLine( argc, argv );
00202 
00203         // By default - use Kd-tree based search; if user asks for advancing front, disable Kd-tree
00204         // algorithm
00205         kdtreeSearch = opts.numOptSet( "advfront,a" ) == 0;
00206 
00207         // negate the option
00208         fNoBubble = !fBubble;
00209 
00210         switch( imeshType )
00211         {
00212             case 0:
00213                 meshType = moab::TempestRemapper::CS;
00214                 break;
00215 
00216             case 1:
00217                 meshType = moab::TempestRemapper::RLL;
00218                 break;
00219 
00220             case 2:
00221                 meshType = moab::TempestRemapper::ICO;
00222                 break;
00223 
00224             case 3:
00225                 meshType = moab::TempestRemapper::OVERLAP_FILES;
00226                 break;
00227 
00228             case 4:
00229                 meshType = moab::TempestRemapper::OVERLAP_MEMORY;
00230                 break;
00231 
00232             case 5:
00233                 meshType = moab::TempestRemapper::OVERLAP_MOAB;
00234                 break;
00235 
00236             default:
00237                 meshType = moab::TempestRemapper::DEFAULT;
00238                 break;
00239         }
00240 
00241         if( meshType > moab::TempestRemapper::ICO )
00242         {
00243             opts.getOptAllArgs( "load,l", inFilenames );
00244             opts.getOptAllArgs( "order,o", disc_orders );
00245             opts.getOptAllArgs( "method,m", disc_methods );
00246             opts.getOptAllArgs( "global_id,i", doftag_names );
00247 
00248             if( disc_orders.size() == 0 ) { disc_orders.resize( 2, 1 ); }
00249 
00250             if( disc_orders.size() == 1 ) { disc_orders.push_back( 1 ); }
00251 
00252             if( disc_methods.size() == 0 ) { disc_methods.resize( 2, "fv" ); }
00253 
00254             if( disc_methods.size() == 1 ) { disc_methods.push_back( "fv" ); }
00255 
00256             if( doftag_names.size() == 0 ) { doftag_names.resize( 2, "GLOBAL_ID" ); }
00257 
00258             if( doftag_names.size() == 1 ) { doftag_names.push_back( "GLOBAL_ID" ); }
00259 
00260             assert( inFilenames.size() == 2 );
00261             assert( disc_orders.size() == 2 );
00262             assert( disc_methods.size() == 2 );
00263         }
00264 
00265         expectedFName.clear();
00266     }
00267 
00268   private:
00269     moab::CpuTimer* timer;
00270     double timer_ops;
00271     std::string opName;
00272 };
00273 
00274 // Forward declare some methods
00275 static moab::ErrorCode CreateTempestMesh( ToolContext&, moab::TempestRemapper& remapper, Mesh* );
00276 inline double sample_slow_harmonic( double dLon, double dLat );
00277 inline double sample_fast_harmonic( double dLon, double dLat );
00278 inline double sample_constant( double dLon, double dLat );
00279 inline double sample_stationary_vortex( double dLon, double dLat );
00280 
00281 int main( int argc, char* argv[] )
00282 {
00283     moab::ErrorCode rval;
00284     NcError error( NcError::verbose_nonfatal );
00285     std::stringstream sstr;
00286 
00287     int proc_id = 0, nprocs = 1;
00288 #ifdef MOAB_HAVE_MPI
00289     MPI_Init( &argc, &argv );
00290     MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
00291     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00292 #endif
00293 
00294     moab::Interface* mbCore = new( std::nothrow ) moab::Core;
00295 
00296     if( NULL == mbCore ) { return 1; }
00297 
00298     ToolContext* runCtx;
00299 #ifdef MOAB_HAVE_MPI
00300     moab::ParallelComm* pcomm = new moab::ParallelComm( mbCore, MPI_COMM_WORLD, 0 );
00301 
00302     runCtx                   = new ToolContext( mbCore, pcomm );
00303     const char* writeOptions = ( nprocs > 1 ? "PARALLEL=WRITE_PART" : "" );
00304 #else
00305     runCtx                   = new ToolContext( mbCore );
00306     const char* writeOptions = "";
00307 #endif
00308     runCtx->ParseCLOptions( argc, argv );
00309 
00310     const double radius_src  = 1.0 /*2.0*acos(-1.0)*/;
00311     const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
00312 
00313     moab::DebugOutput& outputFormatter = runCtx->outputFormatter;
00314 
00315 #ifdef MOAB_HAVE_MPI
00316     moab::TempestRemapper remapper( mbCore, pcomm );
00317 #else
00318     moab::TempestRemapper remapper( mbCore );
00319 #endif
00320     remapper.meshValidate     = true;
00321     remapper.constructEdgeMap = false;
00322     remapper.initialize();
00323 
00324     // Default area_method = lHuiller; Options: Girard, GaussQuadrature (if TR is available)
00325 #ifdef MOAB_HAVE_TEMPESTREMAP
00326     moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::GaussQuadrature );
00327 #else
00328     moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::lHuiller );
00329 #endif
00330 
00331     Mesh* tempest_mesh = new Mesh();
00332     runCtx->timer_push( "create Tempest mesh" );
00333     rval = CreateTempestMesh( *runCtx, remapper, tempest_mesh );MB_CHK_ERR( rval );
00334     runCtx->timer_pop();
00335 
00336     const double epsrel = ReferenceTolerance;  // ReferenceTolerance is defined in Defines.h in tempestremap;
00337                                                // Defines.h is included in SparseMatrix.h
00338                                                // SparseMatrix.h is included in OfflineMap.h
00339                                                // OfflineMap.h is included in TempestOnlineMap.hpp
00340                                                // TempestOnlineMap.hpp is included in this file, and is part of MOAB
00341     // Some constant parameters
00342 
00343     const double boxeps = 1e-6;
00344 
00345     if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MEMORY )
00346     {
00347         // Compute intersections with MOAB
00348         // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
00349         assert( runCtx->meshes.size() == 3 );
00350 
00351 #ifdef MOAB_HAVE_MPI
00352         rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
00353 #endif
00354 
00355         // Load the meshes and validate
00356         rval = remapper.ConvertTempestMesh( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
00357         rval = remapper.ConvertTempestMesh( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
00358         remapper.SetMeshType( moab::Remapper::OverlapMesh, moab::TempestRemapper::OVERLAP_FILES );
00359         rval = remapper.ConvertTempestMesh( moab::Remapper::OverlapMesh );MB_CHK_ERR( rval );
00360         rval = mbCore->write_mesh( "tempest_intersection.h5m", &runCtx->meshsets[2], 1 );MB_CHK_ERR( rval );
00361 
00362         // print verbosely about the problem setting
00363         {
00364             moab::Range rintxverts, rintxelems;
00365             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval );
00366             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval );
00367             outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", rintxverts.size(),
00368                                     rintxelems.size() );
00369 
00370             moab::Range bintxverts, bintxelems;
00371             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval );
00372             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval );
00373             outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", bintxverts.size(),
00374                                     bintxelems.size() );
00375         }
00376 
00377         moab::EntityHandle intxset;  // == remapper.GetMeshSet(moab::Remapper::OverlapMesh);
00378 
00379         // Compute intersections with MOAB
00380         {
00381             // Create the intersection on the sphere object
00382             runCtx->timer_push( "setup the intersector" );
00383 
00384             moab::Intx2MeshOnSphere* mbintx = new moab::Intx2MeshOnSphere( mbCore );
00385             mbintx->set_error_tolerance( epsrel );
00386             mbintx->set_box_error( boxeps );
00387             mbintx->set_radius_source_mesh( radius_src );
00388             mbintx->set_radius_destination_mesh( radius_dest );
00389 #ifdef MOAB_HAVE_MPI
00390             mbintx->set_parallel_comm( pcomm );
00391 #endif
00392             rval = mbintx->FindMaxEdges( runCtx->meshsets[0], runCtx->meshsets[1] );MB_CHK_ERR( rval );
00393 
00394 #ifdef MOAB_HAVE_MPI
00395             moab::Range local_verts;
00396             rval = mbintx->build_processor_euler_boxes( runCtx->meshsets[1], local_verts );MB_CHK_ERR( rval );
00397 
00398             runCtx->timer_pop();
00399 
00400             moab::EntityHandle covering_set;
00401             runCtx->timer_push( "communicate the mesh" );
00402             rval = mbintx->construct_covering_set( runCtx->meshsets[0], covering_set );MB_CHK_ERR( rval );  // lots of communication if mesh is distributed very differently
00403             runCtx->timer_pop();
00404 #else
00405             moab::EntityHandle covering_set = runCtx->meshsets[0];
00406 #endif
00407             // Now let's invoke the MOAB intersection algorithm in parallel with a
00408             // source and target mesh set representing two different decompositions
00409             runCtx->timer_push( "compute intersections with MOAB" );
00410             rval = mbCore->create_meshset( moab::MESHSET_SET, intxset );MB_CHK_SET_ERR( rval, "Can't create new set" );
00411             rval = mbintx->intersect_meshes( covering_set, runCtx->meshsets[1], intxset );MB_CHK_SET_ERR( rval, "Can't compute the intersection of meshes on the sphere" );
00412             runCtx->timer_pop();
00413 
00414             // free the memory
00415             delete mbintx;
00416         }
00417 
00418         {
00419             moab::Range intxelems, intxverts;
00420             rval = mbCore->get_entities_by_dimension( intxset, 2, intxelems );MB_CHK_ERR( rval );
00421             rval = mbCore->get_entities_by_dimension( intxset, 0, intxverts, true );MB_CHK_ERR( rval );
00422             outputFormatter.printf( 0, "The intersection set contains %lu elements and %lu vertices \n",
00423                                     intxelems.size(), intxverts.size() );
00424 
00425             double initial_sarea =
00426                 areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0],
00427                                             radius_src );  // use the target to compute the initial area
00428             double initial_tarea =
00429                 areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1],
00430                                             radius_dest );  // use the target to compute the initial area
00431             double intx_area = areaAdaptor.area_on_sphere( mbCore, intxset, radius_src );
00432 
00433             outputFormatter.printf( 0, "mesh areas: source = %12.10f, target = %12.10f, intersection = %12.10f \n",
00434                                     initial_sarea, initial_tarea, intx_area );
00435             outputFormatter.printf( 0, "relative error w.r.t source = %12.10e, target = %12.10e \n",
00436                                     fabs( intx_area - initial_sarea ) / initial_sarea,
00437                                     fabs( intx_area - initial_tarea ) / initial_tarea );
00438         }
00439 
00440         // Write out our computed intersection file
00441         rval = mbCore->write_mesh( "moab_intersection.h5m", &intxset, 1 );MB_CHK_ERR( rval );
00442 
00443         if( runCtx->computeWeights )
00444         {
00445             runCtx->timer_push( "compute weights with the Tempest meshes" );
00446             // Call to generate an offline map with the tempest meshes
00447             OfflineMap weightMap;
00448             int err = GenerateOfflineMapWithMeshes(
00449                 weightMap, *runCtx->meshes[0], *runCtx->meshes[1], *runCtx->meshes[2], "",
00450                 "",  // std::string strInputMeta, std::string strOutputMeta,
00451                 runCtx->disc_methods[0],
00452                 runCtx->disc_methods[1],                         // std::string strInputType, std::string strOutputType,
00453                 runCtx->disc_orders[0], runCtx->disc_orders[1],  // int nPin=4, int nPout=4,
00454                 runCtx->fNoBubble, true,                         // bool fNoBubble = false, bool fCorrectAreas = false,
00455                 runCtx->ensureMonotonicity                       // int fMonotoneTypeID = 0
00456             );
00457             runCtx->timer_pop();
00458 
00459             std::map< std::string, std::string > mapAttributes;
00460             if( err ) { rval = moab::MB_FAILURE; }
00461             else
00462             {
00463                 weightMap.Write( "outWeights.nc", mapAttributes );
00464             }
00465         }
00466     }
00467     else if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MOAB )
00468     {
00469         // Usage: mpiexec -n 2 tools/mbtempest -t 5 -l mycs_2.h5m -l myico_2.h5m -f myoverlap_2.h5m
00470 #ifdef MOAB_HAVE_MPI
00471         rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
00472 #endif
00473         // print verbosely about the problem setting
00474         {
00475             moab::Range rintxverts, rintxelems;
00476             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval );
00477             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval );
00478             rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[0] );MB_CHK_ERR( rval );
00479             if( runCtx->enforceConvexity )
00480             {
00481                 rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[0], proc_id );MB_CHK_ERR( rval );
00482             }
00483             rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[0], radius_src );MB_CHK_ERR( rval );
00484             if( !proc_id )
00485                 outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n",
00486                                         rintxverts.size(), rintxelems.size() );
00487 
00488             moab::Range bintxverts, bintxelems;
00489             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval );
00490             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval );
00491             rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[1] );MB_CHK_ERR( rval );
00492             if( runCtx->enforceConvexity )
00493             {
00494                 rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[1], proc_id );MB_CHK_ERR( rval );
00495             }
00496             rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[1], radius_dest );MB_CHK_ERR( rval );
00497             if( !proc_id )
00498                 outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n",
00499                                         bintxverts.size(), bintxelems.size() );
00500         }
00501 
00502         // First compute the covering set such that the target elements are fully covered by the
00503         // lcoal source grid
00504         runCtx->timer_push( "construct covering set for intersection" );
00505         rval = remapper.ConstructCoveringSet( epsrel, 1.0, 1.0, boxeps, runCtx->rrmGrids );MB_CHK_ERR( rval );
00506         runCtx->timer_pop();
00507 
00508         // Compute intersections with MOAB with either the Kd-tree or the advancing front algorithm
00509         runCtx->timer_push( "setup and compute mesh intersections" );
00510         rval = remapper.ComputeOverlapMesh( runCtx->kdtreeSearch, false );MB_CHK_ERR( rval );
00511         runCtx->timer_pop();
00512 
00513         // print some diagnostic checks to see if the overlap grid resolved the input meshes
00514         // correctly
00515         {
00516             double local_areas[3],
00517                 global_areas[3];  // Array for Initial area, and through Method 1 and Method 2
00518             // local_areas[0] = area_on_sphere_lHuiller ( mbCore, runCtx->meshsets[1], radius_src );
00519             local_areas[0] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0], radius_src );
00520             local_areas[1] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1], radius_dest );
00521             local_areas[2] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[2], radius_src );
00522 
00523 #ifdef MOAB_HAVE_MPI
00524             MPI_Allreduce( &local_areas[0], &global_areas[0], 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
00525 #else
00526             global_areas[0]                 = local_areas[0];
00527             global_areas[1]                 = local_areas[1];
00528             global_areas[2]                 = local_areas[2];
00529 #endif
00530             if( !proc_id )
00531             {
00532                 outputFormatter.printf( 0,
00533                                         "initial area: source mesh = %12.14f, target mesh = "
00534                                         "%12.14f, overlap mesh = %12.14f\n",
00535                                         global_areas[0], global_areas[1], global_areas[2] );
00536                 // outputFormatter.printf ( 0, " area with l'Huiller: %12.14f with Girard:
00537                 // %12.14f\n", global_areas[2], global_areas[3] ); outputFormatter.printf ( 0, "
00538                 // relative difference areas = %12.10e\n", fabs ( global_areas[2] - global_areas[3]
00539                 // ) / global_areas[2] );
00540                 outputFormatter.printf( 0, "relative error w.r.t source = %12.14e, and target = %12.14e\n",
00541                                         fabs( global_areas[0] - global_areas[2] ) / global_areas[0],
00542                                         fabs( global_areas[1] - global_areas[2] ) / global_areas[1] );
00543             }
00544         }
00545 
00546         if( runCtx->intxFilename.size() )
00547         {
00548             moab::EntityHandle writableOverlapSet;
00549             rval = mbCore->create_meshset( moab::MESHSET_SET, writableOverlapSet );MB_CHK_SET_ERR( rval, "Can't create new set" );
00550             moab::EntityHandle meshOverlapSet = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
00551             moab::Range ovEnts;
00552             rval = mbCore->get_entities_by_dimension( meshOverlapSet, 2, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
00553             rval = mbCore->get_entities_by_dimension( meshOverlapSet, 0, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
00554 
00555 #ifdef MOAB_HAVE_MPI
00556             // Do not remove ghosted entities if we still haven't computed weights
00557             // Remove ghosted entities from overlap set before writing the new mesh set to file
00558             if( nprocs > 1 )
00559             {
00560                 moab::Range ghostedEnts;
00561                 rval = remapper.GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
00562                 ovEnts = moab::subtract( ovEnts, ghostedEnts );
00563             }
00564 #endif
00565             rval = mbCore->add_entities( writableOverlapSet, ovEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" );
00566 
00567             size_t lastindex = runCtx->intxFilename.find_last_of( "." );
00568             sstr.str( "" );
00569             sstr << runCtx->intxFilename.substr( 0, lastindex ) << ".h5m";
00570             if( !runCtx->proc_id )
00571                 std::cout << "Writing out the MOAB intersection mesh file to " << sstr.str() << std::endl;
00572 
00573             // Write out our computed intersection file
00574             rval = mbCore->write_file( sstr.str().c_str(), NULL, writeOptions, &writableOverlapSet, 1 );MB_CHK_ERR( rval );
00575         }
00576 
00577         if( runCtx->computeWeights )
00578         {
00579             runCtx->meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
00580             if( !runCtx->proc_id ) std::cout << std::endl;
00581 
00582             runCtx->timer_push( "setup computation of weights" );
00583             // Call to generate the remapping weights with the tempest meshes
00584             moab::TempestOnlineMap* weightMap = new moab::TempestOnlineMap( &remapper );
00585             runCtx->timer_pop();
00586 
00587             runCtx->timer_push( "compute weights with TempestRemap" );
00588 
00589             rval = weightMap->GenerateRemappingWeights(
00590                 runCtx->disc_methods[0],                         // std::string strInputType
00591                 runCtx->disc_methods[1],                         // std::string strOutputType,
00592                 runCtx->disc_orders[0], runCtx->disc_orders[1],  // int nPin, int nPout,
00593                 runCtx->fNoBubble,                               // bool fNoBubble
00594                 runCtx->ensureMonotonicity,                      // int fMonotoneTypeID,
00595                 runCtx->fVolumetric, runCtx->fNoConservation,    // bool fVolumetric, bool fNoConservation,
00596                 !runCtx->fCheck,                                 //  bool fNoCheck,
00597                 runCtx->doftag_names[0],                         // std::string source_tag_name
00598                 runCtx->doftag_names[1],                         // std::string target_tag_name,
00599                 runCtx->fInputConcave,                           // bool fInputConcave
00600                 runCtx->fOutputConcave                           // bool fOutputConcave
00601             );MB_CHK_ERR( rval );
00602             runCtx->timer_pop();
00603 
00604             // Invoke the CheckMap routine on the TempestRemap serial interface directly, if running
00605             // on a single process
00606             if( nprocs == 1 )
00607             {
00608                 const double dNormalTolerance = 1.0E-8;
00609                 const double dStrictTolerance = 1.0E-12;
00610                 weightMap->CheckMap( runCtx->fCheck, runCtx->fCheck, runCtx->fCheck && ( runCtx->ensureMonotonicity ),
00611                                      dNormalTolerance, dStrictTolerance );
00612             }
00613 
00614             if( runCtx->outFilename.size() )
00615             {
00616                 // Write the map file to disk in parallel using either HDF5 or SCRIP interface
00617                 rval = weightMap->WriteParallelMap( runCtx->outFilename.c_str() );MB_CHK_ERR( rval );
00618 
00619                 // Write out the metadata information for the map file
00620                 if( proc_id == 0 )
00621                 {
00622                     size_t lastindex = runCtx->outFilename.find_last_of( "." );
00623                     sstr.str( "" );
00624                     sstr << runCtx->outFilename.substr( 0, lastindex ) << ".meta";
00625 
00626                     std::ofstream metafile( sstr.str() );
00627                     metafile << "Generator = MOAB-TempestRemap (mbtempest) Offline Regridding "
00628                                 "Weight Generator"
00629                              << std::endl;
00630                     metafile << "domain_a = " << runCtx->inFilenames[0] << std::endl;
00631                     metafile << "domain_b = " << runCtx->inFilenames[1] << std::endl;
00632                     metafile << "domain_aNb = "
00633                              << ( runCtx->intxFilename.size() ? runCtx->intxFilename : "outOverlap.h5m" ) << std::endl;
00634                     metafile << "map_aPb = " << runCtx->outFilename << std::endl;
00635                     metafile << "type_src = " << runCtx->disc_methods[0] << std::endl;
00636                     metafile << "np_src = " << runCtx->disc_orders[0] << std::endl;
00637                     metafile << "concave_src = " << ( runCtx->fInputConcave ? "true" : "false" ) << std::endl;
00638                     metafile << "type_dst = " << runCtx->disc_methods[1] << std::endl;
00639                     metafile << "np_dst = " << runCtx->disc_orders[1] << std::endl;
00640                     metafile << "concave_dst = " << ( runCtx->fOutputConcave ? "true" : "false" ) << std::endl;
00641                     metafile << "mono_type = " << runCtx->ensureMonotonicity << std::endl;
00642                     metafile << "bubble = " << ( runCtx->fNoBubble ? "false" : "true" ) << std::endl;
00643                     metafile << "version = "
00644                              << "MOAB v5.1.0+" << std::endl;
00645                     metafile.close();
00646                 }
00647             }
00648 
00649             if( runCtx->verifyWeights )
00650             {
00651                 // Let us pick a sampling test function for solution evaluation
00652                 moab::TempestOnlineMap::sample_function testFunction =
00653                     &sample_stationary_vortex;  // &sample_slow_harmonic;
00654 
00655                 runCtx->timer_push( "describe a solution on source grid" );
00656                 moab::Tag srcAnalyticalFunction;
00657                 rval = weightMap->DefineAnalyticalSolution( srcAnalyticalFunction, "AnalyticalSolnSrcExact",
00658                                                             moab::Remapper::SourceMesh, testFunction );MB_CHK_ERR( rval );
00659                 runCtx->timer_pop();
00660                 // rval = mbCore->write_file ( "srcWithSolnTag.h5m", NULL, writeOptions,
00661                 // &runCtx->meshsets[0], 1 ); MB_CHK_ERR ( rval );
00662 
00663                 runCtx->timer_push( "describe a solution on target grid" );
00664                 moab::Tag tgtAnalyticalFunction;
00665                 moab::Tag tgtProjectedFunction;
00666                 rval = weightMap->DefineAnalyticalSolution( tgtAnalyticalFunction, "AnalyticalSolnTgtExact",
00667                                                             moab::Remapper::TargetMesh, testFunction,
00668                                                             &tgtProjectedFunction, "ProjectedSolnTgt" );MB_CHK_ERR( rval );
00669                 // rval = mbCore->write_file ( "tgtWithSolnTag.h5m", NULL, writeOptions,
00670                 // &runCtx->meshsets[1], 1 ); MB_CHK_ERR ( rval );
00671                 runCtx->timer_pop();
00672 
00673                 runCtx->timer_push( "compute solution projection on target grid" );
00674                 rval = weightMap->ApplyWeights( srcAnalyticalFunction, tgtProjectedFunction );MB_CHK_ERR( rval );
00675                 runCtx->timer_pop();
00676                 rval = mbCore->write_file( "tgtWithSolnTag2.h5m", NULL, writeOptions, &runCtx->meshsets[1], 1 );MB_CHK_ERR( rval );
00677 
00678                 runCtx->timer_push( "compute error metrics against analytical solution on target grid" );
00679                 std::map< std::string, double > errMetrics;
00680                 rval = weightMap->ComputeMetrics( moab::Remapper::TargetMesh, tgtAnalyticalFunction,
00681                                                   tgtProjectedFunction, errMetrics, true );MB_CHK_ERR( rval );
00682                 runCtx->timer_pop();
00683             }
00684 
00685             delete weightMap;
00686         }
00687     }
00688 
00689     // Clean up
00690     remapper.clear();
00691     delete runCtx;
00692     delete mbCore;
00693 
00694 #ifdef MOAB_HAVE_MPI
00695     MPI_Finalize();
00696 #endif
00697     exit( 0 );
00698 }
00699 
00700 static moab::ErrorCode CreateTempestMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
00701 {
00702     moab::ErrorCode rval = moab::MB_SUCCESS;
00703     int err;
00704     moab::DebugOutput& outputFormatter = ctx.outputFormatter;
00705 
00706     if( !ctx.proc_id ) { outputFormatter.printf( 0, "Creating TempestRemap Mesh object ...\n" ); }
00707 
00708     if( ctx.meshType == moab::TempestRemapper::OVERLAP_FILES )
00709     {
00710         // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
00711         err = GenerateOverlapMesh( ctx.inFilenames[0], ctx.inFilenames[1], *tempest_mesh, ctx.outFilename, "NetCDF4",
00712                                    "exact", true );
00713 
00714         if( err ) { rval = moab::MB_FAILURE; }
00715         else
00716         {
00717             ctx.meshes.push_back( tempest_mesh );
00718         }
00719     }
00720     else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MEMORY )
00721     {
00722         // Load the meshes and validate
00723         ctx.meshsets.resize( 3 );
00724         ctx.meshes.resize( 3 );
00725         ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
00726         ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
00727         ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
00728 
00729         // First the source
00730         rval = remapper.LoadMesh( moab::Remapper::SourceMesh, ctx.inFilenames[0], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
00731         ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
00732 
00733         // Next the target
00734         rval = remapper.LoadMesh( moab::Remapper::TargetMesh, ctx.inFilenames[1], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
00735         ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
00736 
00737         // Now let us construct the overlap mesh, by calling TempestRemap interface directly
00738         // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
00739         err = GenerateOverlapWithMeshes( *ctx.meshes[0], *ctx.meshes[1], *tempest_mesh, "" /*ctx.outFilename*/,
00740                                          "NetCDF4", "exact", false );
00741 
00742         if( err ) { rval = moab::MB_FAILURE; }
00743         else
00744         {
00745             remapper.SetMesh( moab::Remapper::OverlapMesh, tempest_mesh );
00746             ctx.meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
00747             // ctx.meshes.push_back(*tempest_mesh);
00748         }
00749     }
00750     else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MOAB )
00751     {
00752         ctx.meshsets.resize( 3 );
00753         ctx.meshes.resize( 3 );
00754         ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
00755         ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
00756         ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
00757 
00758         const double radius_src  = 1.0 /*2.0*acos(-1.0)*/;
00759         const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
00760 
00761         const char* additional_read_opts = ( ctx.n_procs > 1 ? "NO_SET_CONTAINING_PARENTS;" : "" );
00762         std::vector< int > smetadata, tmetadata;
00763         // Load the source mesh and validate
00764         rval = remapper.LoadNativeMesh( ctx.inFilenames[0], ctx.meshsets[0], smetadata, additional_read_opts );MB_CHK_ERR( rval );
00765         if( smetadata.size() )
00766         {
00767             // remapper.SetMeshType( moab::Remapper::SourceMesh,
00768             //                       static_cast< moab::TempestRemapper::TempestMeshType >( smetadata[0] ), &smetadata
00769             //                       );
00770             remapper.SetMeshType( moab::Remapper::SourceMesh,
00771                                   static_cast< moab::TempestRemapper::TempestMeshType >( smetadata[0] ) );
00772         }
00773         // Rescale the radius of both to compute the intersection
00774         rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[0], radius_src );MB_CHK_ERR( rval );
00775         rval = remapper.ConvertMeshToTempest( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
00776         ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
00777 
00778         // Load the target mesh and validate
00779         rval = remapper.LoadNativeMesh( ctx.inFilenames[1], ctx.meshsets[1], tmetadata, additional_read_opts );MB_CHK_ERR( rval );
00780         if( tmetadata.size() )
00781         {
00782             // remapper.SetMeshType( moab::Remapper::TargetMesh,
00783             //                       static_cast< moab::TempestRemapper::TempestMeshType >( tmetadata[0] ), &tmetadata
00784             //                       );
00785             remapper.SetMeshType( moab::Remapper::TargetMesh,
00786                                   static_cast< moab::TempestRemapper::TempestMeshType >( tmetadata[0] ) );
00787         }
00788         rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[1], radius_dest );MB_CHK_ERR( rval );
00789         rval = remapper.ConvertMeshToTempest( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
00790         ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
00791     }
00792     else if( ctx.meshType == moab::TempestRemapper::ICO )
00793     {
00794         err = GenerateICOMesh( *tempest_mesh, ctx.blockSize, ctx.computeDual, ctx.outFilename, "NetCDF4" );
00795 
00796         if( err ) { rval = moab::MB_FAILURE; }
00797         else
00798         {
00799             ctx.meshes.push_back( tempest_mesh );
00800         }
00801     }
00802     else if( ctx.meshType == moab::TempestRemapper::RLL )
00803     {
00804         err = GenerateRLLMesh( *tempest_mesh,                     // Mesh& meshOut,
00805                                ctx.blockSize * 2, ctx.blockSize,  // int nLongitudes, int nLatitudes,
00806                                0.0, 360.0,                        // double dLonBegin, double dLonEnd,
00807                                -90.0, 90.0,                       // double dLatBegin, double dLatEnd,
00808                                false, false, false,  // bool fGlobalCap, bool fFlipLatLon, bool fForceGlobal,
00809                                "" /*ctx.inFilename*/, "",
00810                                "",  // std::string strInputFile, std::string strInputFileLonName, std::string
00811                                     // strInputFileLatName,
00812                                ctx.outFilename, "NetCDF4",  // std::string strOutputFile, std::string strOutputFormat
00813                                true                         // bool fVerbose
00814         );
00815 
00816         if( err ) { rval = moab::MB_FAILURE; }
00817         else
00818         {
00819             ctx.meshes.push_back( tempest_mesh );
00820         }
00821     }
00822     else  // default
00823     {
00824         err = GenerateCSMesh( *tempest_mesh, ctx.blockSize, ctx.outFilename, "NetCDF4" );
00825 
00826         if( err ) { rval = moab::MB_FAILURE; }
00827         else
00828         {
00829             ctx.meshes.push_back( tempest_mesh );
00830         }
00831     }
00832 
00833     if( ctx.meshType != moab::TempestRemapper::OVERLAP_MOAB && !tempest_mesh )
00834     {
00835         std::cout << "Tempest Mesh is not a complete object; Quitting...";
00836         exit( -1 );
00837     }
00838 
00839     return rval;
00840 }
00841 
00842 ///////////////////////////////////////////////
00843 //         Test functions
00844 
00845 double sample_slow_harmonic( double dLon, double dLat )
00846 {
00847     return ( 2.0 + cos( dLat ) * cos( dLat ) * cos( 2.0 * dLon ) );
00848 }
00849 
00850 double sample_fast_harmonic( double dLon, double dLat )
00851 {
00852     return ( 2.0 + pow( sin( 2.0 * dLat ), 16.0 ) * cos( 16.0 * dLon ) );
00853     // return (2.0 + pow(cos(2.0 * dLat), 16.0) * cos(16.0 * dLon));
00854 }
00855 
00856 double sample_constant( double /*dLon*/, double /*dLat*/ )
00857 {
00858     return 1.0;
00859 }
00860 
00861 double sample_stationary_vortex( double dLon, double dLat )
00862 {
00863     const double dLon0 = 0.0;
00864     const double dLat0 = 0.6;
00865     const double dR0   = 3.0;
00866     const double dD    = 5.0;
00867     const double dT    = 6.0;
00868 
00869     ///     Find the rotated longitude and latitude of a point on a sphere
00870     ///     with pole at (dLonC, dLatC).
00871     {
00872         double dSinC = sin( dLat0 );
00873         double dCosC = cos( dLat0 );
00874         double dCosT = cos( dLat );
00875         double dSinT = sin( dLat );
00876 
00877         double dTrm = dCosT * cos( dLon - dLon0 );
00878         double dX   = dSinC * dTrm - dCosC * dSinT;
00879         double dY   = dCosT * sin( dLon - dLon0 );
00880         double dZ   = dSinC * dSinT + dCosC * dTrm;
00881 
00882         dLon = atan2( dY, dX );
00883         if( dLon < 0.0 ) { dLon += 2.0 * M_PI; }
00884         dLat = asin( dZ );
00885     }
00886 
00887     double dRho = dR0 * cos( dLat );
00888     double dVt  = 3.0 * sqrt( 3.0 ) / 2.0 / cosh( dRho ) / cosh( dRho ) * tanh( dRho );
00889 
00890     double dOmega;
00891     if( dRho == 0.0 ) { dOmega = 0.0; }
00892     else
00893     {
00894         dOmega = dVt / dRho;
00895     }
00896 
00897     return ( 1.0 - tanh( dRho / dD * sin( dLon - dOmega * dT ) ) );
00898 }
00899 
00900 ///////////////////////////////////////////////
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines