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