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