40 #ifdef MOAB_HAVE_TEMPESTREMAP
42 #include "FiniteElementTools.h"
43 #include "GaussLobattoQuadrature.h"
56 if( initialize_fsets )
75 MPI_Initialized( &flagInit );
78 assert( m_pcomm !=
nullptr );
79 rank = m_pcomm->rank();
80 size = m_pcomm->size();
85 AnnounceOnlyOutputOnRankZero();
153 std::string inputFilename,
173 MB_CHK_SET_ERR( MB_FAILURE,
"Invalid IntersectionContext context provided" );
180 if( outputEnabled ) std::cout <<
"\nLoading TempestRemap Mesh object from file = " << inputFilename <<
" ...\n";
183 NcError
error( NcError::silent_nonfatal );
188 if( outputEnabled ) std::cout <<
"Loading mesh ...\n";
189 Mesh* mesh =
new Mesh( inputFilename );
190 mesh->RemoveZeroEdges();
191 if( outputEnabled ) std::cout <<
"----------------\n";
196 if( outputEnabled ) std::cout <<
"Validating mesh ...\n";
198 if( outputEnabled ) std::cout <<
"-------------------\n";
204 if( outputEnabled ) std::cout <<
"Constructing edge map on mesh ...\n";
205 mesh->ConstructEdgeMap(
false );
206 if( outputEnabled ) std::cout <<
"---------------------------------\n";
209 if( tempest_mesh ) *tempest_mesh = mesh;
211 catch( Exception& e )
213 std::cout <<
"TempestRemap ERROR: " << e.ToString() <<
"\n";
231 if( outputEnabled ) std::cout <<
"Converting (source) TempestRemap Mesh object to MOAB representation ...\n";
237 if( outputEnabled ) std::cout <<
"Converting (target) TempestRemap Mesh object to MOAB representation ...\n";
243 if( outputEnabled ) std::cout <<
"Converting (overlap) TempestRemap Mesh object to MOAB representation ...\n";
248 MB_CHK_SET_ERR( MB_FAILURE,
"Invalid IntersectionContext context provided" );
252 #define NEW_CONVERT_LOGIC
254 #ifdef NEW_CONVERT_LOGIC
263 const NodeVector& nodes = mesh->nodes;
264 const FaceVector& faces = mesh->faces;
275 std::vector< double* > arrays;
276 std::vector< int > gidsv( nodes.size() );
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 )
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;
287 Range mbverts( startv, startv + nodes.size() - 1 );
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 );
299 if( storeParentInfo )
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 );
322 for(
unsigned ifaces = 0; ifaces < faces.size(); ++ifaces )
324 const Face&
face = faces[ifaces];
325 const unsigned num_v_per_elem =
face.edges.size();
327 for(
unsigned iedges = 0; iedges < num_v_per_elem; ++iedges )
329 conn[iedges] = startv +
face.edges[iedges].node[0];
332 switch( num_v_per_elem )
355 gidse[ifaces] = ifaces + 1;
357 if( storeParentInfo )
359 srcParent[ifaces] = mesh->vecSourceFaceIx[ifaces] + 1;
360 tgtParent[ifaces] = mesh->vecTargetFaceIx[ifaces] + 1;
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 );
371 if( storeParentInfo )
381 if( vertices ) *vertices = mbverts;
396 const NodeVector& nodes = mesh->nodes;
397 const FaceVector& faces = mesh->faces;
400 dbgprint.set_prefix(
"[TempestToMOAB]: " );
408 std::vector< double* > arrays;
409 std::vector< int > gidsv( nodes.size() );
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 )
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;
420 Range mbverts( startv, startv + nodes.size() - 1 );
427 Tag srcParentTag, tgtParentTag;
428 std::vector< int > srcParent, tgtParent;
429 bool storeParentInfo = ( mesh->vecSourceFaceIx.size() > 0 );
431 if( storeParentInfo )
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 )
455 const int iType = faces[ifaces].edges.size();
457 typeNSeqs[iType].push_back( ifaces );
460 for(
unsigned iType = 0; iType < NMAXPOLYEDGES; ++iType )
462 if( !nPolys[iType] )
continue;
464 const unsigned num_v_per_elem = iType;
469 switch( num_v_per_elem )
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" );
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" );
483 dbgprint.printf( 0,
"....Block %d: Polygonal [%u] Elements [%u].\n", iBlock++, 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" );
489 Range mbcells( starte, starte + nPolys[iType] - 1 );
492 if( storeParentInfo )
494 srcParent.resize( mbcells.size(), -1 );
495 tgtParent.resize( mbcells.size(), -1 );
498 std::vector< int > gids( typeNSeqs[iType].
size() );
499 for(
unsigned ifaces = 0, offset = 0; ifaces < typeNSeqs[iType].size(); ++ifaces )
501 const int fIndex = typeNSeqs[iType][ifaces];
502 const Face&
face = faces[fIndex];
504 for(
unsigned iedges = 0; iedges <
face.edges.size(); ++iedges )
506 conn[offset++] = startv +
face.edges[iedges].node[0];
509 if( storeParentInfo )
511 srcParent[ifaces] = mesh->vecSourceFaceIx[fIndex] + 1;
512 tgtParent[ifaces] = mesh->vecTargetFaceIx[fIndex] + 1;
515 gids[ifaces] = typeNSeqs[iType][ifaces] + 1;
522 rval =
iface->update_adjacencies( starte, nPolys[iType], num_v_per_elem, conn );
MB_CHK_SET_ERR( rval,
"Can't update adjacencies" );
528 if( storeParentInfo )
537 if( vertices ) *vertices = mbverts;
553 "Can't convert source mesh to TempestRemap format" );
559 "Can't convert source coverage mesh to TempestRemap format" );
565 "Can't convert target mesh to TempestRemap format" );
582 if( outputEnabled )
dbgprint.printf( 0,
"Converting (source) MOAB to TempestRemap Mesh representation ...\n" );
585 "Can't convert source mesh to Tempest" );
595 dbgprint.printf( 0,
"Converting (covering source) MOAB to TempestRemap Mesh representation ...\n" );
598 "Can't convert convering source mesh to TempestRemap format" );
603 if( outputEnabled )
dbgprint.printf( 0,
"Converting (target) MOAB to TempestRemap Mesh representation ...\n" );
606 "Can't convert target mesh to Tempest" );
612 if( outputEnabled )
dbgprint.printf( 0,
"Converting (overlap) MOAB to TempestRemap Mesh representation ...\n" );
617 MB_CHK_SET_ERR( MB_FAILURE,
"Invalid IntersectionContext context provided" );
631 NodeVector& nodes = mesh->nodes;
632 FaceVector& faces = mesh->faces;
637 const size_t nelems = elems.
size();
640 faces.resize( nelems );
644 if( verts.
size() == 0 )
650 std::map< EntityHandle, int > indxMap;
651 bool useRange =
true;
660 std::vector< int > globIds( nelems );
663 std::vector< size_t > sortedIdx;
666 sortedIdx.resize( nelems );
668 std::iota( sortedIdx.begin(), sortedIdx.end(), 0 );
671 std::sort( sortedIdx.begin(), sortedIdx.end(),
672 [&globIds](
size_t i1,
size_t i2 ) { return globIds[i1] < globIds[i2]; } );
685 while( connectface[nnodesf - 2] == connectface[nnodesf - 1] && nnodesf > 3 )
688 face.edges.resize( nnodesf );
689 for(
int iverts = 0; iverts < nnodesf; ++iverts )
691 int indx = ( useRange ? verts.
index( connectface[iverts] ) : indxMap[connectface[iverts]] );
693 face.SetNode( iverts, indx );
697 unsigned nnodes = verts.
size();
698 nodes.resize( nnodes );
701 std::vector< double > coordx( nnodes ), coordy( nnodes ), coordz( nnodes );
703 for(
unsigned inode = 0; inode < nnodes; ++inode )
705 Node& node = nodes[inode];
706 node.x = coordx[inode];
707 node.y = coordy[inode];
708 node.z = coordz[inode];
714 mesh->RemoveCoincidentNodes();
715 mesh->RemoveZeroEdges();
737 if( std::get< 1 >( a ) == std::get< 1 >( b ) )
738 return std::get< 2 >( a ) < std::get< 2 >( b );
740 return std::get< 1 >( a ) < std::get< 1 >( b );
745 sharedGhostEntities.
clear();
753 rval = m_interface->get_entities_by_dimension( m_overlap_set, 2, allents );
MB_CHK_SET_ERR( rval,
"Getting entities dim 2 failed" );
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 )
762 sharedents.
insert( allents[i] );
764 allents =
subtract( allents, sharedents );
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 );
776 sharedGhostEntities.
merge( sharedents );
794 std::vector< std::array< int, 3 > > sorted_overlap_order( n_overlap_entitites,
795 std::array< int, 3 >( { -1, -1, -1 } ) );
797 Tag srcParentTag, tgtParentTag;
801 m_overlap->vecTargetFaceIx.resize( n_overlap_entitites );
802 m_overlap->vecSourceFaceIx.resize( n_overlap_entitites );
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 );
819 std::vector< int > ghFlags;
823 ghFlags.resize( n_overlap_entitites );
829 std::vector< int > rbids_src( n_overlap_entitites ), rbids_tgt( n_overlap_entitites );
833 for(
size_t ix = 0; ix < n_overlap_entitites; ++ix )
835 std::get< 0 >( sorted_overlap_order[ix] ) = ix;
836 std::get< 1 >( sorted_overlap_order[ix] ) = find_lid( gids_src, rbids_src[ix] );
838 std::get< 2 >( sorted_overlap_order[ix] ) = -1;
840 std::get< 2 >( sorted_overlap_order[ix] ) = find_lid( gids_tgt, rbids_tgt[ix] );
844 std::sort( sorted_overlap_order.begin(), sorted_overlap_order.end(),
IntPairComparator );
846 for(
unsigned ie = 0; ie < n_overlap_entitites; ++ie )
849 m_overlap->vecSourceFaceIx[ie] = std::get< 1 >( sorted_overlap_order[ie] );
850 m_overlap->vecTargetFaceIx[ie] = std::get< 2 >( sorted_overlap_order[ie] );
855 faces.resize( n_overlap_entitites );
861 std::map< EntityHandle, int > indxMap;
870 const unsigned iface = std::get< 0 >( sorted_overlap_order[ifac] );
871 Face&
face = faces[ifac];
879 face.edges.resize( nnodesf );
880 for(
int iverts = 0; iverts < nnodesf; ++iverts )
882 int indx = indxMap[connectface[iverts]];
884 face.SetNode( iverts, indx );
888 sorted_overlap_order.clear();
890 unsigned nnodes = verts.
size();
892 nodes.resize( nnodes );
895 std::vector< double > coordx( nnodes ), coordy( nnodes ), coordz( nnodes );
897 for(
unsigned inode = 0; inode < nnodes; ++inode )
899 Node& node = nodes[inode];
900 node.x = coordx[inode];
901 node.y = coordy[inode];
902 node.z = coordz[inode];
923 const bool fAllParallel,
924 const bool fInputConcave,
925 const bool fOutputConcave )
930 if( is_root &&
size == 1 )
932 this->m_source->CalculateFaceAreas( fInputConcave );
933 this->m_target->CalculateFaceAreas( fOutputConcave );
934 this->m_overlap->Write( strOutputFileName.c_str(), NcFile::Netcdf4 );
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 );
950 this->m_source->CalculateFaceAreas( fInputConcave );
951 this->m_target->CalculateFaceAreas( fOutputConcave );
952 this->m_overlap->Write( strOutputFileName.c_str(), NcFile::Netcdf4 );
987 #ifndef MOAB_HAVE_MPI
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++;
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 );
1023 const std::string& dofTagName,
1026 const int csResolution = std::sqrt( ntot_elements / 6.0 );
1027 if( csResolution * csResolution * 6 != ntot_elements )
return MB_INVALID_SIZE;
1032 if( GenerateCSMesh( csMesh, csResolution,
"",
"NetCDF4" ) )
1034 "Failed to generate CS mesh through TempestRemap" );
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" );
1044 const int ntot_elements,
1047 const std::string dofTagName,
1053 bool created =
false;
1058 int nElements =
static_cast< int >( csMesh.faces.size() );
1063 DataArray3D< int > dataGLLnodes;
1064 dataGLLnodes.Allocate( nP, nP, nElements );
1066 std::map< Node, int > mapNodes;
1067 std::map< Node, moab::EntityHandle > mapLocalMBNodes;
1070 DataArray1D< double > dG;
1071 DataArray1D< double > dW;
1072 GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );
1075 if( secondary_ents )
entities.insert( secondary_ents->
begin(), secondary_ents->
end() );
1077 for(
unsigned iel = 0; iel <
entities.size(); ++iel )
1081 Node elCentroid( elcoords[0], elcoords[1], elcoords[2] );
1082 mapLocalMBNodes.insert( std::pair< Node, moab::EntityHandle >( elCentroid, eh ) );
1091 int* dofIDs =
new int[nP * nP];
1094 for(
int k = 0; k < nElements; k++ )
1096 const Face&
face = csMesh.faces[k];
1097 const NodeVector& nodes = csMesh.nodes;
1099 if(
face.edges.size() != 4 )
1101 _EXCEPTIONT(
"Mesh must only contain quadrilateral elements" );
1105 centroid.x = centroid.y = centroid.z = 0.0;
1106 for(
unsigned l = 0; l <
face.edges.size(); ++l )
1108 centroid.x += nodes[
face[l]].x;
1109 centroid.y += nodes[
face[l]].y;
1110 centroid.z += nodes[
face[l]].z;
1112 const double factor = 1.0 /
face.edges.size();
1113 centroid.x *= factor;
1114 centroid.y *= factor;
1115 centroid.z *= factor;
1117 bool locElem =
false;
1119 if( mapLocalMBNodes.find( centroid ) != mapLocalMBNodes.end() )
1122 current_eh = mapLocalMBNodes[centroid];
1125 for(
int j = 0; j < nP; j++ )
1127 for(
int i = 0; i < nP; i++ )
1142 const double& dAlpha = dG[i];
1143 const double& dBeta = dG[j];
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;
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;
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;
1158 double dR = sqrt( dXc * dXc + dYc * dYc + dZc * dZc );
1161 nodeGLL.x = dXc / dR;
1162 nodeGLL.y = dYc / dR;
1163 nodeGLL.z = dZc / dR;
1166 std::map< Node, int >::const_iterator iter = mapNodes.find( nodeGLL );
1167 if( iter == mapNodes.end() )
1170 int ixNode =
static_cast< int >( mapNodes.size() );
1171 mapNodes.insert( std::pair< Node, int >( nodeGLL, ixNode ) );
1172 dataGLLnodes[j][i][k] = ixNode + 1;
1176 dataGLLnodes[j][i][k] = iter->second + 1;
1179 dofIDs[j * nP + i] = dataGLLnodes[j][i][k];
1191 mapLocalMBNodes.clear();
1206 int nb_ghost_layers )
1220 #ifdef MOAB_HAVE_MPI
1221 mbintx->set_parallel_comm( m_pcomm );
1231 #ifdef MOAB_HAVE_MPI
1240 std::stringstream filename;
1241 filename <<
"covering" <<
rank <<
".h5m";
1243 std::stringstream targetFile;
1244 targetFile <<
"target" <<
rank <<
".h5m";
1255 double tolerance = 1e-6, btolerance = 1e-3;
1263 for(
unsigned ie = 0; ie < targetVerts.
size(); ++ie )
1280 std::vector< moab::EntityHandle > leaf_elems;
1285 if( !leaf_elems.size() )
1293 std::vector< double > centroids( leaf_elems.size() * 3 );
1298 for(
size_t il = 0; il < leaf_elems.size(); ++il )
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 );
1305 if( locdist < dist )
1316 <<
": [Error] - Could not find a minimum distance within the leaf "
1318 << dist << std::endl;
1334 #ifdef MOAB_HAVE_MPI
1347 const bool outputEnabled = ( this->
rank == 0 );
1349 dbgprint.
set_prefix(
"[ComputeOverlapMesh]: " );
1369 bool concaveMeshA =
false, concaveMeshB =
false;
1374 "exact", concaveMeshA, concaveMeshB,
true,
false ) )
1375 MB_CHK_SET_ERR( MB_FAILURE,
"TempestRemap: cannot compute the intersection of meshes on the sphere" );
1382 if( outputEnabled )
dbgprint.printf( 0,
"Computing intersection mesh with the Kd-tree search algorithm" );
1388 dbgprint.printf( 0,
"Computing intersection mesh with the advancing-front propagation algorithm" );
1392 #ifdef MOAB_HAVE_MPI
1396 std::stringstream ffc, fft, ffo;
1397 ffc <<
"cover_" <<
rank <<
".h5m";
1398 fft <<
"target_" <<
rank <<
".h5m";
1399 ffo <<
"intx_" <<
rank <<
".h5m";
1412 std::map< int, int > loc_gid_to_lid_covsrc;
1413 std::vector< int > gids( covEnts.
size(), -1 );
1418 for(
unsigned ie = 0; ie < gids.size(); ++ie )
1420 assert( gids[ie] > 0 );
1421 loc_gid_to_lid_covsrc[gids[ie]] = ie;
1424 Range intxCov, intxCells;
1434 assert( srcParent >= 0 );
1435 intxCov.
insert( covEnts[loc_gid_to_lid_covsrc[srcParent]] );
1458 rval =
mb->get_connectivity( eh, conn, num_nodes );
MB_CHK_ERR( rval );
1471 rval =
mb->get_connectivity( *it, conn, len,
false );
MB_CHK_ERR( rval );
1472 for(
int ie = 0; ie < len; ++ie )
1474 std::vector< EntityHandle > adjacent_entities;
1476 for(
auto ent : adjacent_entities )
1477 notNeededCovCells.
erase( ent );
1482 std::string(
"sourcecoveragemesh_p" + std::to_string(
rank ) +
".h5m" ).c_str(),
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()
1533 #ifdef MOAB_HAVE_MPI
1559 Range boundaryEdges;
1560 MB_CHK_ERR( skinner.find_skin( 0, this->m_target_entities,
false, boundaryEdges ) );
1564 Range boundaryCells;
1574 std::stringstream ffs;
1575 ffs <<
"boundaryCells_0" <<
rank <<
".h5m";
1582 std::set< int > targetBoundaryIds;
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 );
1596 std::set< int > affectedSourceCellsIds;
1597 Tag targetParentTag, sourceParentTag;
1600 for(
Range::iterator it = overlapCells.begin(); it != overlapCells.end(); it++ )
1603 int targetParentID, sourceParentID;
1605 if( targetBoundaryIds.find( targetParentID ) != targetBoundaryIds.end() )
1609 affectedSourceCellsIds.insert( sourceParentID );
1618 std::map< int, EntityHandle > affectedCovCellFromID;
1624 std::set< EntityHandle > affectedCovCells;
1629 for(
Range::iterator it = covCells.begin(); it != covCells.end(); it++ )
1634 if( affectedSourceCellsIds.find( covID ) != affectedSourceCellsIds.end() )
1637 affectedCovCellFromID[covID] = covCell;
1638 affectedCovCells.insert( covCell );
1649 std::map< int, std::set< EntityHandle > > overlapCellsForTask;
1652 std::set< EntityHandle > overlapCellsToSend;
1654 for(
Range::iterator it = overlapCells.begin(); it != overlapCells.end(); it++ )
1660 if( affectedSourceCellsIds.find( sourceParentID ) != affectedSourceCellsIds.end() )
1663 EntityHandle covCell = affectedCovCellFromID[sourceParentID];
1666 overlapCellsForTask[orgTask].insert( intxCell );
1668 overlapCellsToSend.insert( intxCell );
1677 for( std::set< EntityHandle >::iterator it = overlapCellsToSend.begin(); it != overlapCellsToSend.end(); it++ )
1683 if( maxEdges < nnodes ) maxEdges = nnodes;
1689 MPI_Allreduce( &maxEdges, &globalMaxEdges, 1, MPI_INT, MPI_MAX, m_pcomm->comm() );
1691 globalMaxEdges = maxEdges;
1694 if(
is_root ) std::cout <<
"maximum number of edges for polygons to send is " << globalMaxEdges <<
"\n";
1701 for( std::set< EntityHandle >::iterator it = overlapCellsToSend.begin(); it != overlapCellsToSend.end(); it++ )
1706 for( std::set< EntityHandle >::iterator it = affectedCovCells.begin(); it != affectedCovCells.end(); it++ )
1711 std::stringstream ffs2;
1713 ffs2 <<
"affectedCells_" << m_pcomm->rank() <<
".h5m";
1726 std::map< int, std::set< EntityHandle > > verticesOverlapForTask;
1728 std::set< EntityHandle > allVerticesToSend;
1729 std::map< EntityHandle, int > allVerticesToSendMap;
1731 int numOverlapCells = 0;
1732 for( std::map<
int, std::set< EntityHandle > >::iterator it = overlapCellsForTask.begin();
1733 it != overlapCellsForTask.end(); it++ )
1735 int sendToProc = it->first;
1736 std::set< EntityHandle >& overlapCellsToSend2 = it->second;
1738 std::set< EntityHandle > vertices;
1739 for( std::set< EntityHandle >::iterator set_it = overlapCellsToSend2.begin();
1740 set_it != overlapCellsToSend2.end(); ++set_it )
1742 int nnodes_local = 0;
1745 for(
int k = 0; k < nnodes_local; k++ )
1746 vertices.insert( conn1[k] );
1748 verticesOverlapForTask[sendToProc] = vertices;
1749 numVerts += (int)vertices.size();
1750 numOverlapCells += (int)overlapCellsToSend2.size();
1751 allVerticesToSend.insert( vertices.begin(), vertices.end() );
1755 for( std::set< EntityHandle >::iterator vert_it = allVerticesToSend.begin(); vert_it != allVerticesToSend.end();
1759 allVerticesToSendMap[vert] = j;
1765 TLv.initialize( 2, 0, 0, 3, numVerts );
1766 TLv.enableWriteAccess();
1768 for( std::map<
int, std::set< EntityHandle > >::iterator it = verticesOverlapForTask.begin();
1769 it != verticesOverlapForTask.end(); it++ )
1771 int sendToProc = it->first;
1772 std::set< EntityHandle >& vertices = it->second;
1774 for( std::set< EntityHandle >::iterator it2 = vertices.begin(); it2 != vertices.end(); it2++, i++ )
1776 int n = TLv.get_n();
1777 TLv.vi_wr[2 * n] = sendToProc;
1779 int indexInAllVert = allVerticesToSendMap[v];
1780 TLv.vi_wr[2 * n + 1] = indexInAllVert;
1784 TLv.vr_wr[3 * n] = coords[0];
1785 TLv.vr_wr[3 * n + 1] = coords[1];
1786 TLv.vr_wr[3 * n + 2] = coords[2];
1792 int sizeTuple = 4 + globalMaxEdges;
1794 TLc.initialize( sizeTuple, 0, 0, 0,
1797 TLc.enableWriteAccess();
1799 for( std::map<
int, std::set< EntityHandle > >::iterator it = overlapCellsForTask.begin();
1800 it != overlapCellsForTask.end(); it++ )
1802 int sendToProc = it->first;
1803 std::set< EntityHandle >& overlapCellsToSend2 = it->second;
1805 for( std::set< EntityHandle >::iterator it2 = overlapCellsToSend2.begin(); it2 != overlapCellsToSend2.end();
1809 int sourceParentID, targetParentID;
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;
1819 TLc.vi_wr[sizeTuple * n + 3] = nnodes;
1820 for(
int i = 0; i < nnodes; i++ )
1822 int indexVertex = allVerticesToSendMap[conn[i]];
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;
1828 for(
int i = nnodes; i < globalMaxEdges; i++ )
1829 TLc.vi_wr[sizeTuple * n + 4 + i] = 0;
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() );
1845 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
1846 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLc, 0 );
1849 TLc.print_to_file( ff1.str().c_str() );
1850 TLv.print_to_file( ffv.str().c_str() );
1856 buffer.buffer_init( sizeTuple * TLc.get_n() * 2 );
1859 TLc.print_to_file( ff1.str().c_str() );
1870 std::map< int, std::map< int, int > > availVertexIndicesPerProcessor;
1871 int nv = TLv.get_n();
1872 for(
int i = 0; i < nv; i++ )
1875 int orgProc = TLv.vi_rd[2 * i];
1876 int indexVert = TLv.vi_rd[2 * i + 1];
1877 availVertexIndicesPerProcessor[orgProc][indexVert] = i;
1887 int n = TLc.get_n();
1888 std::map< int, int > currentProcsCount;
1892 std::map< int, std::set< int > > sourcesForTasks;
1896 int currentSourceID = TLc.vi_rd[sizeTuple * 0 + 1];
1897 int proc0 = TLc.vi_rd[sizeTuple * 0];
1898 currentProcsCount[proc0] = 1;
1900 for(
int i = 1; i < n; i++ )
1902 int proc = TLc.vi_rd[sizeTuple * i];
1903 int sourceID = TLc.vi_rd[sizeTuple * i + 1];
1904 if( sourceID == currentSourceID )
1906 if( currentProcsCount.find( proc ) == currentProcsCount.end() )
1908 currentProcsCount[proc] = 1;
1911 currentProcsCount[proc]++;
1913 if( sourceID != currentSourceID || ( ( n - 1 ) == i ) )
1917 if( currentProcsCount.size() > 1 )
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();
1925 int procID = it->first;
1926 int numOverCells = it->second;
1927 std::cout <<
" task:" << procID <<
" " << numOverCells <<
" cells\n";
1932 for( std::map< int, int >::iterator it1 = currentProcsCount.begin(); it1 != currentProcsCount.end();
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++ )
1940 int proc2 = it2->first;
1941 if( proc1 != proc2 ) sizeOfTLc2 += it2->second;
1946 if( sourceID != currentSourceID )
1948 currentSourceID = sourceID;
1949 currentProcsCount.clear();
1950 currentProcsCount[proc] = 1;
1959 std::cout <<
" need to initialize TLc2 with " << sizeOfTLc2 <<
" cells\n ";
1963 int sizeTuple2 = 5 + globalMaxEdges;
1966 TLc2.initialize( sizeTuple2, 0, 0, 0, sizeOfTLc2 );
1967 TLc2.enableWriteAccess();
1970 std::map< int, std::set< int > > verticesToSendForProc;
1972 for(
int i = 0; i < n; i++ )
1974 int sourceID = TLc.vi_rd[sizeTuple * i + 1];
1975 if( sourcesForTasks.find( sourceID ) != sourcesForTasks.end() )
1978 std::set< int > procs = sourcesForTasks[sourceID];
1979 if( procs.size() < 2 )
MB_CHK_SET_ERR( MB_FAILURE,
" not enough processes involved with a sourceID cell" );
1981 int orgProc = TLc.vi_rd[sizeTuple * i];
1984 std::map< int, int >& availableVerticesFromThisProc = availVertexIndicesPerProcessor[orgProc];
1985 for( std::set< int >::iterator setIt = procs.begin(); setIt != procs.end(); setIt++ )
1987 int procID = *setIt;
1990 if( procID != orgProc )
1993 int n2 = TLc2.get_n();
1994 if( n2 >= sizeOfTLc2 )
MB_CHK_SET_ERR( MB_FAILURE,
" memory overflow" );
1996 std::set< int >& indexVerticesInTLv = verticesToSendForProc[procID];
1997 TLc2.vi_wr[n2 * sizeTuple2] = procID;
1998 TLc2.vi_wr[n2 * sizeTuple2 + 1] = orgProc;
1999 TLc2.vi_wr[n2 * sizeTuple2 + 2] = sourceID;
2000 TLc2.vi_wr[n2 * sizeTuple2 + 3] = TLc.vi_rd[sizeTuple * i + 2];
2002 int nvert = TLc.vi_rd[sizeTuple * i + 3];
2003 TLc2.vi_wr[n2 * sizeTuple2 + 4] = nvert;
2008 for(
int j = 0; j < nvert; j++ )
2010 int vertexIndex = TLc.vi_rd[i * sizeTuple + 4 + j];
2012 if( availableVerticesFromThisProc.find( vertexIndex ) == availableVerticesFromThisProc.end() )
2014 MB_CHK_SET_ERR( MB_FAILURE,
" vertex index not available from processor" );
2016 TLc2.vi_wr[n2 * sizeTuple2 + 5 + j] = vertexIndex;
2017 int indexInTLv = availVertexIndicesPerProcessor[orgProc][vertexIndex];
2018 indexVerticesInTLv.insert( indexInTLv );
2021 for(
int j = nvert; j < globalMaxEdges; j++ )
2023 TLc2.vi_wr[n2 * sizeTuple2 + 5 + j] = 0;
2037 for( std::map<
int, std::set< int > >::iterator it = verticesToSendForProc.begin();
2038 it != verticesToSendForProc.end(); it++ )
2040 std::set< int >& indexInTLvSet = it->second;
2041 numVerts2 += (int)indexInTLvSet.size();
2043 TLv2.initialize( 3, 0, 0, 3,
2045 TLv2.enableWriteAccess();
2046 for( std::map<
int, std::set< int > >::iterator it = verticesToSendForProc.begin();
2047 it != verticesToSendForProc.end(); it++ )
2049 int sendToProc = it->first;
2050 std::set< int >& indexInTLvSet = it->second;
2052 for( std::set< int >::iterator itSet = indexInTLvSet.begin(); itSet != indexInTLvSet.end(); itSet++ )
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];
2068 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv2, 0 );
2069 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLc2, 0 );
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() );
2087 int nvNew = TLv2.get_n();
2095 std::map< int, std::map< int, EntityHandle > > vertexPerProcAndIndex;
2096 for(
int i = 0; i < nvNew; i++ )
2098 int orgProc = TLv2.vi_rd[3 * i + 1];
2099 int indexInVert = TLv2.vi_rd[3 * i + 2];
2100 vertexPerProcAndIndex[orgProc][indexInVert] = newVerts[i];
2108 int ne = TLc2.get_n();
2109 for(
int i = 0; i < ne; i++ )
2111 int orgProc = TLc2.vi_rd[i * sizeTuple2 + 1];
2112 int sourceID = TLc2.vi_rd[i * sizeTuple2 + 2];
2113 int targetID = TLc2.vi_wr[i * sizeTuple2 + 3];
2114 int nve = TLc2.vi_wr[i * sizeTuple2 + 4];
2115 std::vector< EntityHandle > conn;
2117 for(
int j = 0; j < nve; j++ )
2119 int indexV = TLc2.vi_wr[i * sizeTuple2 + 5 + j];
2120 EntityHandle vh = vertexPerProcAndIndex[orgProc][indexV];
2125 newPolygons.insert( polyNew );
2137 std::stringstream ffs4;
2138 ffs4 <<
"extraIntxCells" <<
rank <<
".h5m";