MOAB: Mesh Oriented datABase  (version 5.4.0)
TempestRemapper.hpp
Go to the documentation of this file.
00001 /*
00002  * =====================================================================================
00003  *
00004  *       Filename:  TempestRemapper.hpp
00005  *
00006  *    Description:  Interface to the TempestRemap library to enable intersection and
00007  *                  high-order conservative remapping of climate solution from
00008  *                  arbitrary resolution of source and target grids on the sphere.
00009  *
00010  *         Author:  Vijay S. Mahadevan (vijaysm), mahadevan@anl.gov
00011  *
00012  * =====================================================================================
00013  */
00014 
00015 #ifndef MB_TEMPESTREMAPPER_HPP
00016 #define MB_TEMPESTREMAPPER_HPP
00017 
00018 #include "moab/Remapping/Remapper.hpp"
00019 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
00020 #include "moab/IntxMesh/IntxUtils.hpp"
00021 
00022 // Tempest includes
00023 #ifdef MOAB_HAVE_TEMPESTREMAP
00024 #include "netcdfcpp.h"
00025 #include "TempestRemapAPI.h"
00026 #else
00027 #error "This tool depends on TempestRemap library. Reconfigure using --with-tempestremap"
00028 #endif
00029 
00030 namespace moab
00031 {
00032 
00033 // Forward declare our friend, the mapper
00034 class TempestOnlineMap;
00035 
00036 class TempestRemapper : public Remapper
00037 {
00038   public:
00039 #ifdef MOAB_HAVE_MPI
00040     TempestRemapper( moab::Interface* mbInt, moab::ParallelComm* pcomm = NULL )
00041         : Remapper( mbInt, pcomm ),
00042 #else
00043     TempestRemapper( moab::Interface* mbInt )
00044         : Remapper( mbInt ),
00045 #endif
00046           meshValidate( false ), constructEdgeMap( false ), m_source_type( DEFAULT ), m_target_type( DEFAULT )
00047     {
00048     }
00049 
00050     virtual ~TempestRemapper();
00051 
00052     // Mesh type with a correspondence to Tempest/Climate formats
00053     enum TempestMeshType
00054     {
00055         DEFAULT        = -1,
00056         CS             = 0,
00057         RLL            = 1,
00058         ICO            = 2,
00059         ICOD           = 3,
00060         OVERLAP_FILES  = 4,
00061         OVERLAP_MEMORY = 5,
00062         OVERLAP_MOAB   = 6
00063     };
00064 
00065     friend class TempestOnlineMap;
00066 
00067   public:
00068     /// <summary>
00069     ///     Initialize the TempestRemapper object internal datastructures including the mesh sets
00070     ///     and TempestRemap mesh references.
00071     /// </summary>
00072     virtual ErrorCode initialize( bool initialize_fsets = true );
00073 
00074     /// <summary>
00075     ///     Deallocate and clear any memory initialized in the TempestRemapper object
00076     /// </summary>
00077     virtual ErrorCode clear();
00078 
00079     /// <summary>
00080     ///     Generate a mesh in memory of given type (CS/RLL/ICO/MPAS(structured)) and store it
00081     ///     under the context specified by the user.
00082     /// </summary>
00083     moab::ErrorCode GenerateMesh( Remapper::IntersectionContext ctx, TempestMeshType type );
00084 
00085     /// <summary>
00086     ///     Load a mesh from disk of given type and store it under the context specified by the
00087     ///     user.
00088     /// </summary>
00089     moab::ErrorCode LoadMesh( Remapper::IntersectionContext ctx, std::string inputFilename, TempestMeshType type );
00090 
00091     /// <summary>
00092     ///     Construct a source covering mesh such that it completely encompasses the target grid in
00093     ///     parallel. This operation is critical to ensure that the parallel advancing-front
00094     ///     intersection algorithm can the intersection mesh only locally without any process
00095     ///     communication.
00096     /// </summary>
00097     moab::ErrorCode ConstructCoveringSet( double tolerance   = 1e-8,
00098                                           double radius_src  = 1.0,
00099                                           double radius_tgt  = 1.0,
00100                                           double boxeps      = 0.1,
00101                                           bool regional_mesh = false );
00102 
00103     /// <summary>
00104     ///     Compute the intersection mesh between the source and target grids that have been
00105     ///     instantiated in the Remapper. This function invokes the parallel advancing-front
00106     ///     intersection algorithm internally for spherical meshes and can handle arbitrary
00107     ///     unstructured grids (CS, RLL, ICO, MPAS) with and without holes.
00108     /// </summary>
00109     moab::ErrorCode ComputeOverlapMesh( bool kdtree_search = true, bool use_tempest = false );
00110 
00111     /* Converters between MOAB and Tempest representations */
00112 
00113     /// <summary>
00114     ///     Convert the TempestRemap mesh object to a corresponding MOAB mesh representation
00115     ///     according to the intersection context.
00116     /// </summary>
00117     moab::ErrorCode ConvertTempestMesh( Remapper::IntersectionContext ctx );
00118 
00119     /// <summary>
00120     ///     Convert the MOAB mesh representation to a corresponding TempestRemap mesh object
00121     ///     according to the intersection context.
00122     /// </summary>
00123     moab::ErrorCode ConvertMeshToTempest( Remapper::IntersectionContext ctx );
00124 
00125     /// <summary>
00126     ///     Get the TempestRemap mesh object according to the intersection context.
00127     /// </summary>
00128     Mesh* GetMesh( Remapper::IntersectionContext ctx );
00129 
00130     /// <summary>
00131     ///     Set the TempestRemap mesh object according to the intersection context.
00132     /// </summary>
00133     void SetMesh( Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite = true );
00134 
00135     void SetMeshSet( Remapper::IntersectionContext ctx /* Remapper::CoveringMesh*/,
00136                      moab::EntityHandle mset,
00137                      moab::Range& entities );
00138     /// <summary>
00139     ///     Get the covering mesh (TempestRemap) object.
00140     /// </summary>
00141     Mesh* GetCoveringMesh();
00142 
00143     /// <summary>
00144     ///     Get the MOAB mesh set corresponding to the intersection context.
00145     /// </summary>
00146     moab::EntityHandle& GetMeshSet( Remapper::IntersectionContext ctx );
00147 
00148     /// <summary>
00149     ///     Const overload. Get the MOAB mesh set corresponding to the intersection context.
00150     /// </summary>
00151     moab::EntityHandle GetMeshSet( Remapper::IntersectionContext ctx ) const;
00152 
00153     /// <summary>
00154     ///     Get the mesh element entities corresponding to the intersection context.
00155     /// </summary>
00156     moab::Range& GetMeshEntities( Remapper::IntersectionContext ctx );
00157 
00158     /// <summary>
00159     ///     Const overload. Get the mesh element entities corresponding to the intersection context.
00160     /// </summary>
00161     const moab::Range& GetMeshEntities( Remapper::IntersectionContext ctx ) const;
00162 
00163     /// <summary>
00164     ///     Get the mesh vertices corresponding to the intersection context. Useful for point-cloud
00165     ///     meshes
00166     /// </summary>
00167     moab::Range& GetMeshVertices( Remapper::IntersectionContext ctx );
00168 
00169     /// <summary>
00170     ///     Const overload. Get the mesh vertices corresponding to the intersection context. Useful
00171     ///     for point-cloud meshes.
00172     /// </summary>
00173     const moab::Range& GetMeshVertices( Remapper::IntersectionContext ctx ) const;
00174 
00175     /// <summary>
00176     ///     Get access to the underlying source covering set if available. Else return the source
00177     ///     set.
00178     /// </summary>
00179     moab::EntityHandle& GetCoveringSet();
00180 
00181     /// <summary>
00182     ///     Set the mesh type corresponding to the intersection context
00183     /// </summary>
00184     void SetMeshType( Remapper::IntersectionContext ctx,
00185                       TempestMeshType type,
00186                       const std::vector< int >* metadata = nullptr );
00187 
00188     /// <summary>
00189     ///     Get the mesh type corresponding to the intersection context
00190     /// </summary>
00191     TempestMeshType GetMeshType( Remapper::IntersectionContext ctx ) const;
00192 
00193     /// <summary>
00194     ///     Get the global ID corresponding to the local entity ID according to the context (source,
00195     ///     target, intersection)
00196     /// </summary>
00197     int GetGlobalID( Remapper::IntersectionContext ctx, int localID );
00198 
00199     /// <summary>
00200     ///     Get the local ID corresponding to the global entity ID according to the context (source,
00201     ///     target, intersection)
00202     /// </summary>
00203     int GetLocalID( Remapper::IntersectionContext ctx, int globalID );
00204 
00205     /// <summary>
00206     ///     Gather the overlap mesh and asssociated source/target data and write it out to disk
00207     ///     using the TempestRemap output interface. This information can then be used with the
00208     ///     "GenerateOfflineMap" tool in TempestRemap as needed.
00209     /// </summary>
00210     moab::ErrorCode WriteTempestIntersectionMesh( std::string strOutputFileName,
00211                                                   const bool fAllParallel,
00212                                                   const bool fInputConcave,
00213                                                   const bool fOutputConcave );
00214 
00215     /// <summary>
00216     ///     Generate the necessary metadata and specifically the GLL node numbering for DoFs for a
00217     ///     CS mesh. This negates the need for running external code like HOMME to output the
00218     ///     numbering needed for computing maps. The functionality is used through the `mbconvert`
00219     ///     tool to compute processor-invariant Global DoF IDs at GLL nodes.
00220     /// </summary>
00221     moab::ErrorCode GenerateCSMeshMetadata( const int ntot_elements,
00222                                             moab::Range& entities,
00223                                             moab::Range* secondary_entities,
00224                                             const std::string& dofTagName,
00225                                             int nP );
00226 
00227     /// <summary>
00228     ///     Generate the necessary metadata for DoF node numbering in a given mesh.
00229     ///     Currently, only the functionality to generate numbering on CS grids is supported.
00230     /// </summary>
00231     moab::ErrorCode GenerateMeshMetadata( Mesh& mesh,
00232                                           const int ntot_elements,
00233                                           moab::Range& entities,
00234                                           moab::Range* secondary_entities,
00235                                           const std::string dofTagName,
00236                                           int nP );
00237 
00238     /// <summary>
00239     ///     Compute the local and global IDs for elements in source/target/coverage meshes.
00240     /// </summary>
00241     moab::ErrorCode ComputeGlobalLocalMaps();
00242 
00243     /// <summary>
00244     ///     Get all the ghosted overlap entities that were accumulated to enable conservation in
00245     /// parallel
00246     /// </summary>
00247     moab::ErrorCode GetOverlapAugmentedEntities( moab::Range& sharedGhostEntities );
00248 
00249 #ifndef MOAB_HAVE_MPI
00250     /// <summary>
00251     ///    Internal method to assign vertex and element global IDs if one does not exist already
00252     /// </summary>
00253     moab::ErrorCode assign_vertex_element_IDs( Tag idtag,
00254                                                EntityHandle this_set,
00255                                                const int dimension = 2,
00256                                                const int start_id  = 1 );
00257 #endif
00258     //  <summary>
00259     ///     Get the masks that could have been defined
00260     /// </summary>
00261     ErrorCode GetIMasks( Remapper::IntersectionContext ctx, std::vector< int >& masks );
00262 
00263   public:               // public members
00264     bool meshValidate;  // Validate the mesh after loading from file
00265 
00266     bool constructEdgeMap;  //  Construct the edge map within the TempestRemap datastructures
00267 
00268     static const bool verbose = true;
00269 
00270   private:
00271     moab::ErrorCode convert_overlap_mesh_sorted_by_source();
00272 
00273     // private methods
00274     moab::ErrorCode load_tempest_mesh_private( std::string inputFilename, Mesh** tempest_mesh );
00275 
00276     moab::ErrorCode convert_mesh_to_tempest_private( Mesh* mesh,
00277                                                      moab::EntityHandle meshset,
00278                                                      moab::Range& entities,
00279                                                      moab::Range* pverts );
00280 
00281     moab::ErrorCode convert_tempest_mesh_private( TempestMeshType type,
00282                                                   Mesh* mesh,
00283                                                   moab::EntityHandle& meshset,
00284                                                   moab::Range& entities,
00285                                                   moab::Range* vertices );
00286 
00287     moab::ErrorCode augment_overlap_set();
00288 
00289     /* Source meshset, mesh and entity references */
00290     Mesh* m_source;
00291     TempestMeshType m_source_type;
00292     moab::Range m_source_entities;
00293     moab::Range m_source_vertices;
00294     moab::EntityHandle m_source_set;
00295     int max_source_edges;
00296     bool point_cloud_source;
00297     std::vector< int > m_source_metadata;
00298 
00299     /* Target meshset, mesh and entity references */
00300     Mesh* m_target;
00301     TempestMeshType m_target_type;
00302     moab::Range m_target_entities;
00303     moab::Range m_target_vertices;
00304     moab::EntityHandle m_target_set;
00305     int max_target_edges;
00306     bool point_cloud_target;
00307     std::vector< int > m_target_metadata;
00308 
00309     /* Overlap meshset, mesh and entity references */
00310     Mesh* m_overlap;
00311     TempestMeshType m_overlap_type;
00312     moab::Range m_overlap_entities;
00313     moab::EntityHandle m_overlap_set;
00314     std::vector< std::pair< int, int > > m_sorted_overlap_order;
00315 
00316     /* Intersection context on a sphere */
00317     moab::Intx2MeshOnSphere* mbintx;
00318 
00319     /* Parallel - migrated mesh that is in the local view */
00320     Mesh* m_covering_source;
00321     moab::EntityHandle m_covering_source_set;
00322     moab::Range m_covering_source_entities;
00323     moab::Range m_covering_source_vertices;
00324 
00325     /* local to glboal and global to local ID maps */
00326     std::map< int, int > gid_to_lid_src, gid_to_lid_covsrc, gid_to_lid_tgt;
00327     std::map< int, int > lid_to_gid_src, lid_to_gid_covsrc, lid_to_gid_tgt;
00328 
00329     IntxAreaUtils::AreaMethod m_area_method;
00330 
00331     bool rrmgrids;
00332     bool is_parallel, is_root;
00333     int rank, size;
00334 };
00335 
00336 // Inline functions
00337 inline Mesh* TempestRemapper::GetMesh( Remapper::IntersectionContext ctx )
00338 {
00339     switch( ctx )
00340     {
00341         case Remapper::SourceMesh:
00342             return m_source;
00343         case Remapper::TargetMesh:
00344             return m_target;
00345         case Remapper::OverlapMesh:
00346             return m_overlap;
00347         case Remapper::CoveringMesh:
00348             return m_covering_source;
00349         case Remapper::DEFAULT:
00350         default:
00351             return NULL;
00352     }
00353 }
00354 
00355 inline void TempestRemapper::SetMesh( Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite )
00356 {
00357     switch( ctx )
00358     {
00359         case Remapper::SourceMesh:
00360             if( !overwrite && m_source ) return;
00361             if( overwrite && m_source ) delete m_source;
00362             m_source = mesh;
00363             break;
00364         case Remapper::TargetMesh:
00365             if( !overwrite && m_target ) return;
00366             if( overwrite && m_target ) delete m_target;
00367             m_target = mesh;
00368             break;
00369         case Remapper::OverlapMesh:
00370             if( !overwrite && m_overlap ) return;
00371             if( overwrite && m_overlap ) delete m_overlap;
00372             m_overlap = mesh;
00373             break;
00374         case Remapper::CoveringMesh:
00375             if( !overwrite && m_covering_source ) return;
00376             if( overwrite && m_covering_source ) delete m_covering_source;
00377             m_covering_source = mesh;
00378             break;
00379         case Remapper::DEFAULT:
00380         default:
00381             break;
00382     }
00383 }
00384 
00385 inline moab::EntityHandle& TempestRemapper::GetMeshSet( Remapper::IntersectionContext ctx )
00386 {
00387     switch( ctx )
00388     {
00389         case Remapper::SourceMesh:
00390             return m_source_set;
00391         case Remapper::TargetMesh:
00392             return m_target_set;
00393         case Remapper::OverlapMesh:
00394             return m_overlap_set;
00395         case Remapper::CoveringMesh:
00396             return m_covering_source_set;
00397         case Remapper::DEFAULT:
00398         default:
00399             MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_set );
00400     }
00401 }
00402 
00403 inline moab::EntityHandle TempestRemapper::GetMeshSet( Remapper::IntersectionContext ctx ) const
00404 {
00405     switch( ctx )
00406     {
00407         case Remapper::SourceMesh:
00408             return m_source_set;
00409         case Remapper::TargetMesh:
00410             return m_target_set;
00411         case Remapper::OverlapMesh:
00412             return m_overlap_set;
00413         case Remapper::CoveringMesh:
00414             return m_covering_source_set;
00415         case Remapper::DEFAULT:
00416         default:
00417             MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_set );
00418     }
00419 }
00420 
00421 inline moab::Range& TempestRemapper::GetMeshEntities( Remapper::IntersectionContext ctx )
00422 {
00423     switch( ctx )
00424     {
00425         case Remapper::SourceMesh:
00426             return m_source_entities;
00427         case Remapper::TargetMesh:
00428             return m_target_entities;
00429         case Remapper::OverlapMesh:
00430             return m_overlap_entities;
00431         case Remapper::CoveringMesh:
00432             return m_covering_source_entities;
00433         case Remapper::DEFAULT:
00434         default:
00435             MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_entities );
00436     }
00437 }
00438 
00439 inline const moab::Range& TempestRemapper::GetMeshEntities( Remapper::IntersectionContext ctx ) const
00440 {
00441     switch( ctx )
00442     {
00443         case Remapper::SourceMesh:
00444             return m_source_entities;
00445         case Remapper::TargetMesh:
00446             return m_target_entities;
00447         case Remapper::OverlapMesh:
00448             return m_overlap_entities;
00449         case Remapper::CoveringMesh:
00450             return m_covering_source_entities;
00451         case Remapper::DEFAULT:
00452         default:
00453             MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_entities );
00454     }
00455 }
00456 
00457 inline moab::Range& TempestRemapper::GetMeshVertices( Remapper::IntersectionContext ctx )
00458 {
00459     switch( ctx )
00460     {
00461         case Remapper::SourceMesh:
00462             return m_source_vertices;
00463         case Remapper::TargetMesh:
00464             return m_target_vertices;
00465         case Remapper::CoveringMesh:
00466             return m_covering_source_vertices;
00467         case Remapper::DEFAULT:
00468         default:
00469             MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_source_vertices );
00470     }
00471 }
00472 
00473 inline const moab::Range& TempestRemapper::GetMeshVertices( Remapper::IntersectionContext ctx ) const
00474 {
00475     switch( ctx )
00476     {
00477         case Remapper::SourceMesh:
00478             return m_source_vertices;
00479         case Remapper::TargetMesh:
00480             return m_target_vertices;
00481         case Remapper::CoveringMesh:
00482             return m_covering_source_vertices;
00483         case Remapper::DEFAULT:
00484         default:
00485             MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_source_vertices );
00486     }
00487 }
00488 
00489 inline void TempestRemapper::SetMeshType( Remapper::IntersectionContext ctx,
00490                                           TempestRemapper::TempestMeshType type,
00491                                           const std::vector< int >* metadata )
00492 {
00493     switch( ctx )
00494     {
00495         case Remapper::SourceMesh:
00496             m_source_type = type;
00497             if( metadata )
00498             {
00499                 m_source_metadata.resize( metadata->size() );
00500                 std::copy( metadata->begin(), metadata->end(), m_source_metadata.begin() );
00501             }
00502             break;
00503         case Remapper::TargetMesh:
00504             m_target_type = type;
00505             if( metadata )
00506             {
00507                 m_target_metadata.resize( metadata->size() );
00508                 std::copy( metadata->begin(), metadata->end(), m_target_metadata.begin() );
00509             }
00510             break;
00511         case Remapper::OverlapMesh:
00512             m_overlap_type = type;
00513             break;
00514         case Remapper::DEFAULT:
00515         default:
00516             break;
00517     }
00518 }
00519 
00520 inline TempestRemapper::TempestMeshType TempestRemapper::GetMeshType( Remapper::IntersectionContext ctx ) const
00521 {
00522     switch( ctx )
00523     {
00524         case Remapper::SourceMesh:
00525             return m_source_type;
00526         case Remapper::TargetMesh:
00527             return m_target_type;
00528         case Remapper::OverlapMesh:
00529             return m_overlap_type;
00530         case Remapper::DEFAULT:
00531         default:
00532             return TempestRemapper::DEFAULT;
00533     }
00534 }
00535 
00536 inline Mesh* TempestRemapper::GetCoveringMesh()
00537 {
00538     return m_covering_source;
00539 }
00540 
00541 inline moab::EntityHandle& TempestRemapper::GetCoveringSet()
00542 {
00543     return m_covering_source_set;
00544 }
00545 
00546 inline int TempestRemapper::GetGlobalID( Remapper::IntersectionContext ctx, int localID )
00547 {
00548     switch( ctx )
00549     {
00550         case Remapper::SourceMesh:
00551             return lid_to_gid_src[localID];
00552         case Remapper::TargetMesh:
00553             return lid_to_gid_tgt[localID];
00554         case Remapper::CoveringMesh:
00555             return lid_to_gid_covsrc[localID];
00556         case Remapper::OverlapMesh:
00557         case Remapper::DEFAULT:
00558         default:
00559             return -1;
00560     }
00561 }
00562 
00563 inline int TempestRemapper::GetLocalID( Remapper::IntersectionContext ctx, int globalID )
00564 {
00565     switch( ctx )
00566     {
00567         case Remapper::SourceMesh:
00568             return gid_to_lid_src[globalID];
00569         case Remapper::TargetMesh:
00570             return gid_to_lid_tgt[globalID];
00571         case Remapper::CoveringMesh:
00572             return gid_to_lid_covsrc[globalID];
00573         case Remapper::DEFAULT:
00574         case Remapper::OverlapMesh:
00575         default:
00576             return -1;
00577     }
00578 }
00579 
00580 }  // namespace moab
00581 
00582 #endif  // MB_TEMPESTREMAPPER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines