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