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