Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
TempestRemapper.cpp
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), [email protected]
11  *
12  * =====================================================================================
13  */
14 
15 #include <string>
16 #include <iostream>
17 #include <cassert>
18 #include <array>
19 #include <numeric> // std::iota
20 #include <algorithm> // std::sort, std::stable_sort
21 #include <unordered_map> // global-id -> local-index lookup
22 
23 #include "DebugOutput.hpp"
24 #include "moab/Remapping/TempestRemapper.hpp"
25 #include "moab/ReadUtilIface.hpp"
26 // needed for higher order mapping, retrieve additional layers of cells with bridge methods
27 #include "moab/MeshTopoUtil.hpp"
28 #include "AEntityFactory.hpp"
29 
30 // Intersection includes
31 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
32 #include "moab/IntxMesh/IntxUtils.hpp"
33 
34 #include "moab/AdaptiveKDTree.hpp"
35 #include "moab/SpatialLocator.hpp"
36 
37 // skinner for augmenting overlap mesh to complete coverage
38 #include "moab/Skinner.hpp"
39 #include "MBParallelConventions.h"
40 
41 #ifdef MOAB_HAVE_TEMPESTREMAP
42 #include "Announce.h"
43 #include "FiniteElementTools.h"
44 #include "GaussLobattoQuadrature.h"
45 #endif
46 
47 // #define VERBOSE
48 
49 namespace moab
50 {
51 
52 ///////////////////////////////////////////////////////////////////////////////////
53 
54 ErrorCode TempestRemapper::initialize( bool initialize_fsets )
55 {
56  if( initialize_fsets )
57  {
58  MB_CHK_SET_ERR( m_interface->create_meshset( moab::MESHSET_SET, m_source_set ), "Can't create new set" );
59  MB_CHK_SET_ERR( m_interface->create_meshset( moab::MESHSET_SET, m_target_set ), "Can't create new set" );
60  MB_CHK_SET_ERR( m_interface->create_meshset( moab::MESHSET_SET, m_overlap_set ), "Can't create new set" );
61  }
62  else
63  {
64  m_source_set = 0;
65  m_target_set = 0;
66  m_overlap_set = 0;
67  }
68 
69  is_parallel = false;
70  is_root = true;
71  rank = 0;
72  size = 1;
73 #ifdef MOAB_HAVE_MPI
74  int flagInit;
75  MPI_Initialized( &flagInit );
76  if( flagInit )
77  {
78  assert( m_pcomm != nullptr );
79  rank = m_pcomm->rank();
80  size = m_pcomm->size();
81  is_root = ( rank == 0 );
82  is_parallel = ( size > 1 );
83  // is_parallel = true;
84  }
85  AnnounceOnlyOutputOnRankZero();
86 #endif
87 
88  m_source = nullptr;
89  m_target = nullptr;
90  m_overlap = nullptr;
91  m_covering_source = nullptr;
92 
93  point_cloud_source = false;
94  point_cloud_target = false;
95 
96  return MB_SUCCESS;
97 }
98 
99 ///////////////////////////////////////////////////////////////////////////////////
100 
101 TempestRemapper::~TempestRemapper()
102 {
103  this->clear();
104 }
105 
106 ErrorCode TempestRemapper::clear()
107 {
108  // destroy all meshes
109  if( m_source )
110  {
111  delete m_source;
112  m_source = nullptr;
113  }
114  if( m_target )
115  {
116  delete m_target;
117  m_target = nullptr;
118  }
119  if( m_overlap )
120  {
121  delete m_overlap;
122  m_overlap = nullptr;
123  }
124  if( m_covering_source && size > 1 )
125  {
126  delete m_covering_source;
127  m_covering_source = nullptr;
128  }
129 
130  point_cloud_source = false;
131  point_cloud_target = false;
132 
133  m_source_entities.clear();
134  m_source_vertices.clear();
135  m_covering_source_entities.clear();
136  m_covering_source_vertices.clear();
137  m_target_entities.clear();
138  m_target_vertices.clear();
139  m_overlap_entities.clear();
140 
141  return MB_SUCCESS;
142 }
143 
144 ///////////////////////////////////////////////////////////////////////////////////
145 #ifdef MOAB_HAVE_NETCDF
146 
147 ErrorCode TempestRemapper::LoadMesh( Remapper::IntersectionContext ctx,
148  std::string inputFilename,
149  TempestMeshType type )
150 {
151  if( ctx == Remapper::SourceMesh )
152  {
153  m_source_type = type;
154  return load_tempest_mesh_private( inputFilename, &m_source );
155  }
156  else if( ctx == Remapper::TargetMesh )
157  {
158  m_target_type = type;
159  return load_tempest_mesh_private( inputFilename, &m_target );
160  }
161  else if( ctx != Remapper::DEFAULT )
162  {
163  m_overlap_type = type;
164  return load_tempest_mesh_private( inputFilename, &m_overlap );
165  }
166  else
167  {
168  MB_CHK_SET_ERR( MB_FAILURE, "Invalid IntersectionContext context provided" );
169  }
170 }
171 
172 ErrorCode TempestRemapper::load_tempest_mesh_private( std::string inputFilename, Mesh** tempest_mesh )
173 {
174  const bool outputEnabled = ( TempestRemapper::verbose && is_root );
175  if( outputEnabled ) std::cout << "\nLoading TempestRemap Mesh object from file = " << inputFilename << " ...\n";
176 
177  {
178  // NcError is a TempestRemap/netcdfcpp type; this whole routine is already inside the
179  // enclosing #ifdef MOAB_HAVE_NETCDF (LoadMesh forwards to TempestRemap's file reader).
180  NcError error( NcError::silent_nonfatal );
181 
182  try
183  {
184  // Load input mesh
185  if( outputEnabled ) std::cout << "Loading mesh ...\n";
186  Mesh* mesh = new Mesh( inputFilename );
187  mesh->RemoveZeroEdges();
188  if( outputEnabled ) std::cout << "----------------\n";
189 
190  // Validate mesh
191  if( meshValidate )
192  {
193  if( outputEnabled ) std::cout << "Validating mesh ...\n";
194  mesh->Validate();
195  if( outputEnabled ) std::cout << "-------------------\n";
196  }
197 
198  // Construct the edge map on the mesh
199  if( constructEdgeMap )
200  {
201  if( outputEnabled ) std::cout << "Constructing edge map on mesh ...\n";
202  mesh->ConstructEdgeMap( false );
203  if( outputEnabled ) std::cout << "---------------------------------\n";
204  }
205 
206  if( tempest_mesh ) *tempest_mesh = mesh;
207  }
208  catch( Exception& e )
209  {
210  std::cout << "TempestRemap ERROR: " << e.ToString() << "\n";
211  return MB_FAILURE;
212  }
213  catch( ... )
214  {
215  return MB_FAILURE;
216  }
217  }
218  return MB_SUCCESS;
219 }
220 #endif
221 ///////////////////////////////////////////////////////////////////////////////////
222 
223 ErrorCode TempestRemapper::ConvertTempestMesh( Remapper::IntersectionContext ctx )
224 {
225  const bool outputEnabled = ( TempestRemapper::verbose && is_root );
226  if( ctx == Remapper::SourceMesh )
227  {
228  if( outputEnabled ) std::cout << "Converting (source) TempestRemap Mesh object to MOAB representation ...\n";
229  return convert_tempest_mesh_private( m_source_type, m_source, m_source_set, m_source_entities,
230  &m_source_vertices );
231  }
232  else if( ctx == Remapper::TargetMesh )
233  {
234  if( outputEnabled ) std::cout << "Converting (target) TempestRemap Mesh object to MOAB representation ...\n";
235  return convert_tempest_mesh_private( m_target_type, m_target, m_target_set, m_target_entities,
236  &m_target_vertices );
237  }
238  else if( ctx != Remapper::DEFAULT )
239  {
240  if( outputEnabled ) std::cout << "Converting (overlap) TempestRemap Mesh object to MOAB representation ...\n";
241  return convert_tempest_mesh_private( m_overlap_type, m_overlap, m_overlap_set, m_overlap_entities, nullptr );
242  }
243  else
244  {
245  MB_CHK_SET_ERR( MB_FAILURE, "Invalid IntersectionContext context provided" );
246  }
247 }
248 
249 #define NEW_CONVERT_LOGIC
250 
251 #ifdef NEW_CONVERT_LOGIC
252 ErrorCode TempestRemapper::convert_tempest_mesh_private( TempestMeshType /*meshType*/,
253  Mesh* mesh,
254  EntityHandle& mesh_set,
255  Range& entities,
256  Range* vertices )
257 {
258  const bool outputEnabled = ( TempestRemapper::verbose && is_root );
259  const NodeVector& nodes = mesh->nodes;
260  const FaceVector& faces = mesh->faces;
261 
262  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
263  dbgprint.set_prefix( "[TempestToMOAB]: " );
264 
265  ReadUtilIface* iface;
266  MB_CHK_SET_ERR( m_interface->query_interface( iface ), "Can't get reader interface" );
267 
268  Tag gidTag = m_interface->globalId_tag();
269 
270  // Set the data for the vertices
271  std::vector< double* > arrays;
272  std::vector< int > gidsv( nodes.size() );
273  EntityHandle startv;
274  MB_CHK_SET_ERR( iface->get_node_coords( 3, nodes.size(), 0, startv, arrays ), "Can't get node coords" );
275  for( unsigned iverts = 0; iverts < nodes.size(); ++iverts )
276  {
277  const Node& node = nodes[iverts];
278  arrays[0][iverts] = node.x;
279  arrays[1][iverts] = node.y;
280  arrays[2][iverts] = node.z;
281  gidsv[iverts] = iverts + 1;
282  }
283  Range mbverts( startv, startv + nodes.size() - 1 );
284  MB_CHK_SET_ERR( m_interface->add_entities( mesh_set, mbverts ), "Can't add entities" );
285  MB_CHK_SET_ERR( m_interface->tag_set_data( gidTag, mbverts, &gidsv[0] ), "Can't set global_id tag" );
286 
287  gidsv.clear();
288  entities.clear();
289 
290  Tag srcParentTag, tgtParentTag;
291  std::vector< int > srcParent( faces.size(), -1 ), tgtParent( faces.size(), -1 );
292  std::vector< int > gidse( faces.size(), -1 );
293  bool storeParentInfo = ( mesh->vecSourceFaceIx.size() > 0 );
294 
295  if( storeParentInfo )
296  {
297  int defaultInt = -1;
298  MB_CHK_SET_ERR( m_interface->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag,
299  MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt ),
300  "can't create positive tag" );
301 
302  MB_CHK_SET_ERR( m_interface->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag,
303  MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt ),
304  "can't create negative tag" );
305  }
306 
307  // Let us first perform a full pass assuming arbitrary polygons. This is especially true for
308  // overlap meshes.
309  // 1. We do a first pass over faces, decipher edge size and group into categories based on
310  // element type
311  // 2. Next we loop over type, and add blocks of elements into MOAB
312  // 3. For each block within the loop, also update the connectivity of elements.
313  {
314  if( outputEnabled )
315  dbgprint.printf( 0, "..Mesh size: Nodes [%zu] Elements [%zu].\n", nodes.size(), faces.size() );
316  std::vector< EntityHandle > mbcells( faces.size() );
317  unsigned ntris = 0, nquads = 0, npolys = 0;
318  std::vector< EntityHandle > conn( 16 );
319 
320  for( unsigned ifaces = 0; ifaces < faces.size(); ++ifaces )
321  {
322  const Face& face = faces[ifaces];
323  const unsigned num_v_per_elem = face.edges.size();
324 
325  for( unsigned iedges = 0; iedges < num_v_per_elem; ++iedges )
326  {
327  conn[iedges] = startv + face.edges[iedges].node[0];
328  }
329 
330  switch( num_v_per_elem )
331  {
332  case 3:
333  // if( outputEnabled )
334  // dbgprint.printf( 0, "....Block %d: Triangular Elements [%u].\n", iBlock++, nPolys[iType] );
335  MB_CHK_SET_ERR( m_interface->create_element( MBTRI, &conn[0], num_v_per_elem, mbcells[ifaces] ),
336  "Can't get element connectivity" );
337  ntris++;
338  break;
339  case 4:
340  // if( outputEnabled )
341  // dbgprint.printf( 0, "....Block %d: Quadrilateral Elements [%u].\n", iBlock++, nPolys[iType] );
342  MB_CHK_SET_ERR( m_interface->create_element( MBQUAD, &conn[0], num_v_per_elem, mbcells[ifaces] ),
343  "Can't get element connectivity" );
344  nquads++;
345  break;
346  default:
347  // if( outputEnabled )
348  // dbgprint.printf( 0, "....Block %d: Polygonal [%u] Elements [%u].\n", iBlock++, iType,
349  // nPolys[iType] );
350  MB_CHK_SET_ERR( m_interface->create_element( MBPOLYGON, &conn[0], num_v_per_elem, mbcells[ifaces] ),
351  "Can't get element connectivity" );
352  npolys++;
353  break;
354  }
355 
356  gidse[ifaces] = ifaces + 1;
357 
358  if( storeParentInfo )
359  {
360  srcParent[ifaces] = mesh->vecSourceFaceIx[ifaces] + 1;
361  tgtParent[ifaces] = mesh->vecTargetFaceIx[ifaces] + 1;
362  }
363  }
364 
365  if( ntris ) dbgprint.printf( 0, "....Triangular Elements [%u].\n", ntris );
366  if( nquads ) dbgprint.printf( 0, "....Quadrangular Elements [%u].\n", nquads );
367  if( npolys ) dbgprint.printf( 0, "....Polygonal Elements [%u].\n", npolys );
368 
369  MB_CHK_SET_ERR( m_interface->add_entities( mesh_set, &mbcells[0], mbcells.size() ), "Could not add entities" );
370 
371  MB_CHK_SET_ERR( m_interface->tag_set_data( gidTag, &mbcells[0], mbcells.size(), &gidse[0] ),
372  "Can't set global_id tag" );
373 #ifdef MOAB_HAVE_MPI
374  MB_CHK_SET_ERR( m_pcomm->assign_global_ids(mesh_set, 2, 1, false, true, false ), "Unable to set global IDs" );
375 #endif
376 
377  if( storeParentInfo )
378  {
379  MB_CHK_SET_ERR( m_interface->tag_set_data( srcParentTag, &mbcells[0], mbcells.size(), &srcParent[0] ),
380  "Can't set tag data" );
381  MB_CHK_SET_ERR( m_interface->tag_set_data( tgtParentTag, &mbcells[0], mbcells.size(), &tgtParent[0] ),
382  "Can't set tag data" );
383  }
384 
385  // insert from mbcells to entities to preserve ordering
386  std::copy( mbcells.begin(), mbcells.end(), range_inserter( entities ) );
387  }
388 
389  if( vertices ) *vertices = mbverts;
390 
391  return MB_SUCCESS;
392 }
393 
394 #else
395 ErrorCode TempestRemapper::convert_tempest_mesh_private( TempestMeshType meshType,
396  Mesh* mesh,
397  EntityHandle& mesh_set,
398  Range& entities,
399  Range* vertices )
400 {
401  ErrorCode rval;
402 
403  const bool outputEnabled = ( TempestRemapper::verbose && is_root );
404  const NodeVector& nodes = mesh->nodes;
405  const FaceVector& faces = mesh->faces;
406 
407  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
408  dbgprint.set_prefix( "[TempestToMOAB]: " );
409 
410  ReadUtilIface* iface;
411  MB_CHK_SET_ERR( m_interface->query_interface( iface ), "Can't get reader interface" );
412 
413  Tag gidTag = m_interface->globalId_tag();
414 
415  // Set the data for the vertices
416  std::vector< double* > arrays;
417  std::vector< int > gidsv( nodes.size() );
418  EntityHandle startv;
419  MB_CHK_SET_ERR( iface->get_node_coords( 3, nodes.size(), 0, startv, arrays ), "Can't get node coords" );
420  for( unsigned iverts = 0; iverts < nodes.size(); ++iverts )
421  {
422  const Node& node = nodes[iverts];
423  arrays[0][iverts] = node.x;
424  arrays[1][iverts] = node.y;
425  arrays[2][iverts] = node.z;
426  gidsv[iverts] = iverts + 1;
427  }
428  Range mbverts( startv, startv + nodes.size() - 1 );
429  MB_CHK_SET_ERR( m_interface->add_entities( mesh_set, mbverts ), "Can't add entities" );
430  MB_CHK_SET_ERR( m_interface->tag_set_data( gidTag, mbverts, &gidsv[0] ), "Can't set global_id tag" );
431 
432  gidsv.clear();
433  entities.clear();
434 
435  Tag srcParentTag, tgtParentTag;
436  std::vector< int > srcParent, tgtParent;
437  bool storeParentInfo = ( mesh->vecSourceFaceIx.size() > 0 );
438 
439  if( storeParentInfo )
440  {
441  int defaultInt = -1;
442  MB_CHK_SET_ERR( m_interface->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag,
443  MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt ),
444  "can't create positive tag" );
445 
446  MB_CHK_SET_ERR( m_interface->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag,
447  MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt ),
448  "can't create negative tag" );
449  }
450 
451  // Let us first perform a full pass assuming arbitrary polygons. This is especially true for
452  // overlap meshes.
453  // 1. We do a first pass over faces, decipher edge size and group into categories based on
454  // element type
455  // 2. Next we loop over type, and add blocks of elements into MOAB
456  // 3. For each block within the loop, also update the connectivity of elements.
457  {
458  if( outputEnabled )
459  dbgprint.printf( 0, "..Mesh size: Nodes [%zu] Elements [%zu].\n", nodes.size(), faces.size() );
460  const int NMAXPOLYEDGES = 15;
461  std::vector< unsigned > nPolys( NMAXPOLYEDGES, 0 );
462  std::vector< std::vector< int > > typeNSeqs( NMAXPOLYEDGES );
463  for( unsigned ifaces = 0; ifaces < faces.size(); ++ifaces )
464  {
465  const int iType = faces[ifaces].edges.size();
466  nPolys[iType]++;
467  typeNSeqs[iType].push_back( ifaces );
468  }
469  int iBlock = 0;
470  for( unsigned iType = 0; iType < NMAXPOLYEDGES; ++iType )
471  {
472  if( !nPolys[iType] ) continue; // Nothing to do
473 
474  const unsigned num_v_per_elem = iType;
475  EntityHandle starte; // Connectivity
476  EntityHandle* conn;
477 
478  // Allocate the connectivity array, depending on the element type
479  switch( num_v_per_elem )
480  {
481  case 3:
482  if( outputEnabled )
483  dbgprint.printf( 0, "....Block %d: Triangular Elements [%u].\n", iBlock++, nPolys[iType] );
484  MB_CHK_SET_ERR( iface->get_element_connect( nPolys[iType], num_v_per_elem, MBTRI, 0, starte, conn ),
485  "Can't get element connectivity" );
486  break;
487  case 4:
488  if( outputEnabled )
489  dbgprint.printf( 0, "....Block %d: Quadrilateral Elements [%u].\n", iBlock++, nPolys[iType] );
490  MB_CHK_SET_ERR( iface->get_element_connect( nPolys[iType], num_v_per_elem, MBQUAD, 0, starte,
491  conn ),
492  "Can't get element connectivity" );
493  break;
494  default:
495  if( outputEnabled )
496  dbgprint.printf( 0, "....Block %d: Polygonal [%u] Elements [%u].\n", iBlock++, iType,
497  nPolys[iType] );
498  MB_CHK_SET_ERR( iface->get_element_connect( nPolys[iType], num_v_per_elem, MBPOLYGON, 0, starte,
499  conn ),
500  "Can't get element connectivity" );
501  break;
502  }
503 
504  Range mbcells( starte, starte + nPolys[iType] - 1 );
505  m_interface->add_entities( mesh_set, mbcells );
506 
507  if( storeParentInfo )
508  {
509  srcParent.resize( mbcells.size(), -1 );
510  tgtParent.resize( mbcells.size(), -1 );
511  }
512 
513  std::vector< int > gids( typeNSeqs[iType].size() );
514  for( unsigned ifaces = 0, offset = 0; ifaces < typeNSeqs[iType].size(); ++ifaces )
515  {
516  const int fIndex = typeNSeqs[iType][ifaces];
517  const Face& face = faces[fIndex];
518  // conn[offset++] = startv + face.edges[0].node[0];
519  for( unsigned iedges = 0; iedges < face.edges.size(); ++iedges )
520  {
521  conn[offset++] = startv + face.edges[iedges].node[0];
522  }
523 
524  if( storeParentInfo )
525  {
526  srcParent[ifaces] = mesh->vecSourceFaceIx[fIndex] + 1;
527  tgtParent[ifaces] = mesh->vecTargetFaceIx[fIndex] + 1;
528  }
529 
530  gids[ifaces] = typeNSeqs[iType][ifaces] + 1;
531  }
532  MB_CHK_SET_ERR( m_interface->tag_set_data( gidTag, mbcells, &gids[0] ), "Can't set global_id tag" );
533 
534  if( meshType == OVERLAP_FILES )
535  {
536  // Now let us update the adjacency data, because some elements are new
537  MB_CHK_SET_ERR( iface->update_adjacencies( starte, nPolys[iType], num_v_per_elem, conn ),
538  "Can't update adjacencies" );
539  // Generate all adj entities dimension 1 and 2 (edges and faces/ tri or qua)
540  Range edges;
541  MB_CHK_SET_ERR( m_interface->get_adjacencies( mbcells, 1, true, edges, Interface::UNION ),
542  "Can't get edges" );
543  }
544 
545  if( storeParentInfo )
546  {
547  MB_CHK_SET_ERR( m_interface->tag_set_data( srcParentTag, mbcells, &srcParent[0] ),
548  "Can't set tag data" );
549  MB_CHK_SET_ERR( m_interface->tag_set_data( tgtParentTag, mbcells, &tgtParent[0] ),
550  "Can't set tag data" );
551  }
552  entities.merge( mbcells );
553  }
554  }
555 
556  if( vertices ) *vertices = mbverts;
557 
558  return MB_SUCCESS;
559 }
560 #endif
561 
562 ///////////////////////////////////////////////////////////////////////////////////
563 
564 ErrorCode TempestRemapper::ConvertAllMeshesToTempest()
565 {
566  // now, let us ensure we have a valid source and target mesh objects
567  // if not, convert the source and target meshes to TempestRemap format
568  if( this->m_source == nullptr )
569  {
570  // Convert the covering mesh to TempestRemap format
571  MB_CHK_SET_ERR( this->ConvertMeshToTempest( moab::Remapper::SourceMesh ),
572  "Can't convert source mesh to TempestRemap format" );
573  }
574  if( this->m_covering_source == nullptr )
575  {
576  // Convert the covering mesh to TempestRemap format
577  MB_CHK_SET_ERR( this->ConvertMeshToTempest( moab::Remapper::CoveringMesh ),
578  "Can't convert source coverage mesh to TempestRemap format" );
579  }
580  if( this->m_target == nullptr )
581  {
582  // Convert the covering mesh to TempestRemap format
583  MB_CHK_SET_ERR( this->ConvertMeshToTempest( moab::Remapper::TargetMesh ),
584  "Can't convert target mesh to TempestRemap format" );
585  }
586 
587  return moab::MB_SUCCESS;
588 }
589 
590 ErrorCode TempestRemapper::ConvertMeshToTempest( Remapper::IntersectionContext ctx )
591 {
592  const bool outputEnabled = ( TempestRemapper::verbose && is_root );
593 
594  // create a debug output object and set a prefix
595  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
596  dbgprint.set_prefix( "[MOABToTempest]: " );
597 
598  if( ctx == Remapper::SourceMesh )
599  {
600  if( !m_source ) m_source = new Mesh();
601  if( outputEnabled ) dbgprint.printf( 0, "Converting (source) MOAB to TempestRemap Mesh representation ...\n" );
602  MB_CHK_SET_ERR( convert_mesh_to_tempest_private( m_source, m_source_set, m_source_entities,
603  &m_source_vertices ),
604  "Can't convert source mesh to Tempest" );
605  if( m_source_entities.size() == 0 && m_source_vertices.size() != 0 )
606  {
607  this->point_cloud_source = true;
608  }
609  }
610  else if( ctx == Remapper::CoveringMesh )
611  {
612  if( !m_covering_source ) m_covering_source = new Mesh();
613  if( outputEnabled )
614  dbgprint.printf( 0, "Converting (covering source) MOAB to TempestRemap Mesh representation ...\n" );
615  MB_CHK_SET_ERR( convert_mesh_to_tempest_private( m_covering_source, m_covering_source_set,
616  m_covering_source_entities, &m_covering_source_vertices ),
617  "Can't convert convering source mesh to TempestRemap format" );
618  }
619  else if( ctx == Remapper::TargetMesh )
620  {
621  if( !m_target ) m_target = new Mesh();
622  if( outputEnabled ) dbgprint.printf( 0, "Converting (target) MOAB to TempestRemap Mesh representation ...\n" );
623  MB_CHK_SET_ERR( convert_mesh_to_tempest_private( m_target, m_target_set, m_target_entities,
624  &m_target_vertices ),
625  "Can't convert target mesh to Tempest" );
626  if( m_target_entities.size() == 0 && m_target_vertices.size() != 0 ) this->point_cloud_target = true;
627  }
628  else if( ctx == Remapper::OverlapMesh ) // Overlap mesh
629  {
630  if( !m_overlap ) m_overlap = new Mesh();
631  if( outputEnabled ) dbgprint.printf( 0, "Converting (overlap) MOAB to TempestRemap Mesh representation ...\n" );
632  MB_CHK_SET_ERR( ConvertOverlapMeshSourceOrdered(), "Can't convert overlap mesh to Tempest" );
633  }
634  else
635  {
636  MB_CHK_SET_ERR( MB_FAILURE, "Invalid IntersectionContext context provided" );
637  }
638 
639  return moab::MB_SUCCESS;
640 }
641 
642 ErrorCode TempestRemapper::convert_mesh_to_tempest_private( Mesh* mesh,
643  EntityHandle mesh_set,
644  moab::Range& elems,
645  moab::Range* pverts )
646 {
647  NodeVector& nodes = mesh->nodes;
648  FaceVector& faces = mesh->faces;
649 
650  elems.clear();
651  MB_CHK_ERR( m_interface->get_entities_by_dimension( mesh_set, 2, elems ) );
652 
653  const size_t nelems = elems.size();
654 
655  // resize the number of elements in Tempest mesh
656  faces.resize( nelems );
657 
658  // let us now get the vertices from all the elements
659  Range verts;
660  MB_CHK_ERR( m_interface->get_connectivity( elems, verts ) );
661  if( verts.size() == 0 )
662  {
663  MB_CHK_ERR( m_interface->get_entities_by_dimension( mesh_set, 0, verts ) );
664  }
665  // assert(verts.size() > 0); // If not, this may be an invalid mesh ! possible for unbalanced loads
666 
667  std::map< EntityHandle, int > indxMap;
668  bool useRange = true;
669  if( verts.compactness() > 0.01 )
670  {
671  int j = 0;
672  for( Range::iterator it = verts.begin(); it != verts.end(); ++it )
673  indxMap[*it] = j++;
674  useRange = false;
675  }
676 
677  std::vector< int > globIds( nelems );
678  moab::Tag gid = m_interface->globalId_tag();
679  MB_CHK_ERR( m_interface->tag_get_data( gid, elems, &globIds[0] ) );
680  std::vector< size_t > sortedIdx;
681  if( offlineWorkflow )
682  {
683  sortedIdx.resize( nelems );
684  // initialize original index locations
685  std::iota( sortedIdx.begin(), sortedIdx.end(), 0 );
686  // sort indexes based on comparing values in v, using std::stable_sort instead of std::sort
687  // to avoid unnecessary index re-orderings when v contains elements of equal values
688  std::sort( sortedIdx.begin(), sortedIdx.end(),
689  [&globIds]( size_t i1, size_t i2 ) { return globIds[i1] < globIds[i2]; } );
690  }
691 
692  for( unsigned iface = 0; iface < nelems; ++iface )
693  {
694  Face& face = faces[iface];
695  EntityHandle ehandle = ( offlineWorkflow ? elems[sortedIdx[iface]] : elems[iface] );
696 
697  // get the connectivity for each edge
698  const EntityHandle* connectface;
699  int nnodesf;
700  MB_CHK_ERR( m_interface->get_connectivity( ehandle, connectface, nnodesf ) );
701  // account for padded polygons
702  while( connectface[nnodesf - 2] == connectface[nnodesf - 1] && nnodesf > 3 )
703  nnodesf--;
704 
705  face.edges.resize( nnodesf );
706  for( int iverts = 0; iverts < nnodesf; ++iverts )
707  {
708  int indx = ( useRange ? verts.index( connectface[iverts] ) : indxMap[connectface[iverts]] );
709  assert( indx >= 0 );
710  face.SetNode( iverts, indx );
711  }
712  }
713 
714  unsigned nnodes = verts.size();
715  nodes.resize( nnodes );
716 
717  // Set the data for the vertices
718  std::vector< double > coordx( nnodes ), coordy( nnodes ), coordz( nnodes );
719  MB_CHK_ERR( m_interface->get_coords( verts, &coordx[0], &coordy[0], &coordz[0] ) );
720  for( unsigned inode = 0; inode < nnodes; ++inode )
721  {
722  Node& node = nodes[inode];
723  node.x = coordx[inode];
724  node.y = coordy[inode];
725  node.z = coordz[inode];
726  }
727  coordx.clear();
728  coordy.clear();
729  coordz.clear();
730 
731  mesh->RemoveCoincidentNodes();
732  mesh->RemoveZeroEdges();
733 
734  // Generate reverse node array and edge map
735  if( constructEdgeMap ) mesh->ConstructEdgeMap( false );
736  // mesh->ConstructReverseNodeArray();
737 
738  // mesh->Validate();
739 
740  if( pverts )
741  {
742  pverts->clear();
743  *pverts = verts;
744  }
745  verts.clear();
746 
747  return MB_SUCCESS;
748 }
749 
750 ///////////////////////////////////////////////////////////////////////////////////
751 
752 bool IntPairComparator( const std::array< int, 3 >& a, const std::array< int, 3 >& b )
753 {
754  if( std::get< 1 >( a ) == std::get< 1 >( b ) )
755  return std::get< 2 >( a ) < std::get< 2 >( b );
756  else
757  return std::get< 1 >( a ) < std::get< 1 >( b );
758 }
759 
760 moab::ErrorCode moab::TempestRemapper::GetOverlapAugmentedEntities( moab::Range& sharedGhostEntities )
761 {
762  sharedGhostEntities.clear();
763 #ifdef MOAB_HAVE_MPI
764 
765  // Remove entities in the intersection mesh that are part of the ghosted overlap
766  if( is_parallel )
767  {
768  moab::Range allents;
769  MB_CHK_SET_ERR( m_interface->get_entities_by_dimension( m_overlap_set, 2, allents ),
770  "Getting entities dim 2 failed" );
771 
772  moab::Range sharedents;
773  moab::Tag ghostTag;
774  std::vector< int > ghFlags( allents.size() );
775  MB_CHK_ERR( m_interface->tag_get_handle( "ORIG_PROC", ghostTag ) );
776  MB_CHK_ERR( m_interface->tag_get_data( ghostTag, allents, &ghFlags[0] ) );
777  for( unsigned i = 0; i < allents.size(); ++i )
778  if( ghFlags[i] >= 0 ) // it means it is a ghost overlap element
779  sharedents.insert( allents[i] ); // this should not participate in smat!
780 
781  allents = subtract( allents, sharedents );
782 
783  // Get connectivity from all ghosted elements and filter out
784  // the vertices that are not owned
785  moab::Range ownedverts, sharedverts;
786  MB_CHK_SET_ERR( m_interface->get_connectivity( allents, ownedverts ), "Deleting entities dim 0 failed" );
787  MB_CHK_SET_ERR( m_interface->get_connectivity( sharedents, sharedverts ), "Deleting entities dim 0 failed" );
788  sharedverts = subtract( sharedverts, ownedverts );
789  // MB_CHK_SET_ERR( m_interface->remove_entities(m_overlap_set, sharedents), // "Deleting entities dim 2 failed" ); MB_CHK_SET_ERR( m_interface->remove_entities(m_overlap_set,
790  // sharedverts), "Deleting entities dim 0 failed" );
791 
792  sharedGhostEntities.merge( sharedents );
793  // sharedGhostEntities.merge(sharedverts);
794  }
795 #endif
796  return moab::MB_SUCCESS;
797 }
798 
799 ErrorCode TempestRemapper::ConvertOverlapMeshSourceOrdered()
800 {
801  m_overlap_entities.clear();
802  MB_CHK_ERR( m_interface->get_entities_by_dimension( m_overlap_set, 2, m_overlap_entities ) );
803  size_t n_overlap_entitites = m_overlap_entities.size();
804 
805  // Allocate for the overlap mesh
806  if( !m_overlap ) m_overlap = new Mesh();
807 
808  std::vector< std::array< int, 3 > > sorted_overlap_order( n_overlap_entitites,
809  std::array< int, 3 >( { -1, -1, -1 } ) );
810  {
811  Tag srcParentTag, tgtParentTag;
812  MB_CHK_ERR( m_interface->tag_get_handle( "SourceParent", srcParentTag ) );
813  MB_CHK_ERR( m_interface->tag_get_handle( "TargetParent", tgtParentTag ) );
814  // Overlap mesh: resize the source and target connection arrays
815  m_overlap->vecTargetFaceIx.resize( n_overlap_entitites );
816  m_overlap->vecSourceFaceIx.resize( n_overlap_entitites );
817 
818  // We need a global to local numbering
819  Tag gidtag = m_interface->globalId_tag();
820  std::vector< int > gids_src( m_covering_source_entities.size(), -1 ), gids_tgt( m_target_entities.size(), -1 );
821  MB_CHK_ERR( m_interface->tag_get_data( gidtag, m_covering_source_entities, gids_src.data() ) );
822  MB_CHK_ERR( m_interface->tag_get_data( gidtag, m_target_entities, gids_tgt.data() ) );
823 
824  // Global-id -> local-index lookup.
825  //
826  // This used to be a std::find() over the whole gid vector for every overlap
827  // element, making the loop below O(n_overlap * n_source). On a 288k-cell RLL
828  // source that single std::find accounted for 68% of total runtime (13029 of
829  // 19243 samples). Build the index once instead, so each lookup is O(1) and
830  // the loop becomes O(n_overlap + n_source).
831  //
832  // Insert with first-occurrence-wins so the result matches what std::find
833  // returned when a global id appears more than once.
834  auto build_index = []( const std::vector< int >& gids ) {
835  std::unordered_map< int, int > index;
836  index.reserve( gids.size() * 2 );
837  for( size_t i = 0; i < gids.size(); i++ )
838  index.emplace( gids[i], static_cast< int >( i ) );
839  return index;
840  };
841 
842  std::unordered_map< int, int > index_src, index_tgt;
843  if( !offlineWorkflow )
844  {
845  index_src = build_index( gids_src );
846  index_tgt = build_index( gids_tgt );
847  }
848 
849  auto find_lid = [this]( const std::unordered_map< int, int >& index, int gid ) -> int {
850  if( offlineWorkflow ) return gid - 1;
851  auto it = index.find( gid );
852  return ( it != index.end() ? it->second : -1 );
853  };
854 
855  std::vector< int > ghFlags;
856  if( is_parallel )
857  {
858  Tag ghostTag;
859  ghFlags.resize( n_overlap_entitites );
860  MB_CHK_SET_ERR( m_interface->tag_get_handle( "ORIG_PROC", ghostTag ), "Could not find ORIG_PROC tag" );
861  MB_CHK_ERR( m_interface->tag_get_data( ghostTag, m_overlap_entities, ghFlags.data() ) );
862  }
863 
864  // Overlap mesh: resize the source and target connection arrays
865  std::vector< int > rbids_src( n_overlap_entitites ), rbids_tgt( n_overlap_entitites );
866  MB_CHK_ERR( m_interface->tag_get_data( srcParentTag, m_overlap_entities, rbids_src.data() ) );
867  MB_CHK_ERR( m_interface->tag_get_data( tgtParentTag, m_overlap_entities, rbids_tgt.data() ) );
868 
869  for( size_t ix = 0; ix < n_overlap_entitites; ++ix )
870  {
871  std::get< 0 >( sorted_overlap_order[ix] ) = ix;
872  std::get< 1 >( sorted_overlap_order[ix] ) = find_lid( index_src, rbids_src[ix] );
873  assert( std::get< 1 >( sorted_overlap_order[ix] ) >= 0 );
874  if( is_parallel && ghFlags[ix] >= 0 )
875  {
876  // Ghost overlap element: its target cell lives on the rank that owns this element.
877  // Weight contributions for this element are computed on that owning rank (where the
878  // element has ORIG_PROC=-1 and a valid local target face index). Using -1 here
879  // ensures the element is skipped on this rank to avoid double-counting.
880  std::get< 2 >( sorted_overlap_order[ix] ) = -1;
881  }
882  else
883  std::get< 2 >( sorted_overlap_order[ix] ) = find_lid( index_tgt, rbids_tgt[ix] );
884  }
885  // now sort the overlap elements such that they are ordered by source parent first
886  // and then target parent next
887  std::sort( sorted_overlap_order.begin(), sorted_overlap_order.end(), IntPairComparator );
888 
889  for( unsigned ie = 0; ie < n_overlap_entitites; ++ie )
890  {
891  m_overlap->vecSourceFaceIx[ie] = std::get< 1 >( sorted_overlap_order[ie] );
892  m_overlap->vecTargetFaceIx[ie] = std::get< 2 >( sorted_overlap_order[ie] );
893  }
894  }
895 
896  FaceVector& faces = m_overlap->faces;
897  faces.resize( n_overlap_entitites );
898 
899  Range verts;
900  // let us now get the vertices from all the elements
901  MB_CHK_ERR( m_interface->get_connectivity( m_overlap_entities, verts ) );
902 
903  std::map< EntityHandle, int > indxMap;
904  {
905  int j = 0;
906  for( Range::iterator it = verts.begin(); it != verts.end(); ++it )
907  indxMap[*it] = j++;
908  }
909 
910  for( unsigned ifac = 0; ifac < m_overlap_entities.size(); ++ifac )
911  {
912  const unsigned iface = std::get< 0 >( sorted_overlap_order[ifac] );
913  Face& face = faces[ifac];
914  EntityHandle ehandle = m_overlap_entities[iface];
915 
916  // get the connectivity for each edge
917  const EntityHandle* connectface;
918  int nnodesf;
919  MB_CHK_ERR( m_interface->get_connectivity( ehandle, connectface, nnodesf ) );
920 
921  face.edges.resize( nnodesf );
922  for( int iverts = 0; iverts < nnodesf; ++iverts )
923  {
924  int indx = indxMap[connectface[iverts]];
925  assert( indx >= 0 );
926  face.SetNode( iverts, indx );
927  }
928  }
929  indxMap.clear();
930  sorted_overlap_order.clear();
931 
932  unsigned nnodes = verts.size();
933  NodeVector& nodes = m_overlap->nodes;
934  nodes.resize( nnodes );
935 
936  // Set the data for the vertices
937  std::vector< double > coordx( nnodes ), coordy( nnodes ), coordz( nnodes );
938  MB_CHK_ERR( m_interface->get_coords( verts, &coordx[0], &coordy[0], &coordz[0] ) );
939  for( unsigned inode = 0; inode < nnodes; ++inode )
940  {
941  Node& node = nodes[inode];
942  node.x = coordx[inode];
943  node.y = coordy[inode];
944  node.z = coordz[inode];
945  }
946  coordx.clear();
947  coordy.clear();
948  coordz.clear();
949  verts.clear();
950 
951  // ideally, we should do these in MOAB after intersction computation
952  // m_overlap->RemoveZeroEdges();
953  // m_overlap->RemoveCoincidentNodes( false );
954 
955  // Generate reverse node array and edge map
956  // if ( constructEdgeMap ) m_overlap->ConstructEdgeMap(false);
957  // m_overlap->ConstructReverseNodeArray();
958 
959  // For now, comment out validation; not needed
960  // m_overlap->Validate();
961 
962  // all done. return success
963  return MB_SUCCESS;
964 }
965 
966 ///////////////////////////////////////////////////////////////////////////////////
967 
968 moab::ErrorCode moab::TempestRemapper::WriteTempestIntersectionMesh( std::string strOutputFileName,
969  const bool fAllParallel,
970  const bool fInputConcave,
971  const bool fOutputConcave )
972 {
973 
974  // Let us alos write out the TempestRemap equivalent so that we can do some verification checks
975  if( fAllParallel )
976  {
977  if( is_root && size == 1 )
978  {
979  this->m_source->CalculateFaceAreas( fInputConcave );
980  this->m_target->CalculateFaceAreas( fOutputConcave );
981  }
982  else
983  {
984  // Perform reduction and write from root processor
985  this->m_source->CalculateFaceAreas( fInputConcave );
986  this->m_covering_source->CalculateFaceAreas( fInputConcave );
987  this->m_target->CalculateFaceAreas( fOutputConcave );
988  }
989  }
990  else
991  {
992  this->m_source->CalculateFaceAreas( fInputConcave );
993  this->m_target->CalculateFaceAreas( fOutputConcave );
994  }
995 
996  // The overlap mesh is written by TempestRemap's own NetCDF writer (Mesh::Write). This is
997  // forwarded to libTempestRemap and cannot go through the MBNcDispatch layer; guard it in
998  // one place and skip with a notice when NetCDF is unavailable.
999 #ifdef MOAB_HAVE_NETCDF
1000  this->m_overlap->Write( strOutputFileName.c_str(), NcFile::Netcdf4 );
1001 #else
1002  if( is_root ) std::cout << "NetCDF is not configured. Intersection mesh write will be skipped...\n";
1003 #endif
1004 
1005  return moab::MB_SUCCESS;
1006 }
1007 
1008 void TempestRemapper::SetMeshSet( Remapper::IntersectionContext ctx /* Remapper::CoveringMesh*/,
1009  moab::EntityHandle mset,
1010  moab::Range* entities )
1011 {
1012 
1013  if( ctx == Remapper::SourceMesh ) // should not be used
1014  {
1015  m_source_set = mset;
1016  if( entities ) m_source_entities = *entities;
1017  }
1018  else if( ctx == Remapper::TargetMesh )
1019  {
1020  m_target_set = mset;
1021  if( entities ) m_target_entities = *entities;
1022  }
1023  else if( ctx == Remapper::CoveringMesh )
1024  {
1025  m_covering_source_set = mset;
1026  if( entities ) m_covering_source_entities = *entities;
1027  }
1028  else
1029  {
1030  // nothing to do really..
1031  return;
1032  }
1033 }
1034 
1035 ///////////////////////////////////////////////////////////////////////////////////
1036 
1037 #ifndef MOAB_HAVE_MPI
1038 ErrorCode TempestRemapper::assign_vertex_element_IDs( Tag idtag,
1039  EntityHandle this_set,
1040  const int dimension,
1041  const int start_id )
1042 {
1043  assert( idtag );
1044 
1045  ErrorCode rval;
1046  Range entities;
1047  MB_CHK_SET_ERR( m_interface->get_entities_by_dimension( this_set, dimension, entities ), "Failed to get entities" );
1048 
1049  if( entities.size() == 0 ) return moab::MB_SUCCESS;
1050 
1051  int idoffset = start_id;
1052  std::vector< int > gid( entities.size() );
1053  for( unsigned i = 0; i < entities.size(); ++i )
1054  gid[i] = idoffset++;
1055 
1056  MB_CHK_ERR( m_interface->tag_set_data( idtag, entities, &gid[0] ) );
1057 
1058  return moab::MB_SUCCESS;
1059 }
1060 #endif
1061 
1062 ///////////////////////////////////////////////////////////////////////////////
1063 
1064 // Create a custom comparator for Nodes
1065 bool operator<( Node const& lhs, Node const& rhs )
1066 {
1067  return std::pow( lhs.x - rhs.x, 2.0 ) + std::pow( lhs.y - rhs.y, 2.0 ) + std::pow( lhs.z - rhs.z, 2.0 );
1068 }
1069 
1070 ErrorCode TempestRemapper::GenerateCSMeshMetadata( const int ntot_elements,
1071  moab::Range& ents,
1072  moab::Range* secondary_ents,
1073  const std::string& dofTagName,
1074  int nP )
1075 {
1076  const int csResolution = std::sqrt( ntot_elements / 6.0 );
1077  if( csResolution * csResolution * 6 != ntot_elements ) return MB_INVALID_SIZE;
1078 
1079  // Create a temporary Cubed-Sphere mesh
1080  // NOTE: This will not work for RRM grids. Need to run HOMME for that case anyway
1081  Mesh csMesh;
1082  if( GenerateCSMesh( csMesh, csResolution, "", "NetCDF4" ) )
1083  MB_CHK_SET_ERR( moab::MB_FAILURE, // unsuccessful call
1084  "Failed to generate CS mesh through TempestRemap" );
1085 
1086  // let us now generate the mesh metadata
1087  if( this->GenerateMeshMetadata( csMesh, ntot_elements, ents, secondary_ents, dofTagName, nP ) )
1088  MB_CHK_SET_ERR( moab::MB_FAILURE, "Failed in call to GenerateMeshMetadata" ); // unsuccessful call
1089 
1090  return moab::MB_SUCCESS;
1091 }
1092 
1093 ErrorCode TempestRemapper::GenerateMeshMetadata( Mesh& csMesh,
1094  const int ntot_elements,
1095  moab::Range& ents,
1096  moab::Range* secondary_ents,
1097  const std::string& dofTagName,
1098  int nP )
1099 {
1100  Tag dofTag;
1101  bool created = false;
1102  MB_CHK_SET_ERR( m_interface->tag_get_handle( dofTagName.c_str(), nP * nP, MB_TYPE_INTEGER, dofTag,
1103  MB_TAG_DENSE | MB_TAG_CREAT, 0, &created ),
1104  "Failed creating DoF tag" );
1105 
1106  // Number of Faces
1107  int nElements = static_cast< int >( csMesh.faces.size() );
1108 
1109  if( nElements != ntot_elements ) return MB_INVALID_SIZE;
1110 
1111  // Initialize data structures
1112  DataArray3D< int > dataGLLnodes;
1113  dataGLLnodes.Allocate( nP, nP, nElements );
1114 
1115  std::map< Node, int > mapNodes;
1116  std::map< Node, moab::EntityHandle > mapLocalMBNodes;
1117 
1118  // GLL Quadrature nodes
1119  DataArray1D< double > dG;
1120  DataArray1D< double > dW;
1121  GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );
1122 
1123  moab::Range entities( ents );
1124  if( secondary_ents ) entities.insert( secondary_ents->begin(), secondary_ents->end() );
1125  double elcoords[3];
1126  for( unsigned iel = 0; iel < entities.size(); ++iel )
1127  {
1128  EntityHandle eh = entities[iel];
1129  MB_CHK_SET_ERR( m_interface->get_coords( &eh, 1, elcoords ), "failed to get element coordinates" );
1130  Node elCentroid( elcoords[0], elcoords[1], elcoords[2] );
1131  mapLocalMBNodes.insert( std::pair< Node, moab::EntityHandle >( elCentroid, eh ) );
1132  }
1133 
1134  // Build a Kd-tree for local mesh (nearest neighbor searches)
1135  // Loop over all elements in CS-Mesh
1136  // Then find if current centroid is in an element
1137  // If yes - then let us compute the DoF numbering and set to tag data
1138  // If no - then compute DoF numbering BUT DO NOT SET to tag data
1139  // continue
1140  int* dofIDs = new int[nP * nP];
1141 
1142  // Write metadata
1143  for( int k = 0; k < nElements; k++ )
1144  {
1145  const Face& face = csMesh.faces[k];
1146  const NodeVector& nodes = csMesh.nodes;
1147 
1148  if( face.edges.size() != 4 )
1149  {
1150  _EXCEPTIONT( "Mesh must only contain quadrilateral elements" );
1151  }
1152 
1153  Node centroid;
1154  centroid.x = centroid.y = centroid.z = 0.0;
1155  for( unsigned l = 0; l < face.edges.size(); ++l )
1156  {
1157  centroid.x += nodes[face[l]].x;
1158  centroid.y += nodes[face[l]].y;
1159  centroid.z += nodes[face[l]].z;
1160  }
1161  const double factor = 1.0 / face.edges.size();
1162  centroid.x *= factor;
1163  centroid.y *= factor;
1164  centroid.z *= factor;
1165 
1166  bool locElem = false;
1167  EntityHandle current_eh;
1168  if( mapLocalMBNodes.find( centroid ) != mapLocalMBNodes.end() )
1169  {
1170  locElem = true;
1171  current_eh = mapLocalMBNodes[centroid];
1172  }
1173 
1174  for( int j = 0; j < nP; j++ )
1175  {
1176  for( int i = 0; i < nP; i++ )
1177  {
1178  // Get local map vectors
1179  Node nodeGLL;
1180  Node dDx1G;
1181  Node dDx2G;
1182 
1183  // ApplyLocalMap(
1184  // face,
1185  // nodevec,
1186  // dG[i],
1187  // dG[j],
1188  // nodeGLL,
1189  // dDx1G,
1190  // dDx2G);
1191  const double& dAlpha = dG[i];
1192  const double& dBeta = dG[j];
1193 
1194  // Calculate nodal locations on the plane
1195  double dXc = nodes[face[0]].x * ( 1.0 - dAlpha ) * ( 1.0 - dBeta ) +
1196  nodes[face[1]].x * dAlpha * ( 1.0 - dBeta ) + nodes[face[2]].x * dAlpha * dBeta +
1197  nodes[face[3]].x * ( 1.0 - dAlpha ) * dBeta;
1198 
1199  double dYc = nodes[face[0]].y * ( 1.0 - dAlpha ) * ( 1.0 - dBeta ) +
1200  nodes[face[1]].y * dAlpha * ( 1.0 - dBeta ) + nodes[face[2]].y * dAlpha * dBeta +
1201  nodes[face[3]].y * ( 1.0 - dAlpha ) * dBeta;
1202 
1203  double dZc = nodes[face[0]].z * ( 1.0 - dAlpha ) * ( 1.0 - dBeta ) +
1204  nodes[face[1]].z * dAlpha * ( 1.0 - dBeta ) + nodes[face[2]].z * dAlpha * dBeta +
1205  nodes[face[3]].z * ( 1.0 - dAlpha ) * dBeta;
1206 
1207  double dR = sqrt( dXc * dXc + dYc * dYc + dZc * dZc );
1208 
1209  // Mapped node location
1210  nodeGLL.x = dXc / dR;
1211  nodeGLL.y = dYc / dR;
1212  nodeGLL.z = dZc / dR;
1213 
1214  // Determine if this is a unique Node
1215  std::map< Node, int >::const_iterator iter = mapNodes.find( nodeGLL );
1216  if( iter == mapNodes.end() )
1217  {
1218  // Insert new unique node into map
1219  int ixNode = static_cast< int >( mapNodes.size() );
1220  mapNodes.insert( std::pair< Node, int >( nodeGLL, ixNode ) );
1221  dataGLLnodes[j][i][k] = ixNode + 1;
1222  }
1223  else
1224  {
1225  dataGLLnodes[j][i][k] = iter->second + 1;
1226  }
1227 
1228  dofIDs[j * nP + i] = dataGLLnodes[j][i][k];
1229  }
1230  }
1231 
1232  if( locElem )
1233  {
1234  MB_CHK_SET_ERR( m_interface->tag_set_data( dofTag, &current_eh, 1, dofIDs ),
1235  "Failed to tag_set_data for DoFs" );
1236  }
1237  }
1238 
1239  // clear memory
1240  delete[] dofIDs;
1241  mapLocalMBNodes.clear();
1242  mapNodes.clear();
1243 
1244  return moab::MB_SUCCESS;
1245 }
1246 
1247 ///////////////////////////////////////////////////////////////////////////////////
1248 
1249 //#define MOAB_DBG
1250 ErrorCode TempestRemapper::validate_global_ids_private( EntityHandle mesh_set,
1251  int dimension,
1252  const char* mesh_name,
1253  std::string& error_message )
1254 {
1255  Range ents;
1256  MB_CHK_ERR( m_interface->get_entities_by_dimension( mesh_set, dimension, ents ) );
1257  if( ents.empty() ) return MB_SUCCESS; // nothing to check (e.g. empty partition)
1258 
1259  const char* what = ( 0 == dimension ? "vertex" : "element" );
1260 
1261  std::vector< int > gids( ents.size(), -1 );
1262  MB_CHK_ERR( m_interface->tag_get_data( m_interface->globalId_tag(), ents, gids.data() ) );
1263 
1264  // 1) every id must be strictly positive. A missing tag reads back as the
1265  // dense default of -1, which is the common failure mode.
1266  size_t n_bad = 0;
1267  int first_bad = 0;
1268  for( size_t i = 0; i < gids.size(); ++i )
1269  if( gids[i] <= 0 )
1270  {
1271  if( !n_bad ) first_bad = gids[i];
1272  ++n_bad;
1273  }
1274 
1275  if( n_bad )
1276  {
1277  std::ostringstream os;
1278  os << mesh_name << " mesh has " << n_bad << " of " << ents.size() << " " << what
1279  << " entities with a non-positive GLOBAL_ID (first bad value = " << first_bad << ")";
1280  if( n_bad == ents.size() )
1281  os << "; the GLOBAL_ID tag is most likely absent, so every entry is the tag default";
1282  error_message = os.str();
1283  return MB_FAILURE;
1284  }
1285 
1286  // 2) ids must be locally unique. Duplicates make entity identification
1287  // ambiguous during migration and corrupt the map rows/columns.
1288  std::vector< int > sorted( gids );
1289  std::sort( sorted.begin(), sorted.end() );
1290  std::vector< int >::iterator dup = std::adjacent_find( sorted.begin(), sorted.end() );
1291  if( dup != sorted.end() )
1292  {
1293  const size_t n_unique = std::distance( sorted.begin(), std::unique( sorted.begin(), sorted.end() ) );
1294  std::ostringstream os;
1295  os << mesh_name << " mesh has duplicate " << what << " GLOBAL_IDs: " << ents.size() << " entities but only "
1296  << n_unique << " distinct ids (e.g. id " << *dup << " appears more than once)";
1297  error_message = os.str();
1298  return MB_FAILURE;
1299  }
1300 
1301  return MB_SUCCESS;
1302 }
1303 
1304 ErrorCode TempestRemapper::ValidateGlobalIds( bool throw_error )
1305 {
1306  // Element ids are always required: they become the rows and columns of the map.
1307  //
1308  // Vertex ids are only required in parallel, where construct_covering_set() keys on
1309  // them to identify the corners of migrated source cells -- unnumbered vertices
1310  // collapse every corner to one handle and produce degenerate cells. Nothing
1311  // migrates in serial, so the vertex ids are never consulted and demanding them
1312  // rejects meshes that work perfectly well: in-memory NetCDF domain meshes, for
1313  // instance, carry element ids but no vertex ids.
1314  const bool check_vertices = is_parallel;
1315 
1316  const EntityHandle sets[2] = { m_source_set, m_target_set };
1317  const char* names[2] = { "source", "target" };
1318 
1319  std::string message;
1320  int local_bad = 0;
1321  for( int im = 0; im < 2 && !local_bad; ++im )
1322  {
1323  if( !sets[im] ) continue;
1324  for( int dim = ( check_vertices ? 0 : 2 ); dim <= 2; dim += 2 ) // vertices (0) and faces (2)
1325  {
1326  if( MB_SUCCESS != validate_global_ids_private( sets[im], dim, names[im], message ) )
1327  {
1328  local_bad = 1;
1329  break;
1330  }
1331  }
1332  }
1333 
1334  int global_bad = local_bad;
1335 #ifdef MOAB_HAVE_MPI
1336  if( m_pcomm )
1337  {
1338  // Collective: a mesh is only valid if it is valid on every rank, and every
1339  // rank must agree so that they fail (or continue) together.
1340  MPI_Allreduce( &local_bad, &global_bad, 1, MPI_INT, MPI_MAX, m_pcomm->comm() );
1341  }
1342 #endif
1343 
1344  if( global_bad )
1345  {
1346  if( local_bad )
1347  std::cout << "[ERROR] rank " << rank << ": " << message << std::endl;
1348  else if( is_root )
1349  std::cout << "[ERROR] invalid GLOBAL_IDs detected on another rank" << std::endl;
1350 
1351  if( throw_error )
1352  MB_SET_ERR( MB_FAILURE, "Invalid GLOBAL_ID tags on input meshes; refusing to proceed. "
1353  "Ensure both meshes carry unique, strictly positive GLOBAL_IDs on "
1354  "vertices and elements (e.g. re-partition with 'mbpart -j', or read "
1355  "SCRIP/NetCDF sources in parallel so ids are assigned)." );
1356  return MB_FAILURE;
1357  }
1358 
1359  return MB_SUCCESS;
1360 }
1361 
1362 ErrorCode TempestRemapper::ConstructCoveringSet( double tolerance,
1363  double radius_src,
1364  double radius_tgt,
1365  double boxeps,
1366  bool regional_mesh,
1367  bool gnomonic,
1368  int nb_ghost_layers )
1369 {
1370  if( nb_ghost_layers >= 1 ) gnomonic = false;
1371  rrmgrids = regional_mesh;
1372  moab::Range local_verts;
1373 
1374  // Fail fast on meshes whose GLOBAL_IDs cannot support the coverage migration or
1375  // the map assembly. Doing this before any expensive work turns what used to be a
1376  // silent corruption (collapsed cells, bogus map rows) into an immediate diagnosis.
1377  MB_CHK_ERR( ValidateGlobalIds() );
1378 
1379  // Initialize intersection context
1380  // Use the remapper-wide area formula rather than a hardcoded one, so that a single
1381  // choice (see SetAreaMethod) governs intersection, coverage and orientation fixups.
1382  mbintx = new moab::Intx2MeshOnSphere( m_interface, m_area_method );
1383 
1384  mbintx->set_error_tolerance( tolerance );
1385  mbintx->set_radius_source_mesh( radius_src );
1386  mbintx->set_radius_destination_mesh( radius_tgt );
1387  mbintx->set_box_error( boxeps );
1388 #ifdef MOAB_HAVE_MPI
1389  mbintx->set_parallel_comm( m_pcomm );
1390 #endif
1391 
1392  // compute the maxiumum edges in elements comprising source and target mesh
1393  MB_CHK_ERR( mbintx->FindMaxEdges( m_source_set, m_target_set ) );
1394 
1395  this->max_source_edges = mbintx->max_edges_1;
1396  this->max_target_edges = mbintx->max_edges_2;
1397 
1398  // Note: lots of communication possible, if mesh is distributed very differently
1399 #ifdef MOAB_HAVE_MPI
1400  if( is_parallel )
1401  {
1402  MB_CHK_ERR( mbintx->build_processor_euler_boxes( m_target_set, local_verts, gnomonic ) );
1403 
1404  MB_CHK_SET_ERR( m_interface->create_meshset( moab::MESHSET_SET, m_covering_source_set ),
1405  "Can't create new set" );
1406 
1407  MB_CHK_ERR( mbintx->construct_covering_set( m_source_set, m_covering_source_set, gnomonic, nb_ghost_layers ) );
1408 #ifdef MOAB_DBG
1409  std::stringstream filename;
1410  filename << "covering" << rank << ".h5m";
1411  MB_CHK_ERR( m_interface->write_file( filename.str().c_str(), 0, 0, &m_covering_source_set, 1 ) );
1412  std::stringstream targetFile;
1413  targetFile << "target" << rank << ".h5m";
1414  MB_CHK_ERR( m_interface->write_file( targetFile.str().c_str(), 0, 0, &m_target_set, 1 ) );
1415 #endif
1416  }
1417  else
1418  {
1419 #endif
1420  // Serial: every source cell is trivially in the coverage set.
1421  //
1422  // A "regional mesh" shortcut used to live here, picking for each target vertex
1423  // the source cell with the nearest centroid. Nearest centroid is not the same
1424  // relation as overlaps, so it dropped most of the covering set and produced a
1425  // silently wrong map (row sums down to 0.004). Holes are handled where they
1426  // belong instead -- in the overlap-based area correction; see IsRegionalMesh().
1427  m_covering_source_set = m_source_set;
1428  m_covering_source = m_source;
1429  m_covering_source_entities = m_source_entities; // this is a tempest mesh object; careful about
1430  // incrementing the reference?
1431  m_covering_source_vertices = m_source_vertices; // this is a tempest mesh object; careful about
1432  // incrementing the reference?
1433 #ifdef MOAB_HAVE_MPI
1434  }
1435 #endif
1436 
1437  // Convert the source, target and coverage meshes to TempestRemap format
1438  MB_CHK_ERR( ConvertAllMeshesToTempest() );
1439 
1440  return moab::MB_SUCCESS;
1441 }
1442 #undef MOAB_DBG
1443 //#define MOAB_DBG
1444 
1445 ErrorCode TempestRemapper::ComputeOverlapMesh( bool kdtree_search, bool use_tempest )
1446 {
1447  const bool outputEnabled = ( this->rank == 0 );
1448  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
1449  dbgprint.set_prefix( "[ComputeOverlapMesh]: " );
1450 
1451  //
1452  // Create the intersection on the sphere object and set up necessary parameters
1453  //
1454  // First, split based on whether to use TempestRemap or MOAB for intersection
1455  // If TempestRemap,
1456  // 1) Check for valid Mesh and pointers to objects for source/target
1457  // 2) Invoke GenerateOverlapWithMeshes routine from Tempest library
1458  // If MOAB,
1459  // 1) Check for valid source and target meshsets (and entities)
1460  // 2) Build processor bounding boxes and construct a covering set
1461  // 3) Perform intersection between the source (covering) and target entities
1462  if( use_tempest )
1463  {
1464  // Now let us construct the overlap mesh, by calling TempestRemap interface directly
1465  // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
1466  assert( m_covering_source != nullptr );
1467  assert( m_target != nullptr );
1468  if( m_overlap != nullptr ) delete m_overlap;
1469  bool concaveMeshA = false, concaveMeshB = false;
1470  // we have reset the overlap mesh - allocate now
1471  m_overlap = new Mesh();
1472  // Generate the overlap mesh using TempestRemap
1473  if( GenerateOverlapWithMeshes( *m_covering_source, *m_target, *m_overlap, "" /*outFilename*/, "Netcdf4",
1474  "exact", concaveMeshA, concaveMeshB, true, false ) )
1475  MB_CHK_SET_ERR( MB_FAILURE, "TempestRemap: cannot compute the intersection of meshes on the sphere" );
1476  }
1477  else
1478  {
1479  // Now perform the actual parallel intersection between the source and the target meshes
1480  if( kdtree_search )
1481  {
1482  if( outputEnabled ) dbgprint.printf( 0, "Computing intersection mesh with the Kd-tree search algorithm" );
1483  MB_CHK_SET_ERR( mbintx->intersect_meshes_kdtree( m_covering_source_set, m_target_set, m_overlap_set ),
1484  "Can't compute the intersection of meshes on the sphere with kd-tree" );
1485  }
1486  else
1487  {
1488  if( outputEnabled )
1489  dbgprint.printf( 0, "Computing intersection mesh with the advancing-front propagation algorithm" );
1490  MB_CHK_SET_ERR( mbintx->intersect_meshes( m_covering_source_set, m_target_set, m_overlap_set ),
1491  "Can't compute the intersection of meshes on the sphere" );
1492  }
1493 
1494 #ifdef MOAB_HAVE_MPI
1495  if( is_parallel )
1496  {
1497 #ifdef VERBOSE
1498  std::stringstream ffc, fft, ffo;
1499  ffc << "cover_" << rank << ".h5m";
1500  fft << "target_" << rank << ".h5m";
1501  ffo << "intx_" << rank << ".h5m";
1502  MB_CHK_ERR( m_interface->write_mesh( ffc.str().c_str(), &m_covering_source_set, 1 ) );
1503  MB_CHK_ERR( m_interface->write_mesh( fft.str().c_str(), &m_target_set, 1 ) );
1504  MB_CHK_ERR( m_interface->write_mesh( ffo.str().c_str(), &m_overlap_set, 1 ) );
1505 #endif
1506  // because we do not want to work with elements in coverage set that do not participate
1507  // in intersection, remove them from the coverage set we will not delete them yet, just
1508  // remove from the set !
1509  if( !point_cloud_target )
1510  {
1511  Range covEnts;
1512  MB_CHK_ERR( m_interface->get_entities_by_dimension( m_covering_source_set, 2, covEnts ) );
1513 
1514  std::map< int, int > loc_gid_to_lid_covsrc;
1515  std::vector< int > gids( covEnts.size(), -1 );
1516 
1517  Tag gidtag = m_interface->globalId_tag();
1518  MB_CHK_ERR( m_interface->tag_get_data( gidtag, covEnts, gids.data() ) );
1519 
1520  for( unsigned ie = 0; ie < gids.size(); ++ie )
1521  {
1522  assert( gids[ie] > 0 );
1523  loc_gid_to_lid_covsrc[gids[ie]] = ie;
1524  }
1525 
1526  Range intxCov, intxCells;
1527  Tag srcParentTag;
1528  MB_CHK_ERR( m_interface->tag_get_handle( "SourceParent", srcParentTag ) );
1529  MB_CHK_ERR( m_interface->get_entities_by_dimension( m_overlap_set, 2, intxCells ) );
1530  for( Range::iterator it = intxCells.begin(); it != intxCells.end(); ++it )
1531  {
1532  EntityHandle intxCell = *it;
1533  int srcParent = -1;
1534  MB_CHK_ERR( m_interface->tag_get_data( srcParentTag, &intxCell, 1, &srcParent ) );
1535 
1536  assert( srcParent >= 0 );
1537  intxCov.insert( covEnts[loc_gid_to_lid_covsrc[srcParent]] );
1538  }
1539 
1540  Range notNeededCovCells = moab::subtract( covEnts, intxCov );
1541 
1542  // now let us get only the covering entities that participate in intersection mesh
1543  covEnts = moab::subtract( covEnts, notNeededCovCells );
1544 
1545  // in order for getting 1-ring neighborhood, we need to be sure that the adjacencies are updated (created)
1546  if( false )
1547  {
1548  // update all adjacency list
1549  Core* mb = dynamic_cast< Core* >( m_interface );
1550  AEntityFactory* adj_fact = mb->a_entity_factory();
1551  if( !adj_fact->vert_elem_adjacencies() )
1552  adj_fact->create_vert_elem_adjacencies();
1553  else
1554  {
1555  for( Range::iterator it = covEnts.begin(); it != covEnts.end(); ++it )
1556  {
1557  EntityHandle eh = *it;
1558  const EntityHandle* conn = nullptr;
1559  int num_nodes = 0;
1560  MB_CHK_ERR( mb->get_connectivity( eh, conn, num_nodes ) );
1561  adj_fact->notify_create_entity( eh, conn, num_nodes );
1562  }
1563  }
1564 
1565  // next, for elements on the edge of the partition, get one ring adjacencies
1566  Skinner skinner( mb );
1567  Range skin;
1568  MB_CHK_SET_ERR( skinner.find_skin( m_covering_source_set, covEnts, false, skin ),
1569  "Unable to find skin" );
1570  for( Range::iterator it = skin.begin(); it != skin.end(); ++it )
1571  {
1572  const EntityHandle* conn = nullptr;
1573  int len = 0;
1574  MB_CHK_ERR( mb->get_connectivity( *it, conn, len, false ) );
1575  for( int ie = 0; ie < len; ++ie )
1576  {
1577  std::vector< EntityHandle > adjacent_entities;
1578  MB_CHK_ERR( adj_fact->get_adjacencies( conn[ie], 2, false, adjacent_entities ) );
1579  for( auto ent : adjacent_entities )
1580  notNeededCovCells.erase( ent ); // ent is part of the 1-ring neighborhood
1581  }
1582  }
1583 
1584  MB_CHK_ERR( m_interface->write_mesh(
1585  std::string( "sourcecoveragemesh_p" + std::to_string( rank ) + ".h5m" ).c_str(),
1586  &m_covering_source_set, 1 ) );
1587  }
1588 
1589  // remove now from coverage set the cells that are not needed
1590  // MB_CHK_ERR( m_interface->remove_entities( m_covering_source_set, notNeededCovCells ) );
1591 
1592  // Need to loop over covEnts now and ensure at least N-rings are available dependign on whether bilinear (1) or
1593  // high order FV (p) methods are being used for map generation. For bilinear/FV(1): need 1 ring, and for FV(p)
1594  // need p=ring neighborhood to recover exact conservation and consistency wrt serial/parallel.
1595 #ifdef VERBOSE
1596  std::cout << " total participating elements in the covering set: " << intxCov.size() << "\n";
1597  std::cout << " remove from coverage set elements that are not intersected: " << notNeededCovCells.size()
1598  << "\n";
1599 #endif
1600  if( size > 1 )
1601  {
1602  // some source elements cover multiple target partitions; the conservation logic
1603  // requires to know all overlap elements for a source element; they need to be
1604  // communicated from the other target partitions
1605  //
1606  // so first we have to identify source (coverage) elements that cover multiple
1607  // target partitions
1608  //
1609  // we will then mark the source, we will need to migrate the overlap elements
1610  // that cover this to the original source for the source element; then
1611  // distribute the overlap elements to all processors that have the coverage mesh
1612  // used
1613  MB_CHK_ERR( AugmentOverlapSet() );
1614  }
1615  }
1616  }
1617 #endif
1618 
1619  // Fix any inconsistencies in the overlap mesh
1620  {
1621  IntxAreaUtils areaAdaptor( m_area_method );
1622  MB_CHK_ERR( IntxUtils::fix_degenerate_quads( m_interface, m_overlap_set ) );
1623  MB_CHK_ERR( areaAdaptor.positive_orientation( m_interface, m_overlap_set, 1.0 /*radius*/ ) );
1624  }
1625 
1626  // free the memory
1627  delete mbintx;
1628  }
1629 
1630  // Now, let us convert the overlap mesh to MOAB format so that we have a consistent interface
1631  MB_CHK_SET_ERR( ConvertOverlapMeshSourceOrdered(), "Can't convert overlap TempestRemap mesh to MOAB format" );
1632 
1633  return MB_SUCCESS;
1634 }
1635 
1636 #ifdef MOAB_HAVE_MPI
1637 
1638 // this function is called only in parallel
1639 ///////////////////////////////////////////////////////////////////////////////////
1640 ErrorCode TempestRemapper::AugmentOverlapSet()
1641 {
1642  /*
1643  * overall strategy:
1644  *
1645  * 1) collect all boundary target cells on the current task, affected by the partition boundary;
1646  * note: not only partition boundary, we need all boundary (all coastal lines) and partition
1647  * boundary targetBoundaryIds is the set of target boundary cell IDs
1648  *
1649  * 2) collect all source cells that are intersecting boundary cells (call them
1650  * affectedSourceCellsIds)
1651  *
1652  * 3) collect overlap, that is accumulate all overlap cells that have source target in
1653  * affectedSourceCellsIds
1654  */
1655  // first, get all edges on the partition boundary, on the target mesh, then all the target
1656  // elements that border the partition boundary
1657  Skinner skinner( m_interface );
1658 
1659  // now let us find all the boundary edges
1660  Range boundaryEdges;
1661  MB_CHK_ERR( skinner.find_skin( 0, this->m_target_entities, false, boundaryEdges ) );
1662 
1663  // filter boundary edges that are on partition boundary, not on boundary
1664  // find all cells adjacent to these boundary edges, from target set
1665  Range boundaryCells; // these will be filtered from target_set
1666  MB_CHK_ERR( m_interface->get_adjacencies( boundaryEdges, 2, false, boundaryCells, Interface::UNION ) );
1667  boundaryCells = intersect( boundaryCells, this->m_target_entities );
1668 
1669 #ifdef MOAB_DBG
1670  EntityHandle tmpSet;
1671  MB_CHK_SET_ERR( m_interface->create_meshset( MESHSET_SET, tmpSet ), "Can't create temporary set" );
1672  // add the boundary set and edges, and save it to a file
1673  MB_CHK_SET_ERR( m_interface->add_entities( tmpSet, boundaryCells ), "Can't add entities" );
1674  MB_CHK_SET_ERR( m_interface->add_entities( tmpSet, boundaryEdges ), "Can't add edges" );
1675  std::stringstream ffs;
1676  ffs << "boundaryCells_0" << rank << ".h5m";
1677  MB_CHK_ERR( m_interface->write_mesh( ffs.str().c_str(), &tmpSet, 1 ) );
1678 #endif
1679 
1680  // now that we have the boundary cells, see which overlap polys have these as parents;
1681  // find the ids of the boundary cells;
1682  Tag gid = m_interface->globalId_tag();
1683  std::set< int > targetBoundaryIds;
1684  for( Range::iterator it = boundaryCells.begin(); it != boundaryCells.end(); ++it )
1685  {
1686  int tid;
1687  EntityHandle targetCell = *it;
1688  MB_CHK_SET_ERR( m_interface->tag_get_data( gid, &targetCell, 1, &tid ),
1689  "Can't get global id tag on target cell" );
1690  if( tid < 0 ) std::cout << " incorrect id for a target cell\n";
1691  targetBoundaryIds.insert( tid );
1692  }
1693 
1694  Range overlapCells;
1695  MB_CHK_ERR( m_interface->get_entities_by_dimension( m_overlap_set, 2, overlapCells ) );
1696 
1697  std::set< int > affectedSourceCellsIds;
1698  Tag targetParentTag, sourceParentTag; // do not use blue/red, as it is more confusing
1699  MB_CHK_ERR( m_interface->tag_get_handle( "TargetParent", targetParentTag ) );
1700  MB_CHK_ERR( m_interface->tag_get_handle( "SourceParent", sourceParentTag ) );
1701  for( Range::iterator it = overlapCells.begin(); it != overlapCells.end(); ++it )
1702  {
1703  EntityHandle intxCell = *it;
1704  int targetParentID, sourceParentID;
1705  MB_CHK_ERR( m_interface->tag_get_data( targetParentTag, &intxCell, 1, &targetParentID ) );
1706  if( targetBoundaryIds.find( targetParentID ) != targetBoundaryIds.end() )
1707  {
1708  // this means that the source element is affected
1709  MB_CHK_ERR( m_interface->tag_get_data( sourceParentTag, &intxCell, 1, &sourceParentID ) );
1710  affectedSourceCellsIds.insert( sourceParentID );
1711  }
1712  }
1713 
1714  // Now find all source cells affected, based on their id;
1715  // ( we do not yet have a global to local mapping for covering source mesh )
1716  // So create a map from source cell id to the eh;
1717  // it is needed to find out the original owning processor
1718  // this one came from, so to know where to send the overlap elements
1719  std::map< int, EntityHandle > affectedCovCellFromID;
1720 
1721  // use std::set<EntityHandle> instead of moab::Range for collecting cells,
1722  // either on coverage or target or intersection cells
1723  // their overlap cells will be sent to their original task, then distributed to all
1724  // other processes that might need them to compute conservation
1725  std::set< EntityHandle > affectedCovCells;
1726 
1727  Range covCells;
1728  MB_CHK_ERR( m_interface->get_entities_by_dimension( m_covering_source_set, 2, covCells ) );
1729  // loop thru all cov cells, to find the ones with global ids in affectedSourceCellsIds
1730  for( Range::iterator it = covCells.begin(); it != covCells.end(); ++it )
1731  {
1732  EntityHandle covCell = *it; //
1733  int covID;
1734  MB_CHK_ERR( m_interface->tag_get_data( gid, &covCell, 1, &covID ) );
1735  if( affectedSourceCellsIds.find( covID ) != affectedSourceCellsIds.end() )
1736  {
1737  // this source cell is affected;
1738  affectedCovCellFromID[covID] = covCell;
1739  affectedCovCells.insert( covCell );
1740  }
1741  }
1742 
1743  // now loop again over all overlap cells, to see if their source parent is "affected"
1744  // store in ranges the overlap cells that need to be sent to original task of the source cell
1745  // from there, they will be redistributed to the tasks that need that coverage cell
1746  Tag sendProcTag;
1747  MB_CHK_ERR( m_interface->tag_get_handle( "sending_processor", 1, MB_TYPE_INTEGER, sendProcTag ) );
1748 
1749  // basically a map from original processor task to the set of overlap cells to be sent there
1750  std::map< int, std::set< EntityHandle > > overlapCellsForTask;
1751  // this set will contain all intx cells that will need to be sent ( a union of above sets ,
1752  // that are organized per task on the above map )
1753  std::set< EntityHandle > overlapCellsToSend;
1754 
1755  for( Range::iterator it = overlapCells.begin(); it != overlapCells.end(); ++it )
1756  {
1757  EntityHandle intxCell = *it;
1758  int sourceParentID;
1759  MB_CHK_ERR( m_interface->tag_get_data( sourceParentTag, &intxCell, 1, &sourceParentID ) );
1760 
1761  if( affectedSourceCellsIds.find( sourceParentID ) != affectedSourceCellsIds.end() )
1762  {
1763  int orgTask = -1; // the original task that this source cell came from
1764  EntityHandle covCell = affectedCovCellFromID[sourceParentID];
1765  MB_CHK_ERR( m_interface->tag_get_data( sendProcTag, &covCell, 1, &orgTask ) );
1766  // put the overlap cell in corresponding range (set<EntityHandle>)
1767  overlapCellsForTask[orgTask].insert( intxCell );
1768  // also put it in this range, for debugging mostly
1769  overlapCellsToSend.insert( intxCell );
1770  }
1771  }
1772 
1773  // now prepare to send; will use crystal router, as the buffers in ParComm are prepared only
1774  // for neighbors; coverage mesh was also migrated with crystal router, so here we go again :(
1775  // find out the maximum number of edges of the polygons needed to be sent
1776  // we could we conservative and use a big number, or the number from intx, if we store it then?
1777  int maxEdges = 0;
1778  for( std::set< EntityHandle >::iterator it = overlapCellsToSend.begin(); it != overlapCellsToSend.end(); ++it )
1779  {
1780  EntityHandle intxCell = *it;
1781  int nnodes;
1782  const EntityHandle* conn;
1783  MB_CHK_ERR( m_interface->get_connectivity( intxCell, conn, nnodes ) );
1784  if( maxEdges < nnodes ) maxEdges = nnodes;
1785  }
1786 
1787  // find the maximum among processes in intersection
1788  int globalMaxEdges;
1789  if( m_pcomm )
1790  MPI_Allreduce( &maxEdges, &globalMaxEdges, 1, MPI_INT, MPI_MAX, m_pcomm->comm() );
1791  else
1792  globalMaxEdges = maxEdges;
1793 
1794 #ifdef MOAB_DBG
1795  if( is_root ) std::cout << "maximum number of edges for polygons to send is " << globalMaxEdges << "\n";
1796 #endif
1797 
1798 #ifdef MOAB_DBG
1799  EntityHandle tmpSet2;
1800  MB_CHK_SET_ERR( m_interface->create_meshset( MESHSET_SET, tmpSet2 ), "Can't create temporary set2" );
1801  // add the affected source and overlap elements
1802  for( std::set< EntityHandle >::iterator it = overlapCellsToSend.begin(); it != overlapCellsToSend.end(); ++it )
1803  {
1804  EntityHandle intxCell = *it;
1805  MB_CHK_SET_ERR( m_interface->add_entities( tmpSet2, &intxCell, 1 ), "Can't add entities" );
1806  }
1807  for( std::set< EntityHandle >::iterator it = affectedCovCells.begin(); it != affectedCovCells.end(); ++it )
1808  {
1809  EntityHandle covCell = *it;
1810  MB_CHK_SET_ERR( m_interface->add_entities( tmpSet2, &covCell, 1 ), "Can't add entities" );
1811  }
1812  std::stringstream ffs2;
1813  // these will contain coverage cells and intx cells on the boundary
1814  ffs2 << "affectedCells_" << m_pcomm->rank() << ".h5m";
1815  MB_CHK_ERR( m_interface->write_mesh( ffs2.str().c_str(), &tmpSet2, 1 ) );
1816 #endif
1817  // form tuple lists to send vertices and cells;
1818  // the problem is that the lists of vertices will need to have other information, like the
1819  // processor it comes from, and its index in that list; we may have to duplicate vertices, but
1820  // we do not care much; we will not duplicate overlap elements, just the vertices, as they may
1821  // come from different cells and different processes each vertex will have a local index and a
1822  // processor task it is coming from
1823 
1824  // look through the std::set's to be sent to other processes, and form the vertex tuples and
1825  // cell tuples
1826  //
1827  std::map< int, std::set< EntityHandle > > verticesOverlapForTask;
1828  // Range allVerticesToSend;
1829  std::set< EntityHandle > allVerticesToSend;
1830  std::map< EntityHandle, int > allVerticesToSendMap;
1831  int numVerts = 0;
1832  int numOverlapCells = 0;
1833  for( std::map< int, std::set< EntityHandle > >::iterator it = overlapCellsForTask.begin();
1834  it != overlapCellsForTask.end(); ++it )
1835  {
1836  int sendToProc = it->first;
1837  const std::set< EntityHandle >& overlapCellsToSend2 = it->second; // organize vertices in std::set per processor
1838  // Range vertices;
1839  std::set< EntityHandle > vertices; // collect all vertices connected to overlapCellsToSend2
1840  for( std::set< EntityHandle >::iterator set_it = overlapCellsToSend2.begin();
1841  set_it != overlapCellsToSend2.end(); ++set_it )
1842  {
1843  int nnodes_local = 0;
1844  const EntityHandle* conn1 = nullptr;
1845  MB_CHK_ERR( m_interface->get_connectivity( *set_it, conn1, nnodes_local ) );
1846  for( int k = 0; k < nnodes_local; k++ )
1847  vertices.insert( conn1[k] );
1848  }
1849  verticesOverlapForTask[sendToProc] = vertices;
1850  numVerts += (int)vertices.size();
1851  numOverlapCells += (int)overlapCellsToSend2.size();
1852  allVerticesToSend.insert( vertices.begin(), vertices.end() );
1853  }
1854  // build the index map, from entity handle to index in all vert set
1855  int j = 0;
1856  for( std::set< EntityHandle >::iterator vert_it = allVerticesToSend.begin(); vert_it != allVerticesToSend.end();
1857  ++vert_it, ++j )
1858  {
1859  EntityHandle vert = *vert_it;
1860  allVerticesToSendMap[vert] = j;
1861  }
1862 
1863  // first send vertices in a tuple list, then send overlap cells, according to requests
1864  // overlap cells need to send info about the blue and red parent tags, too
1865  TupleList TLv; //
1866  TLv.initialize( 2, 0, 0, 3, numVerts ); // to proc, index in all range, DP points
1867  TLv.enableWriteAccess();
1868 
1869  for( std::map< int, std::set< EntityHandle > >::iterator it = verticesOverlapForTask.begin();
1870  it != verticesOverlapForTask.end(); ++it )
1871  {
1872  int sendToProc = it->first;
1873  const std::set< EntityHandle >& vertices = it->second;
1874  int i = 0;
1875  for( std::set< EntityHandle >::iterator it2 = vertices.begin(); it2 != vertices.end(); ++it2, ++i )
1876  {
1877  int n = TLv.get_n();
1878  TLv.vi_wr[2 * n] = sendToProc; // send to processor
1879  EntityHandle v = *it2;
1880  int indexInAllVert = allVerticesToSendMap[v];
1881  TLv.vi_wr[2 * n + 1] = indexInAllVert; // will be orgProc, to differentiate indices
1882  // of vertices sent to "sentToProc"
1883  double coords[3];
1884  MB_CHK_ERR( m_interface->get_coords( &v, 1, coords ) );
1885  TLv.vr_wr[3 * n] = coords[0]; // departure position, of the node local_verts[i]
1886  TLv.vr_wr[3 * n + 1] = coords[1];
1887  TLv.vr_wr[3 * n + 2] = coords[2];
1888  TLv.inc_n();
1889  }
1890  }
1891 
1892  TupleList TLc;
1893  int sizeTuple = 4 + globalMaxEdges;
1894  // total number of overlap cells to send
1895  TLc.initialize( sizeTuple, 0, 0, 0,
1896  numOverlapCells ); // to proc, blue parent ID, red parent ID, nvert,
1897  // connectivity[globalMaxEdges] (global ID v), local eh)
1898  TLc.enableWriteAccess();
1899 
1900  for( std::map< int, std::set< EntityHandle > >::iterator it = overlapCellsForTask.begin();
1901  it != overlapCellsForTask.end(); ++it )
1902  {
1903  int sendToProc = it->first;
1904  const std::set< EntityHandle >& overlapCellsToSend2 = it->second;
1905  // send also the target and source parents for these overlap cells
1906  for( std::set< EntityHandle >::const_iterator it2 = overlapCellsToSend2.begin();
1907  it2 != overlapCellsToSend2.end(); ++it2 )
1908  {
1909  EntityHandle intxCell = *it2;
1910  int sourceParentID, targetParentID;
1911  MB_CHK_ERR( m_interface->tag_get_data( targetParentTag, &intxCell, 1, &targetParentID ) );
1912  MB_CHK_ERR( m_interface->tag_get_data( sourceParentTag, &intxCell, 1, &sourceParentID ) );
1913  int n = TLc.get_n();
1914  TLc.vi_wr[sizeTuple * n] = sendToProc;
1915  TLc.vi_wr[sizeTuple * n + 1] = sourceParentID;
1916  TLc.vi_wr[sizeTuple * n + 2] = targetParentID;
1917  int nnodes;
1918  const EntityHandle* conn = nullptr;
1919  MB_CHK_ERR( m_interface->get_connectivity( intxCell, conn, nnodes ) );
1920  TLc.vi_wr[sizeTuple * n + 3] = nnodes;
1921  for( int i = 0; i < nnodes; i++ )
1922  {
1923  int indexVertex = allVerticesToSendMap[conn[i]];
1924  ; // the vertex index will be now unique per original proc
1925  if( -1 == indexVertex ) MB_CHK_SET_ERR( MB_FAILURE, "Can't find vertex in range of vertices to send" );
1926  TLc.vi_wr[sizeTuple * n + 4 + i] = indexVertex;
1927  }
1928  // fill the rest with 0, just because we do not like uninitialized data
1929  for( int i = nnodes; i < globalMaxEdges; i++ )
1930  TLc.vi_wr[sizeTuple * n + 4 + i] = 0;
1931 
1932  TLc.inc_n();
1933  }
1934  }
1935 
1936  // send first the vertices and overlap cells to original task for coverage cells
1937  // now we are done populating the tuples; route them to the appropriate processors
1938 #ifdef MOAB_DBG
1939  std::stringstream ff1;
1940  ff1 << "TLc_" << rank << ".txt";
1941  TLc.print_to_file( ff1.str().c_str() );
1942  std::stringstream ffv;
1943  ffv << "TLv_" << rank << ".txt";
1944  TLv.print_to_file( ffv.str().c_str() );
1945 #endif
1946  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
1947  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLc, 0 );
1948 
1949 #ifdef MOAB_DBG
1950  TLc.print_to_file( ff1.str().c_str() ); // will append to existing file
1951  TLv.print_to_file( ffv.str().c_str() );
1952 #endif
1953  // first phase of transfer complete
1954  // now look at TLc, and sort by the source parent (index 1)
1955 
1956  TupleList::buffer buffer;
1957  buffer.buffer_init( sizeTuple * TLc.get_n() * 2 ); // allocate memory for sorting !! double
1958  TLc.sort( 1, &buffer );
1959 #ifdef MOAB_DBG
1960  TLc.print_to_file( ff1.str().c_str() );
1961 #endif
1962 
1963  // will keep a map with vertices per processor that will need to be used in TLv2;
1964  // so, availVertexIndicesPerProcessor[proc] is a map from vertex indices that are available from
1965  // this processor to the index in the local TLv; the used vertices will have to be sent to the
1966  // tasks that need them
1967 
1968  // connectivity of a cell is given by sending proc and index in original list of vertices from
1969  // that proc
1970 
1971  std::map< int, std::map< int, int > > availVertexIndicesPerProcessor;
1972  int nv = TLv.get_n();
1973  for( int i = 0; i < nv; i++ )
1974  {
1975  // int proc=TLv.vi_rd[3*i]; // it is coming from this processor
1976  int orgProc = TLv.vi_rd[2 * i]; // this is the original processor, for index vertex consideration
1977  int indexVert = TLv.vi_rd[2 * i + 1];
1978  availVertexIndicesPerProcessor[orgProc][indexVert] = i;
1979  }
1980 
1981  // now we have sorted the incoming overlap elements by the source element;
1982  // if we have overlap elements for one source coming from 2 or more processes, we need to send
1983  // back to the processes that do not have that overlap cell;
1984 
1985  // form new TLc2, TLv2, that will be distributed to necessary processes
1986  // first count source elements that are "spread" over multiple processes
1987  // TLc is ordered now by source ID; loop over them
1988  int n = TLc.get_n(); // total number of overlap elements received on current task;
1989  std::map< int, int > currentProcsCount;
1990  // form a map from proc to sets of vertex indices that will be sent using TLv2
1991  // will form a map between a source cell ID and tasks/targets that are partially overlapped by
1992  // these sources
1993  std::map< int, std::set< int > > sourcesForTasks;
1994  int sizeOfTLc2 = 0; // only increase when we will have to send data
1995  if( n > 0 )
1996  {
1997  int currentSourceID = TLc.vi_rd[sizeTuple * 0 + 1]; // we have written sizeTuple*0 for "clarity"
1998  int proc0 = TLc.vi_rd[sizeTuple * 0];
1999  currentProcsCount[proc0] = 1; //
2000 
2001  for( int i = 1; i < n; i++ )
2002  {
2003  int proc = TLc.vi_rd[sizeTuple * i];
2004  int sourceID = TLc.vi_rd[sizeTuple * i + 1];
2005  if( sourceID == currentSourceID )
2006  {
2007  if( currentProcsCount.find( proc ) == currentProcsCount.end() )
2008  {
2009  currentProcsCount[proc] = 1;
2010  }
2011  else
2012  currentProcsCount[proc]++;
2013  }
2014  if( sourceID != currentSourceID || ( ( n - 1 ) == i ) ) // we study the current source if we reach the last
2015  {
2016  // we have found a new source id, need to reset the proc counts, and establish if we
2017  // need to send data
2018  if( currentProcsCount.size() > 1 )
2019  {
2020 #ifdef VERBOSE
2021  std::cout << " source element " << currentSourceID << " intersects with "
2022  << currentProcsCount.size() << " target partitions\n";
2023  for( std::map< int, int >::iterator it = currentProcsCount.begin(); it != currentProcsCount.end();
2024  ++it )
2025  {
2026  int procID = it->first;
2027  int numOverCells = it->second;
2028  std::cout << " task:" << procID << " " << numOverCells << " cells\n";
2029  }
2030 
2031 #endif
2032  // estimate what we need to send
2033  for( std::map< int, int >::iterator it1 = currentProcsCount.begin(); it1 != currentProcsCount.end();
2034  ++it1 )
2035  {
2036  int proc1 = it1->first;
2037  sourcesForTasks[currentSourceID].insert( proc1 );
2038  for( std::map< int, int >::iterator it2 = currentProcsCount.begin();
2039  it2 != currentProcsCount.end(); ++it2 )
2040  {
2041  int proc2 = it2->first;
2042  if( proc1 != proc2 ) sizeOfTLc2 += it2->second;
2043  }
2044  }
2045  // mark vertices in TLv tuple that need to be sent
2046  }
2047  if( sourceID != currentSourceID ) // maybe we are not at the end, so continue on
2048  {
2049  currentSourceID = sourceID;
2050  currentProcsCount.clear();
2051  currentProcsCount[proc] = 1;
2052  }
2053  }
2054  }
2055  }
2056  // begin a loop to send the needed cells to the processes; also mark the vertices that need to
2057  // be sent, put them in a set
2058 
2059 #ifdef MOAB_DBG
2060  std::cout << " need to initialize TLc2 with " << sizeOfTLc2 << " cells\n ";
2061 #endif
2062 
2063  TupleList TLc2;
2064  int sizeTuple2 = 5 + globalMaxEdges; // send to, original proc for intx cell, source parent id,
2065  // target parent id,
2066  // number of vertices, then connectivity in terms of indices in vertex lists from original proc
2067  TLc2.initialize( sizeTuple2, 0, 0, 0, sizeOfTLc2 );
2068  TLc2.enableWriteAccess();
2069  // loop again through TLc, and select intx cells that have the problem sources;
2070 
2071  std::map< int, std::set< int > > verticesToSendForProc; // will look at indices in the TLv list
2072  // will form for each processor, the index list from TLv
2073  for( int i = 0; i < n; i++ )
2074  {
2075  int sourceID = TLc.vi_rd[sizeTuple * i + 1];
2076  if( sourcesForTasks.find( sourceID ) != sourcesForTasks.end() )
2077  {
2078  // it means this intx cell needs to be sent to any proc that is not "original" to it
2079  std::set< int > procs = sourcesForTasks[sourceID]; // set of processors involved with this source
2080  if( procs.size() < 2 ) MB_CHK_SET_ERR( MB_FAILURE, " not enough processes involved with a sourceID cell" );
2081 
2082  int orgProc = TLc.vi_rd[sizeTuple * i]; // this intx cell was sent from this orgProc, originally
2083  // will need to be sent to all other procs from above set; also, need to mark the vertex
2084  // indices for that proc, and check that they are available to populate TLv2
2085  std::map< int, int >& availableVerticesFromThisProc = availVertexIndicesPerProcessor[orgProc];
2086  for( std::set< int >::iterator setIt = procs.begin(); setIt != procs.end(); ++setIt )
2087  {
2088  int procID = *setIt;
2089  // send this cell to the other processors, not to orgProc this cell is coming from
2090 
2091  if( procID != orgProc )
2092  {
2093  // send the cell to this processor;
2094  int n2 = TLc2.get_n();
2095  if( n2 >= sizeOfTLc2 ) MB_CHK_SET_ERR( MB_FAILURE, " memory overflow" );
2096  //
2097  std::set< int >& indexVerticesInTLv = verticesToSendForProc[procID];
2098  TLc2.vi_wr[n2 * sizeTuple2] = procID; // send to
2099  TLc2.vi_wr[n2 * sizeTuple2 + 1] = orgProc; // this cell is coming from here
2100  TLc2.vi_wr[n2 * sizeTuple2 + 2] = sourceID; // source parent of the intx cell
2101  TLc2.vi_wr[n2 * sizeTuple2 + 3] = TLc.vi_rd[sizeTuple * i + 2]; // target parent of the intx cell
2102  // number of vertices of the intx cell
2103  int nvert = TLc.vi_rd[sizeTuple * i + 3];
2104  TLc2.vi_wr[n2 * sizeTuple2 + 4] = nvert;
2105  // now loop through the connectivity, and make sure the vertices are available;
2106  // mark them, to populate later the TLv2 tuple list
2107 
2108  // just copy the vertices, including 0 ones
2109  for( int j = 0; j < nvert; j++ )
2110  {
2111  int vertexIndex = TLc.vi_rd[i * sizeTuple + 4 + j];
2112  // is this vertex available from org proc?
2113  if( availableVerticesFromThisProc.find( vertexIndex ) == availableVerticesFromThisProc.end() )
2114  {
2115  MB_CHK_SET_ERR( MB_FAILURE, " vertex index not available from processor" );
2116  }
2117  TLc2.vi_wr[n2 * sizeTuple2 + 5 + j] = vertexIndex;
2118  int indexInTLv = availVertexIndicesPerProcessor[orgProc][vertexIndex];
2119  indexVerticesInTLv.insert( indexInTLv );
2120  }
2121 
2122  for( int j = nvert; j < globalMaxEdges; j++ )
2123  {
2124  TLc2.vi_wr[n2 * sizeTuple2 + 5 + j] = 0; // or mark them 0
2125  }
2126  TLc2.inc_n();
2127  }
2128  }
2129  }
2130  }
2131 
2132  // now we have to populate TLv2, with original source proc, index of vertex, and coordinates
2133  // from TLv use the verticesToSendForProc sets from above, and the map from index in proc to the
2134  // index in TLv
2135  TupleList TLv2;
2136  int numVerts2 = 0;
2137  // how many vertices to send?
2138  for( std::map< int, std::set< int > >::iterator it = verticesToSendForProc.begin();
2139  it != verticesToSendForProc.end(); ++it )
2140  {
2141  const std::set< int >& indexInTLvSet = it->second;
2142  numVerts2 += (int)indexInTLvSet.size();
2143  }
2144  TLv2.initialize( 3, 0, 0, 3,
2145  numVerts2 ); // send to, original proc, index in original proc, and 3 coords
2146  TLv2.enableWriteAccess();
2147  for( std::map< int, std::set< int > >::iterator it = verticesToSendForProc.begin();
2148  it != verticesToSendForProc.end(); ++it )
2149  {
2150  int sendToProc = it->first;
2151  const std::set< int >& indexInTLvSet = it->second;
2152  // now, look at indices in TLv, to find out the original proc, and the index in that list
2153  for( std::set< int >::iterator itSet = indexInTLvSet.begin(); itSet != indexInTLvSet.end(); ++itSet )
2154  {
2155  int indexInTLv = *itSet;
2156  int orgProc = TLv.vi_rd[2 * indexInTLv];
2157  int indexVertexInOrgProc = TLv.vi_rd[2 * indexInTLv + 1];
2158  int nv2 = TLv2.get_n();
2159  TLv2.vi_wr[3 * nv2] = sendToProc;
2160  TLv2.vi_wr[3 * nv2 + 1] = orgProc;
2161  TLv2.vi_wr[3 * nv2 + 2] = indexVertexInOrgProc;
2162  for( int j = 0; j < 3; j++ )
2163  TLv2.vr_wr[3 * nv2 + j] =
2164  TLv.vr_rd[3 * indexInTLv + j]; // departure position, of the node local_verts[i]
2165  TLv2.inc_n();
2166  }
2167  }
2168  // now, finally, transfer the vertices and the intx cells;
2169  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv2, 0 );
2170  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLc2, 0 );
2171  // now, look at vertices from TLv2, and create them
2172  // we should have in TLv2 only vertices with orgProc different from current task
2173 #ifdef MOAB_DBG
2174  std::stringstream ff2;
2175  ff2 << "TLc2_" << rank << ".txt";
2176  TLc2.print_to_file( ff2.str().c_str() );
2177  std::stringstream ffv2;
2178  ffv2 << "TLv2_" << rank << ".txt";
2179  TLv2.print_to_file( ffv2.str().c_str() );
2180 #endif
2181  // first create vertices, and make a map from origin processor, and index, to entity handle
2182  // (index in TLv2 )
2183  Tag ghostTag;
2184  int orig_proc = -1;
2185  MB_CHK_ERR( m_interface->tag_get_handle( "ORIG_PROC", 1, MB_TYPE_INTEGER, ghostTag, MB_TAG_DENSE | MB_TAG_CREAT,
2186  &orig_proc ) );
2187 
2188  int nvNew = TLv2.get_n();
2189  // if number of vertices to be created is 0, it means there is no need of ghost intx cells,
2190  // because everything matched perfectly (it can happen in manufactured cases)
2191  if( 0 == nvNew ) return MB_SUCCESS;
2192  // create a vertex h for each coordinate
2193  Range newVerts;
2194  MB_CHK_ERR( m_interface->create_vertices( &( TLv2.vr_rd[0] ), nvNew, newVerts ) );
2195  // now create a map from index , org proc, to actual entity handle corresponding to it
2196  std::map< int, std::map< int, EntityHandle > > vertexPerProcAndIndex;
2197  for( int i = 0; i < nvNew; i++ )
2198  {
2199  int orgProc = TLv2.vi_rd[3 * i + 1];
2200  int indexInVert = TLv2.vi_rd[3 * i + 2];
2201  vertexPerProcAndIndex[orgProc][indexInVert] = newVerts[i];
2202  }
2203 
2204  // new polygons will receive a dense tag, with default value -1, with the processor task they
2205  // originally belonged to
2206 
2207  // now form the needed cells, in order
2208  Range newPolygons;
2209  int ne = TLc2.get_n();
2210  for( int i = 0; i < ne; i++ )
2211  {
2212  int orgProc = TLc2.vi_rd[i * sizeTuple2 + 1]; // this cell is coming from here, originally
2213  int sourceID = TLc2.vi_rd[i * sizeTuple2 + 2]; // source parent of the intx cell
2214  int targetID = TLc2.vi_wr[i * sizeTuple2 + 3]; // target parent of intx cell
2215  int nve = TLc2.vi_wr[i * sizeTuple2 + 4]; // number of vertices for the polygon
2216  std::vector< EntityHandle > conn;
2217  conn.resize( nve );
2218  for( int j = 0; j < nve; j++ )
2219  {
2220  int indexV = TLc2.vi_wr[i * sizeTuple2 + 5 + j];
2221  EntityHandle vh = vertexPerProcAndIndex[orgProc][indexV];
2222  conn[j] = vh;
2223  }
2224  EntityHandle polyNew;
2225  MB_CHK_ERR( m_interface->create_element( MBPOLYGON, &conn[0], nve, polyNew ) );
2226  newPolygons.insert( polyNew );
2227  MB_CHK_ERR( m_interface->tag_set_data( targetParentTag, &polyNew, 1, &targetID ) );
2228  MB_CHK_ERR( m_interface->tag_set_data( sourceParentTag, &polyNew, 1, &sourceID ) );
2229  MB_CHK_ERR( m_interface->tag_set_data( ghostTag, &polyNew, 1, &orgProc ) );
2230  }
2231 
2232 #ifdef MOAB_DBG
2233  EntityHandle tmpSet3;
2234  MB_CHK_SET_ERR( m_interface->create_meshset( MESHSET_SET, tmpSet3 ), "Can't create temporary set3" );
2235  // add the boundary set and edges, and save it to a file
2236  MB_CHK_SET_ERR( m_interface->add_entities( tmpSet3, newPolygons ), "Can't add entities" );
2237 
2238  std::stringstream ffs4;
2239  ffs4 << "extraIntxCells" << rank << ".h5m";
2240  MB_CHK_ERR( m_interface->write_mesh( ffs4.str().c_str(), &tmpSet3, 1 ) );
2241 #endif
2242 
2243  // add the new polygons to the overlap set
2244  // these will be ghosted, so will participate in conservation only
2245  MB_CHK_ERR( m_interface->add_entities( m_overlap_set, newPolygons ) );
2246  return MB_SUCCESS;
2247 }
2248 #endif
2249 
2250 #undef MOAB_DBG
2251 
2252 ErrorCode TempestRemapper::GetIMasks( Remapper::IntersectionContext ctx, std::vector< int >& masks )
2253 {
2254  Tag maskTag;
2255  // it should have been created already, if not, we might have a problem
2256  int def_val = 1;
2257  MB_CHK_SET_ERR( m_interface->tag_get_handle( "GRID_IMASK", 1, MB_TYPE_INTEGER, maskTag, MB_TAG_DENSE | MB_TAG_CREAT,
2258  &def_val ),
2259  "Trouble creating GRID_IMASK tag" );
2260 
2261  switch( ctx )
2262  {
2263  case Remapper::SourceMesh: {
2264  if( point_cloud_source )
2265  {
2266  masks.resize( m_source_vertices.size() );
2267  MB_CHK_SET_ERR( m_interface->tag_get_data( maskTag, m_source_vertices, &masks[0] ),
2268  "Trouble getting GRID_IMASK tag" );
2269  }
2270  else
2271  {
2272  masks.resize( m_source_entities.size() );
2273  MB_CHK_SET_ERR( m_interface->tag_get_data( maskTag, m_source_entities, &masks[0] ),
2274  "Trouble getting GRID_IMASK tag" );
2275  }
2276  return MB_SUCCESS;
2277  }
2278  case Remapper::TargetMesh: {
2279  if( point_cloud_target )
2280  {
2281  masks.resize( m_target_vertices.size() );
2282  MB_CHK_SET_ERR( m_interface->tag_get_data( maskTag, m_target_vertices, &masks[0] ),
2283  "Trouble getting GRID_IMASK tag" );
2284  }
2285  else
2286  {
2287  masks.resize( m_target_entities.size() );
2288  MB_CHK_SET_ERR( m_interface->tag_get_data( maskTag, m_target_entities, &masks[0] ),
2289  "Trouble getting GRID_IMASK tag" );
2290  }
2291  return MB_SUCCESS;
2292  }
2293  case Remapper::CoveringMesh:
2294  case Remapper::OverlapMesh:
2295  default:
2296  return MB_SUCCESS;
2297  }
2298 }
2299 
2300 } // namespace moab