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