Loading [MathJax]/extensions/tex2jax.js
Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TempestRemapper.hpp
Go to the documentation of this file.
1 /* 2  * ===================================================================================== 3  * 4  * Filename: TempestRemapper.hpp 5  * 6  * Description: Interface to the TempestRemap library to enable intersection and 7  * high-order conservative remapping of climate solution from 8  * arbitrary resolution of source and target grids on the sphere. 9  * 10  * Author: Vijay S. Mahadevan (vijaysm), mahadevan@anl.gov 11  * 12  * ===================================================================================== 13  */ 14  15 #ifndef MB_TEMPESTREMAPPER_HPP 16 #define MB_TEMPESTREMAPPER_HPP 17  18 #include "moab/Remapping/Remapper.hpp" 19 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp" 20 #include "moab/IntxMesh/IntxUtils.hpp" 21  22 // Tempest includes 23 #ifdef MOAB_HAVE_TEMPESTREMAP 24 #include "netcdfcpp.h" 25 #include "TempestRemapAPI.h" 26 #else 27 #error "This tool depends on TempestRemap library. Reconfigure using --with-tempestremap" 28 #endif 29  30 namespace moab 31 { 32  33 // Forward declare our friend, the mapper 34 class TempestOnlineMap; 35  36 class TempestRemapper : public Remapper 37 { 38  public: 39 #ifdef MOAB_HAVE_MPI 40  TempestRemapper( moab::Interface* mbInt, moab::ParallelComm* pcomm = NULL, bool offlineMode = false ) 41  : Remapper( mbInt, pcomm ), 42 #else 43  TempestRemapper( moab::Interface* mbInt, bool offlineMode = false ) 44  : Remapper( mbInt ), 45 #endif 46  offlineWorkflow( offlineMode ), meshValidate( false ), constructEdgeMap( false ), m_source_type( DEFAULT ), 47  m_target_type( DEFAULT ) 48  { 49  } 50  51  virtual ~TempestRemapper(); 52  53  // Mesh type with a correspondence to Tempest/Climate formats 54  enum TempestMeshType 55  { 56  DEFAULT = -1, 57  CS = 0, 58  RLL = 1, 59  ICO = 2, 60  ICOD = 3, 61  OVERLAP_FILES = 4, 62  OVERLAP_MEMORY = 5, 63  OVERLAP_MOAB = 6 64  }; 65  66  friend class TempestOnlineMap; 67  68  public: 69  /// <summary> 70  /// Initialize the TempestRemapper object internal datastructures including the mesh sets 71  /// and TempestRemap mesh references. 72  /// </summary> 73  virtual ErrorCode initialize( bool initialize_fsets = true ); 74  75  /// <summary> 76  /// Deallocate and clear any memory initialized in the TempestRemapper object 77  /// </summary> 78  virtual ErrorCode clear(); 79  80  /// <summary> 81  /// Generate a mesh in memory of given type (CS/RLL/ICO/MPAS(structured)) and store it 82  /// under the context specified by the user. 83  /// </summary> 84  moab::ErrorCode GenerateMesh( Remapper::IntersectionContext ctx, TempestMeshType type ); 85  86  /// <summary> 87  /// Load a mesh from disk of given type and store it under the context specified by the 88  /// user. 89  /// </summary> 90  moab::ErrorCode LoadMesh( Remapper::IntersectionContext ctx, std::string inputFilename, TempestMeshType type ); 91  92  /// <summary> 93  /// Construct a source covering mesh such that it completely encompasses the target grid in 94  /// parallel. This operation is critical to ensure that the parallel advancing-front 95  /// intersection algorithm can compute the intersection mesh only locally without any process 96  /// communication. 97  /// </summary> 98  moab::ErrorCode ConstructCoveringSet( double tolerance = 1e-8, 99  double radius_src = 1.0, 100  double radius_tgt = 1.0, 101  double boxeps = 0.1, 102  bool regional_mesh = false, 103  bool gnomonic = true, 104  int nb_ghost_layers = 0 ); 105  106  /// <summary> 107  /// Compute the intersection mesh between the source and target grids that have been 108  /// instantiated in the Remapper. This function invokes the parallel advancing-front 109  /// intersection algorithm internally for spherical meshes and can handle arbitrary 110  /// unstructured grids (CS, RLL, ICO, MPAS) with and without holes. 111  /// </summary> 112  moab::ErrorCode ComputeOverlapMesh( bool kdtree_search = true, bool use_tempest = false); 113  114  /* Converters between MOAB and Tempest representations */ 115  116  /// <summary> 117  /// Convert the TempestRemap mesh object to a corresponding MOAB mesh representation 118  /// according to the intersection context. 119  /// </summary> 120  moab::ErrorCode ConvertTempestMesh( Remapper::IntersectionContext ctx ); 121  122  /// <summary> 123  /// Convert the MOAB mesh representation to a corresponding TempestRemap mesh object 124  /// according to the intersection context. 125  /// </summary> 126  moab::ErrorCode ConvertMeshToTempest( Remapper::IntersectionContext ctx ); 127  128  /// <summary> 129  /// Get the TempestRemap mesh object according to the intersection context. 130  /// </summary> 131  Mesh* GetMesh( Remapper::IntersectionContext ctx ); 132  133  /// <summary> 134  /// Set the TempestRemap mesh object according to the intersection context. 135  /// </summary> 136  void SetMesh( Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite = true ); 137  138  void SetMeshSet( Remapper::IntersectionContext ctx /* Remapper::CoveringMesh*/, 139  moab::EntityHandle mset, 140  moab::Range* entities = nullptr ); 141  /// <summary> 142  /// Get the covering mesh (TempestRemap) object. 143  /// </summary> 144  Mesh* GetCoveringMesh(); 145  146  /// <summary> 147  /// Get the MOAB mesh set corresponding to the intersection context. 148  /// </summary> 149  moab::EntityHandle& GetMeshSet( Remapper::IntersectionContext ctx ); 150  151  /// <summary> 152  /// Const overload. Get the MOAB mesh set corresponding to the intersection context. 153  /// </summary> 154  moab::EntityHandle GetMeshSet( Remapper::IntersectionContext ctx ) const; 155  156  /// <summary> 157  /// Get the mesh element entities corresponding to the intersection context. 158  /// </summary> 159  moab::Range& GetMeshEntities( Remapper::IntersectionContext ctx ); 160  161  /// <summary> 162  /// Const overload. Get the mesh element entities corresponding to the intersection context. 163  /// </summary> 164  const moab::Range& GetMeshEntities( Remapper::IntersectionContext ctx ) const; 165  166  /// <summary> 167  /// Get the mesh vertices corresponding to the intersection context. Useful for point-cloud 168  /// meshes 169  /// </summary> 170  moab::Range& GetMeshVertices( Remapper::IntersectionContext ctx ); 171  172  /// <summary> 173  /// Const overload. Get the mesh vertices corresponding to the intersection context. Useful 174  /// for point-cloud meshes. 175  /// </summary> 176  const moab::Range& GetMeshVertices( Remapper::IntersectionContext ctx ) const; 177  178  /// <summary> 179  /// Get access to the underlying source covering set if available. Else return the source 180  /// set. 181  /// </summary> 182  moab::EntityHandle& GetCoveringSet(); 183  184  /// <summary> 185  /// Set the mesh type corresponding to the intersection context 186  /// </summary> 187  void SetMeshType( Remapper::IntersectionContext ctx, const std::vector< int >& metadata ); 188  189  /// <summary> 190  /// reconstruct mesh, used now only for IO; need a better solution maybe 191  /// </summary> 192  void ResetMeshSet( Remapper::IntersectionContext ctx, moab::EntityHandle meshSet ); 193  194  /// <summary> 195  /// Get the mesh type corresponding to the intersection context 196  /// </summary> 197  TempestMeshType GetMeshType( Remapper::IntersectionContext ctx ) const; 198  199  /// <summary> 200  /// Get the global ID corresponding to the local entity ID according to the context (source, 201  /// target, intersection) 202  /// </summary> 203  int GetGlobalID( Remapper::IntersectionContext ctx, int localID ); 204  205  /// <summary> 206  /// Get the local ID corresponding to the global entity ID according to the context (source, 207  /// target, intersection) 208  /// </summary> 209  int GetLocalID( Remapper::IntersectionContext ctx, int globalID ); 210  211  /// <summary> 212  /// Gather the overlap mesh and asssociated source/target data and write it out to disk 213  /// using the TempestRemap output interface. This information can then be used with the 214  /// "GenerateOfflineMap" tool in TempestRemap as needed. 215  /// </summary> 216  moab::ErrorCode WriteTempestIntersectionMesh( std::string strOutputFileName, 217  const bool fAllParallel, 218  const bool fInputConcave, 219  const bool fOutputConcave ); 220  221  /// <summary> 222  /// Generate the necessary metadata and specifically the GLL node numbering for DoFs for a 223  /// CS mesh. This negates the need for running external code like HOMME to output the 224  /// numbering needed for computing maps. The functionality is used through the `mbconvert` 225  /// tool to compute processor-invariant Global DoF IDs at GLL nodes. 226  /// </summary> 227  moab::ErrorCode GenerateCSMeshMetadata( const int ntot_elements, 228  moab::Range& entities, 229  moab::Range* secondary_entities, 230  const std::string& dofTagName, 231  int nP ); 232  233  /// <summary> 234  /// Generate the necessary metadata for DoF node numbering in a given mesh. 235  /// Currently, only the functionality to generate numbering on CS grids is supported. 236  /// </summary> 237  moab::ErrorCode GenerateMeshMetadata( Mesh& mesh, 238  const int ntot_elements, 239  moab::Range& entities, 240  moab::Range* secondary_entities, 241  const std::string dofTagName, 242  int nP ); 243  244  /// <summary> 245  /// Compute the local and global IDs for elements in source/target/coverage meshes. 246  /// </summary> 247  moab::ErrorCode ComputeGlobalLocalMaps(); 248  249  /// <summary> 250  /// Get all the ghosted overlap entities that were accumulated to enable conservation in 251  /// parallel 252  /// </summary> 253  moab::ErrorCode GetOverlapAugmentedEntities( moab::Range& sharedGhostEntities ); 254  255 #ifndef MOAB_HAVE_MPI 256  /// <summary> 257  /// Internal method to assign vertex and element global IDs if one does not exist already 258  /// </summary> 259  moab::ErrorCode assign_vertex_element_IDs( Tag idtag, 260  EntityHandle this_set, 261  const int dimension = 2, 262  const int start_id = 1 ); 263 #endif 264  // <summary> 265  /// Get the masks that could have been defined 266  /// </summary> 267  ErrorCode GetIMasks( Remapper::IntersectionContext ctx, std::vector< int >& masks ); 268  269  public: // public members 270  const bool offlineWorkflow; // check whether we are in an offline workflow context (mbtempest) 271  bool meshValidate; // Validate the mesh after loading from file 272  273  bool constructEdgeMap; // Construct the edge map within the TempestRemap datastructures 274  275  static const bool verbose = true; 276  277  private: 278  moab::ErrorCode convert_overlap_mesh_sorted_by_source(); 279  280  // private methods 281  moab::ErrorCode load_tempest_mesh_private( std::string inputFilename, Mesh** tempest_mesh ); 282  283  moab::ErrorCode convert_mesh_to_tempest_private( Mesh* mesh, 284  moab::EntityHandle meshset, 285  moab::Range& entities, 286  moab::Range* pverts ); 287  288  moab::ErrorCode convert_tempest_mesh_private( TempestMeshType type, 289  Mesh* mesh, 290  moab::EntityHandle& meshset, 291  moab::Range& entities, 292  moab::Range* vertices ); 293  294  moab::ErrorCode augment_overlap_set(); 295  296  /* Source meshset, mesh and entity references */ 297  Mesh* m_source; 298  TempestMeshType m_source_type; 299  moab::Range m_source_entities; 300  moab::Range m_source_vertices; 301  moab::EntityHandle m_source_set; 302  int max_source_edges; 303  bool point_cloud_source; 304  std::vector< int > m_source_metadata; 305  306  /* Target meshset, mesh and entity references */ 307  Mesh* m_target; 308  TempestMeshType m_target_type; 309  moab::Range m_target_entities; 310  moab::Range m_target_vertices; 311  moab::EntityHandle m_target_set; 312  int max_target_edges; 313  bool point_cloud_target; 314  std::vector< int > m_target_metadata; 315  316  /* Overlap meshset, mesh and entity references */ 317  Mesh* m_overlap; 318  TempestMeshType m_overlap_type; 319  moab::Range m_overlap_entities; 320  moab::EntityHandle m_overlap_set; 321  std::vector< std::pair< int, int > > m_sorted_overlap_order; 322  323  /* Intersection context on a sphere */ 324  moab::Intx2MeshOnSphere* mbintx; 325  326  /* Parallel - migrated mesh that is in the local view */ 327  Mesh* m_covering_source; 328  moab::EntityHandle m_covering_source_set; 329  moab::Range m_covering_source_entities; 330  moab::Range m_covering_source_vertices; 331  332  /* local to glboal and global to local ID maps */ 333  std::map< int, int > gid_to_lid_src, gid_to_lid_covsrc, gid_to_lid_tgt; 334  std::map< int, int > lid_to_gid_src, lid_to_gid_covsrc, lid_to_gid_tgt; 335  336  IntxAreaUtils::AreaMethod m_area_method; 337  338  bool rrmgrids; 339  bool is_parallel, is_root; 340  int rank, size; 341 }; 342  343 // Inline functions 344 inline Mesh* TempestRemapper::GetMesh( Remapper::IntersectionContext ctx ) 345 { 346  switch( ctx ) 347  { 348  case Remapper::SourceMesh: 349  return m_source; 350  case Remapper::TargetMesh: 351  return m_target; 352  case Remapper::OverlapMesh: 353  return m_overlap; 354  case Remapper::CoveringMesh: 355  return m_covering_source; 356  case Remapper::DEFAULT: 357  default: 358  return NULL; 359  } 360 } 361  362 inline void TempestRemapper::SetMesh( Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite ) 363 { 364  switch( ctx ) 365  { 366  case Remapper::SourceMesh: 367  if( !overwrite && m_source ) return; 368  if( overwrite && m_source ) delete m_source; 369  m_source = mesh; 370  break; 371  case Remapper::TargetMesh: 372  if( !overwrite && m_target ) return; 373  if( overwrite && m_target ) delete m_target; 374  m_target = mesh; 375  break; 376  case Remapper::OverlapMesh: 377  if( !overwrite && m_overlap ) return; 378  if( overwrite && m_overlap ) delete m_overlap; 379  m_overlap = mesh; 380  break; 381  case Remapper::CoveringMesh: 382  if( !overwrite && m_covering_source ) return; 383  if( overwrite && m_covering_source ) delete m_covering_source; 384  m_covering_source = mesh; 385  break; 386  case Remapper::DEFAULT: 387  default: 388  break; 389  } 390 } 391  392 inline void TempestRemapper::ResetMeshSet( Remapper::IntersectionContext ctx, moab::EntityHandle meshSet ) 393 { 394  switch( ctx ) 395  { 396  case Remapper::SourceMesh: 397  delete m_source; 398  m_source = new Mesh; 399  m_source_set = meshSet; 400  convert_mesh_to_tempest_private( m_source, m_source_set, m_source_entities, &m_source_vertices ); 401  m_source->CalculateFaceAreas( false ); // fInputConcave is false ? 402  break; 403  case Remapper::TargetMesh: 404  // not needed yet 405  break; 406  case Remapper::OverlapMesh: 407  // not needed yet 408  break; 409  case Remapper::CoveringMesh: 410  // not needed yet 411  break; 412  case Remapper::DEFAULT: 413  default: 414  break; 415  } 416 } 417  418  419 inline moab::EntityHandle& TempestRemapper::GetMeshSet( Remapper::IntersectionContext ctx ) 420 { 421  switch( ctx ) 422  { 423  case Remapper::SourceMesh: 424  return m_source_set; 425  case Remapper::TargetMesh: 426  return m_target_set; 427  case Remapper::OverlapMesh: 428  return m_overlap_set; 429  case Remapper::CoveringMesh: 430  return m_covering_source_set; 431  case Remapper::DEFAULT: 432  default: 433  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_set ); 434  } 435 } 436  437 inline moab::EntityHandle TempestRemapper::GetMeshSet( Remapper::IntersectionContext ctx ) const 438 { 439  switch( ctx ) 440  { 441  case Remapper::SourceMesh: 442  return m_source_set; 443  case Remapper::TargetMesh: 444  return m_target_set; 445  case Remapper::OverlapMesh: 446  return m_overlap_set; 447  case Remapper::CoveringMesh: 448  return m_covering_source_set; 449  case Remapper::DEFAULT: 450  default: 451  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_set ); 452  } 453 } 454  455 inline moab::Range& TempestRemapper::GetMeshEntities( Remapper::IntersectionContext ctx ) 456 { 457  switch( ctx ) 458  { 459  case Remapper::SourceMesh: 460  return m_source_entities; 461  case Remapper::TargetMesh: 462  return m_target_entities; 463  case Remapper::OverlapMesh: 464  return m_overlap_entities; 465  case Remapper::CoveringMesh: 466  return m_covering_source_entities; 467  case Remapper::DEFAULT: 468  default: 469  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_entities ); 470  } 471 } 472  473 inline const moab::Range& TempestRemapper::GetMeshEntities( Remapper::IntersectionContext ctx ) const 474 { 475  switch( ctx ) 476  { 477  case Remapper::SourceMesh: 478  return m_source_entities; 479  case Remapper::TargetMesh: 480  return m_target_entities; 481  case Remapper::OverlapMesh: 482  return m_overlap_entities; 483  case Remapper::CoveringMesh: 484  return m_covering_source_entities; 485  case Remapper::DEFAULT: 486  default: 487  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_entities ); 488  } 489 } 490  491 inline moab::Range& TempestRemapper::GetMeshVertices( Remapper::IntersectionContext ctx ) 492 { 493  switch( ctx ) 494  { 495  case Remapper::SourceMesh: 496  return m_source_vertices; 497  case Remapper::TargetMesh: 498  return m_target_vertices; 499  case Remapper::CoveringMesh: 500  return m_covering_source_vertices; 501  case Remapper::DEFAULT: 502  default: 503  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_source_vertices ); 504  } 505 } 506  507 inline const moab::Range& TempestRemapper::GetMeshVertices( Remapper::IntersectionContext ctx ) const 508 { 509  switch( ctx ) 510  { 511  case Remapper::SourceMesh: 512  return m_source_vertices; 513  case Remapper::TargetMesh: 514  return m_target_vertices; 515  case Remapper::CoveringMesh: 516  return m_covering_source_vertices; 517  case Remapper::DEFAULT: 518  default: 519  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_source_vertices ); 520  } 521 } 522  523 inline void TempestRemapper::SetMeshType( Remapper::IntersectionContext ctx, const std::vector< int >& metadata ) 524 { 525  switch( ctx ) 526  { 527  case Remapper::SourceMesh: 528  m_source_type = static_cast< moab::TempestRemapper::TempestMeshType >( metadata[0] ); 529  if( metadata[0] == 1 ) // RLL mesh 530  { 531  m_source_metadata.resize( 2 ); 532  m_source_metadata[0] = metadata[1]; 533  m_source_metadata[1] = metadata[2]; 534  } 535  else 536  { 537  m_source_metadata.resize( 1 ); 538  m_source_metadata[0] = metadata[1]; 539  } 540  break; 541  case Remapper::TargetMesh: 542  m_target_type = static_cast< moab::TempestRemapper::TempestMeshType >( metadata[0] ); 543  if( metadata[0] == 1 ) // RLL mesh 544  { 545  m_target_metadata.resize( 2 ); 546  m_target_metadata[0] = metadata[1]; 547  m_target_metadata[1] = metadata[2]; 548  } 549  else 550  { 551  m_target_metadata.resize( 1 ); 552  m_target_metadata[0] = metadata[1]; 553  } 554  break; 555  case Remapper::OverlapMesh: 556  m_overlap_type = moab::TempestRemapper::OVERLAP_FILES; 557  default: 558  break; 559  } 560 } 561  562 inline TempestRemapper::TempestMeshType TempestRemapper::GetMeshType( Remapper::IntersectionContext ctx ) const 563 { 564  switch( ctx ) 565  { 566  case Remapper::SourceMesh: 567  return m_source_type; 568  case Remapper::TargetMesh: 569  return m_target_type; 570  case Remapper::OverlapMesh: 571  return m_overlap_type; 572  case Remapper::DEFAULT: 573  default: 574  return TempestRemapper::DEFAULT; 575  } 576 } 577  578 inline Mesh* TempestRemapper::GetCoveringMesh() 579 { 580  return m_covering_source; 581 } 582  583 inline moab::EntityHandle& TempestRemapper::GetCoveringSet() 584 { 585  return m_covering_source_set; 586 } 587  588 inline int TempestRemapper::GetGlobalID( Remapper::IntersectionContext ctx, int localID ) 589 { 590  switch( ctx ) 591  { 592  case Remapper::SourceMesh: 593  return lid_to_gid_src[localID]; 594  case Remapper::TargetMesh: 595  return lid_to_gid_tgt[localID]; 596  case Remapper::CoveringMesh: 597  return lid_to_gid_covsrc[localID]; 598  case Remapper::OverlapMesh: 599  case Remapper::DEFAULT: 600  default: 601  return -1; 602  } 603 } 604  605 inline int TempestRemapper::GetLocalID( Remapper::IntersectionContext ctx, int globalID ) 606 { 607  switch( ctx ) 608  { 609  case Remapper::SourceMesh: 610  return gid_to_lid_src[globalID]; 611  case Remapper::TargetMesh: 612  return gid_to_lid_tgt[globalID]; 613  case Remapper::CoveringMesh: 614  return gid_to_lid_covsrc[globalID]; 615  case Remapper::DEFAULT: 616  case Remapper::OverlapMesh: 617  default: 618  return -1; 619  } 620 } 621  622 } // namespace moab 623  624 #endif // MB_TEMPESTREMAPPER_HPP