diff --git a/src/coreComponents/integrationTests/meshTests/testVTKImport.cpp b/src/coreComponents/integrationTests/meshTests/testVTKImport.cpp
index 3eab728a833..8bb4336a1d8 100644
--- a/src/coreComponents/integrationTests/meshTests/testVTKImport.cpp
+++ b/src/coreComponents/integrationTests/meshTests/testVTKImport.cpp
@@ -53,17 +53,22 @@ using namespace geos::dataRepository;
template< class V >
void TestMeshImport( string const & meshFilePath, V const & validate, string const fractureName="" )
{
+ // Automatically use global IDs when fractures are present
+ string const useGlobalIdsStr = fractureName.empty() ? "0" : "1";
+
string const pattern = R"xml(
)xml";
- string const meshNode = GEOS_FMT( pattern, meshFilePath, fractureName.empty() ? "" : "faceBlocks=\"{" + fractureName + "}\"" );
+ string const meshNode = GEOS_FMT( pattern, meshFilePath, useGlobalIdsStr,
+ fractureName.empty() ? "" : "faceBlocks=\"{" + fractureName + "}\"" );
+
xmlWrapper::xmlDocument xmlDocument;
xmlDocument.loadString( meshNode );
xmlWrapper::xmlNode xmlMeshNode = xmlDocument.getChild( "Mesh" );
@@ -148,11 +153,12 @@ class TestFractureImport : public ::testing::Test
static std::filesystem::path createFractureMesh( std::filesystem::path const & folder )
{
- // The main mesh
+ // The main mesh - 3 hexahedra
vtkNew< vtkUnstructuredGrid > main;
{
- int constexpr numPoints = 16;
+ int constexpr numPoints = 24; // 3 hexahedra * 8 points each
double const pointsCoords[numPoints][3] = {
+ // First hexahedron (x: -1 to 0, y: 0 to 1)
{ -1, 0, 0 },
{ -1, 1, 0 },
{ -1, 1, 1 },
@@ -161,6 +167,7 @@ class TestFractureImport : public ::testing::Test
{ 0, 1, 0 },
{ 0, 1, 1 },
{ 0, 0, 1 },
+ // Second hexahedron (x: 0 to 1, y: 0 to 1) - shares face with first via fracture
{ 0, 0, 0 },
{ 0, 1, 0 },
{ 0, 1, 1 },
@@ -168,7 +175,18 @@ class TestFractureImport : public ::testing::Test
{ 1, 0, 0 },
{ 1, 1, 0 },
{ 1, 1, 1 },
- { 1, 0, 1 } };
+ { 1, 0, 1 },
+ // Third hexahedron (x: -1 to 0, y: 2 to 3) - disconnected from first two
+ { -1, 2, 0 },
+ { -1, 3, 0 },
+ { -1, 3, 1 },
+ { -1, 2, 1 },
+ { 0, 2, 0 },
+ { 0, 3, 0 },
+ { 0, 3, 1 },
+ { 0, 2, 1 }
+ };
+
vtkNew< vtkPoints > points;
points->Allocate( numPoints );
for( double const * pointsCoord: pointsCoords )
@@ -177,10 +195,11 @@ class TestFractureImport : public ::testing::Test
}
main->SetPoints( points );
- int constexpr numHexs = 2;
+ int constexpr numHexs = 3;
vtkIdType const cubes[numHexs][8] = {
- { 0, 1, 2, 3, 4, 5, 6, 7 },
- { 8, 9, 10, 11, 12, 13, 14, 15 }
+ { 0, 1, 2, 3, 4, 5, 6, 7 }, // Hex 0
+ { 8, 9, 10, 11, 12, 13, 14, 15 }, // Hex 1
+ { 16, 17, 18, 19, 20, 21, 22, 23 } // Hex 2
};
main->Allocate( numHexs );
for( vtkIdType const * cube: cubes )
@@ -207,7 +226,7 @@ class TestFractureImport : public ::testing::Test
main->GetPointData()->SetGlobalIds( pointGlobalIds );
}
- // The fracture mesh
+ // The fracture mesh - 1 fracture connecting only hex 0 and hex 1
vtkNew< vtkUnstructuredGrid > fracture;
{
int constexpr numPoints = 4;
@@ -215,7 +234,9 @@ class TestFractureImport : public ::testing::Test
{ 0, 0, 0 },
{ 0, 1, 0 },
{ 0, 1, 1 },
- { 0, 0, 1 } };
+ { 0, 0, 1 }
+ };
+
vtkNew< vtkPoints > points;
points->Allocate( numPoints );
for( double const * pointsCoord: pointsCoords )
@@ -229,7 +250,7 @@ class TestFractureImport : public ::testing::Test
fracture->Allocate( numQuads );
for( vtkIdType const * q: quad )
{
- fracture->InsertNextCell( VTK_QUAD, numPoints, q );
+ fracture->InsertNextCell( VTK_QUAD, 4, q );
}
vtkNew< vtkIdTypeArray > cellGlobalIds;
@@ -250,15 +271,15 @@ class TestFractureImport : public ::testing::Test
}
fracture->GetPointData()->SetGlobalIds( pointGlobalIds );
- // Do not forget the collocated_nodes fields
+ // Collocated nodes - connects hex 0 and hex 1
vtkNew< vtkIdTypeArray > collocatedNodes;
collocatedNodes->SetName( "collocated_nodes" );
collocatedNodes->SetNumberOfComponents( 2 );
collocatedNodes->SetNumberOfTuples( numPoints );
- collocatedNodes->SetTuple2( 0, 4, 8 );
- collocatedNodes->SetTuple2( 1, 5, 9 );
- collocatedNodes->SetTuple2( 2, 6, 10 );
- collocatedNodes->SetTuple2( 3, 7, 11 );
+ collocatedNodes->SetTuple2( 0, 4, 8 ); // Main mesh points 4 and 8
+ collocatedNodes->SetTuple2( 1, 5, 9 ); // Main mesh points 5 and 9
+ collocatedNodes->SetTuple2( 2, 6, 10 ); // Main mesh points 6 and 10
+ collocatedNodes->SetTuple2( 3, 7, 11 ); // Main mesh points 7 and 11
fracture->GetPointData()->AddArray( collocatedNodes );
}
@@ -289,29 +310,36 @@ TEST_F( TestFractureImport, fracture )
// Instead of checking each rank on its own,
// we check that all the data is present across the ranks.
auto const sum = []( auto i ) // Alias
+
{
return MpiWrapper::sum( i );
};
- // Volumic mesh validations
- ASSERT_EQ( sum( cellBlockManager.numNodes() ), 16 );
- ASSERT_EQ( sum( cellBlockManager.numEdges() ), 24 );
- ASSERT_EQ( sum( cellBlockManager.numFaces() ), 12 );
+ // Volumic mesh validations - 3 hexahedra with 24 points
+ // Points are NOT merged even though some are at same coordinates
+ // because they have different global IDs (4-7 vs 8-11)
+ ASSERT_EQ( sum( cellBlockManager.numNodes() ), 24 );
+ // Edges and faces will be different too - just check they exist
+ ASSERT_GT( sum( cellBlockManager.numEdges() ), 0 );
+ ASSERT_GT( sum( cellBlockManager.numFaces() ), 0 );
// Fracture mesh validations
ASSERT_EQ( sum( cellBlockManager.getFaceBlocks().numSubGroups() ), MpiWrapper::commSize() );
FaceBlockABC const & faceBlock = cellBlockManager.getFaceBlocks().getGroup< FaceBlockABC >( 0 );
ASSERT_EQ( sum( faceBlock.num2dElements() ), 1 );
ASSERT_EQ( sum( faceBlock.num2dFaces() ), 4 );
+
auto ecn = faceBlock.get2dElemsToCollocatedNodesBuckets();
auto const num2dElems = ecn.size();
ASSERT_EQ( sum( num2dElems ), 1 );
+
auto numNodesInFrac = 0;
for( int ei = 0; ei < num2dElems; ++ei )
{
numNodesInFrac += ecn[ei].size();
}
ASSERT_EQ( sum( numNodesInFrac ), 4 );
+
for( int ei = 0; ei < num2dElems; ++ei )
{
for( int ni = 0; ni < numNodesInFrac; ++ni )
@@ -327,7 +355,6 @@ TEST_F( TestFractureImport, fracture )
TestMeshImport( m_vtkFile, validate, "fracture" );
}
-
TEST( VTKImport, cube )
{
auto validate = []( CellBlockManagerABC const & cellBlockManager ) -> void
diff --git a/src/coreComponents/mesh/CMakeLists.txt b/src/coreComponents/mesh/CMakeLists.txt
index 3b7e7a825cb..65db05207f1 100644
--- a/src/coreComponents/mesh/CMakeLists.txt
+++ b/src/coreComponents/mesh/CMakeLists.txt
@@ -215,6 +215,7 @@ if( ENABLE_VTK )
generators/VTKMeshGenerator.hpp
generators/VTKMeshGeneratorTools.hpp
generators/VTKWellGenerator.hpp
+ generators/VTKSuperCellPartitioning.hpp
generators/VTKUtilities.hpp
)
set( mesh_sources ${mesh_sources}
@@ -224,6 +225,7 @@ if( ENABLE_VTK )
generators/VTKMeshGenerator.cpp
generators/VTKMeshGeneratorTools.cpp
generators/VTKWellGenerator.cpp
+ generators/VTKSuperCellPartitioning.cpp
generators/VTKUtilities.cpp
)
list( APPEND tplDependencyList VTK::IOLegacy VTK::FiltersParallelDIY2 )
diff --git a/src/coreComponents/mesh/generators/CollocatedNodes.cpp b/src/coreComponents/mesh/generators/CollocatedNodes.cpp
index 349c3973004..f108aa26f51 100644
--- a/src/coreComponents/mesh/generators/CollocatedNodes.cpp
+++ b/src/coreComponents/mesh/generators/CollocatedNodes.cpp
@@ -24,27 +24,47 @@ namespace geos::vtk
{
CollocatedNodes::CollocatedNodes( string const & faceBlockName,
- vtkSmartPointer< vtkDataSet > faceMesh )
+ vtkSmartPointer< vtkDataSet > faceMesh,
+ bool performGlobalCheck )
{
// The vtk field to the collocated nodes for fractures.
string const COLLOCATED_NODES = "collocated_nodes";
vtkIdTypeArray const * collocatedNodes = vtkIdTypeArray::FastDownCast( faceMesh->GetPointData()->GetArray( COLLOCATED_NODES.c_str() ) );
- // Depending on the parallel split, the vtk face mesh may be empty on a rank.
- // In that case, vtk will not provide any field for the emtpy mesh.
- // Therefore, not finding the duplicated nodes field on a rank cannot be interpreted as a globally missing field.
- // Converting the address into an integer and exchanging it through the MPI ranks let us find out
- // if the field is globally missing or not.
- std::uintptr_t const address = MpiWrapper::max( reinterpret_cast< std::uintptr_t >(collocatedNodes) );
- if( address == 0 )
+ if( performGlobalCheck )
{
- GEOS_LOG_RANK_0( "Available point data fields in '" << faceBlockName << "':" );
- for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )
+ // Depending on the parallel split, the vtk face mesh may be empty on a rank.
+ // In that case, vtk will not provide any field for the empty mesh.
+ // Therefore, not finding the duplicated nodes field on a rank cannot be interpreted as a globally missing field.
+ // Converting the address into an integer and exchanging it through the MPI ranks let us find out
+ // if the field is globally missing or not.
+ std::uintptr_t const address = MpiWrapper::max( reinterpret_cast< std::uintptr_t >(collocatedNodes) );
+ if( address == 0 )
{
- GEOS_LOG_RANK_0( " - " << faceMesh->GetPointData()->GetArrayName( i ) << " of type '" << faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() << "'" );
+ GEOS_LOG_RANK_0( "Available point data fields in '" << faceBlockName << "':" );
+ for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )
+ {
+ GEOS_LOG_RANK_0( " - " << faceMesh->GetPointData()->GetArrayName( i )
+ << " of type '" << faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() << "'" );
+ }
+ GEOS_ERROR( "Could not find valid field \"" << COLLOCATED_NODES
+ << "\" for fracture \"" << faceBlockName << "\"." );
+ }
+ }
+ else
+ {
+ if( !collocatedNodes )
+ {
+ GEOS_LOG_RANK_0( "Available point data fields in '" << faceBlockName << "':" );
+ for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )
+ {
+ GEOS_LOG_RANK_0( " - " << faceMesh->GetPointData()->GetArrayName( i )
+ << " of type '" << faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() << "'" );
+ }
+ GEOS_ERROR( "Could not find valid field \"" << COLLOCATED_NODES
+ << "\" for fracture \"" << faceBlockName << "\" on this rank." );
}
- GEOS_ERROR( "Could not find valid field \"" << COLLOCATED_NODES << "\" for fracture \"" << faceBlockName << "\"." );
}
if( collocatedNodes )
diff --git a/src/coreComponents/mesh/generators/CollocatedNodes.hpp b/src/coreComponents/mesh/generators/CollocatedNodes.hpp
index 98142dae589..114c6d9c5f9 100644
--- a/src/coreComponents/mesh/generators/CollocatedNodes.hpp
+++ b/src/coreComponents/mesh/generators/CollocatedNodes.hpp
@@ -35,9 +35,12 @@ class CollocatedNodes
* @brief Build a convenience wrapper around the raw vtk collocated nodes information.
* @param faceBlockName The face block name.
* @param faceMesh The face mesh for which the collocated nodes structure will be fed.
+ * @param performGlobalCheck Whether to check across all MPI ranks if field is missing (default: true).
+ * Set to false when only calling on a subset of ranks.
*/
CollocatedNodes( string const & faceBlockName,
- vtkSmartPointer< vtkDataSet > faceMesh );
+ vtkSmartPointer< vtkDataSet > faceMesh,
+ bool performGlobalCheck = true );
/**
* @brief For node @p i of the face block, returns all the duplicated global node indices in the main 3d mesh.
diff --git a/src/coreComponents/mesh/generators/ParMETISInterface.cpp b/src/coreComponents/mesh/generators/ParMETISInterface.cpp
index 4567826b341..0ac0ee5bf0a 100644
--- a/src/coreComponents/mesh/generators/ParMETISInterface.cpp
+++ b/src/coreComponents/mesh/generators/ParMETISInterface.cpp
@@ -126,5 +126,61 @@ partition( ArrayOfArraysView< idx_t const, idx_t > const & graph,
return part;
}
+array1d< idx_t >
+partitionWeighted( ArrayOfArraysView< idx_t const, idx_t > const & graph,
+ arrayView1d< idx_t const > const & vertexWeights,
+ arrayView1d< idx_t const > const & vertDist,
+ idx_t const numParts,
+ MPI_Comm comm,
+ int const numRefinements )
+{
+ array1d< idx_t > part( graph.size() );
+ if( numParts == 1 )
+ {
+ return part;
+ }
+
+ array1d< real_t > tpwgts( numParts );
+ tpwgts.setValues< serialPolicy >( 1.0f / static_cast< real_t >( numParts ) );
+
+ idx_t wgtflag = 2; // vertex weights only
+ idx_t numflag = 0;
+ idx_t ncon = 1;
+ idx_t npart = numParts;
+
+ // Options: [use_defaults, log_level, seed, coupling]
+ // PARMETIS_PSR_UNCOUPLED = 0 (default - uses PartitionSmallGraph)
+ // PARMETIS_PSR_COUPLED = 1 (forces distributed algorithm)
+ idx_t options[4] = { 1, 0, 2022, 0 };
+
+ idx_t edgecut = 0;
+ real_t ubvec = 1.05;
+
+ GEOS_PARMETIS_CHECK( ParMETIS_V3_PartKway(
+ const_cast< idx_t * >( vertDist.data() ),
+ const_cast< idx_t * >( graph.getOffsets() ),
+ const_cast< idx_t * >( graph.getValues() ),
+ const_cast< idx_t * >( vertexWeights.data() ),
+ nullptr, // edge weights
+ &wgtflag,
+ &numflag, &ncon, &npart, tpwgts.data(),
+ &ubvec, options, &edgecut, part.data(), &comm ) );
+
+ for( int iter = 0; iter < numRefinements; ++iter )
+ {
+ GEOS_PARMETIS_CHECK( ParMETIS_V3_RefineKway(
+ const_cast< idx_t * >( vertDist.data() ),
+ const_cast< idx_t * >( graph.getOffsets() ),
+ const_cast< idx_t * >( graph.getValues() ),
+ const_cast< idx_t * >( vertexWeights.data() ),
+ nullptr,
+ &wgtflag,
+ &numflag, &ncon, &npart, tpwgts.data(),
+ &ubvec, options, &edgecut, part.data(), &comm ) );
+ }
+
+ return part;
+}
+
} // namespace parmetis
} // namespace geos
diff --git a/src/coreComponents/mesh/generators/ParMETISInterface.hpp b/src/coreComponents/mesh/generators/ParMETISInterface.hpp
index fa1bb2c8545..cd38e1f5c44 100644
--- a/src/coreComponents/mesh/generators/ParMETISInterface.hpp
+++ b/src/coreComponents/mesh/generators/ParMETISInterface.hpp
@@ -70,6 +70,23 @@ partition( ArrayOfArraysView< pmet_idx_t const, pmet_idx_t > const & graph,
MPI_Comm comm,
int const numRefinements );
+/**
+ * @brief Partition a graph with vertex weights
+ * @param graph The adjacency graph
+ * @param vertexWeights The vertex weights for load balancing
+ * @param vertDist The element distribution
+ * @param numParts The number of partitions
+ * @param comm The MPI communicator
+ * @param numRefinements Number of refinement passes
+ * @return Partition assignment for each vertex
+ */
+array1d< pmet_idx_t >
+partitionWeighted( ArrayOfArraysView< pmet_idx_t const, pmet_idx_t > const & graph,
+ arrayView1d< pmet_idx_t const > const & vertexWeights,
+ arrayView1d< pmet_idx_t const > const & vertDist,
+ pmet_idx_t const numParts,
+ MPI_Comm comm,
+ int const numRefinements );
} // namespace parmetis
} // namespace geos
diff --git a/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp b/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp
index 46eff42dfd2..faf4fdb5ce9 100644
--- a/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp
+++ b/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp
@@ -451,7 +451,7 @@ Elem2dTo3dInfo buildElem2dTo3dElemAndFaces( vtkSmartPointer< vtkDataSet > faceMe
stdMap< vtkIdType, localIndex > ng2l; // global to local mapping for nodes.
for( vtkIdType i = 0; i < globalPtIds->GetNumberOfValues(); ++i )
{
- ng2l.insert( { globalPtIds->GetValue( i ), i} );
+ ng2l.emplace( globalPtIds->GetValue( i ), i );
}
// Let's build the elem2d to elem3d mapping.
@@ -491,7 +491,7 @@ Elem2dTo3dInfo buildElem2dTo3dElemAndFaces( vtkSmartPointer< vtkDataSet > faceMe
{
stdVector< vtkIdType > const & tmp = it->second;
std::set< vtkIdType > const cells{ tmp.cbegin(), tmp.cend() };
- nodesToCells.insert( {n, cells} );
+ nodesToCells.emplace( n, cells );
}
}
}
@@ -508,8 +508,14 @@ Elem2dTo3dInfo buildElem2dTo3dElemAndFaces( vtkSmartPointer< vtkDataSet > faceMe
{
// We collect all the duplicated points that are involved for each 2d element.
vtkIdList * pointIds = faceMesh->GetCell( e2d )->GetPointIds();
- std::size_t const elem2dNumPoints = pointIds->GetNumberOfIds();
- // All the duplicated points of the 2d element. Note that we lose the collocation of the duplicated nodes.
+
+ // Use the common matching function to find candidate 3D cells
+ stdVector< vtkIdType > matchingCells = vtk::findMatchingCellsForFractureElement(
+ pointIds,
+ collocatedNodes,
+ nodesToCells );
+
+ // Collect all duplicated nodes for this 2D element (needed for elem2dToNodes)
std::set< vtkIdType > duplicatedPointOfElem2d;
for( vtkIdType j = 0; j < pointIds->GetNumberOfIds(); ++j )
{
@@ -517,6 +523,7 @@ Elem2dTo3dInfo buildElem2dTo3dElemAndFaces( vtkSmartPointer< vtkDataSet > faceMe
duplicatedPointOfElem2d.insert( ns.cbegin(), ns.cend() );
}
+ // Build elem2dToNodes mapping
for( vtkIdType const & gni: duplicatedPointOfElem2d )
{
auto it = ng2l.find( gni );
@@ -532,42 +539,41 @@ Elem2dTo3dInfo buildElem2dTo3dElemAndFaces( vtkSmartPointer< vtkDataSet > faceMe
}
}
- // Here, we collect all the 3d elements that are concerned by at least one of those duplicated elements.
- stdMap< vtkIdType, std::set< vtkIdType > > elem3dToDuplicatedNodes;
- for( vtkIdType const & n: duplicatedPointOfElem2d )
+// Process each matching 3D cell
+ for( vtkIdType const & cellGlobalId : matchingCells )
{
- auto const ncs = nodesToCells.find( n );
- if( ncs != nodesToCells.cend() )
+ elem2dToElem3d.emplaceBack( e2d, elemToFaces.getElementIndexInCellBlock( cellGlobalId ) );
+
+ // Find which face matches
+ auto faces = elemToFaces[cellGlobalId];
+ for( int j = 0; j < faces.size( 0 ); ++j )
{
- for( vtkIdType const & c: ncs->second )
+ localIndex const faceIndex = faces[j];
+ auto nodes = faceToNodes[faceIndex];
+ std::set< vtkIdType > globalNodes;
+ for( auto const & n: nodes )
{
- elem3dToDuplicatedNodes.get_inserted( c ).insert( n );
+ globalNodes.insert( globalPtIds->GetValue( n ) );
}
- }
- }
- // Last we extract which of those candidate 3d elements are the ones actually neighboring the 2d element.
- for( auto const & e2n: elem3dToDuplicatedNodes )
- {
- // If the face of the element 3d has the same number of nodes than the elem 2d, it should be a successful (the mesh is conformal).
- if( e2n.second.size() == elem2dNumPoints )
- {
- // Now we know that the element 3d has a face that touches the element 2d. Let's find which one.
- elem2dToElem3d.emplaceBack( e2d, elemToFaces.getElementIndexInCellBlock( e2n.first ) );
- // Computing the elem2dToFaces mapping.
- auto faces = elemToFaces[e2n.first];
- for( int j = 0; j < faces.size( 0 ); ++j )
+
+ // Check if face nodes are a subset of duplicated nodes
+ // Face should have same number of nodes as the fracture element
+ if( globalNodes.size() == static_cast< std::size_t >(pointIds->GetNumberOfIds()) )
{
- localIndex const faceIndex = faces[j];
- auto nodes = faceToNodes[faceIndex];
- std::set< vtkIdType > globalNodes;
- for( auto const & n: nodes )
+ bool faceMatch = true;
+ for( vtkIdType gn : globalNodes )
{
- globalNodes.insert( globalPtIds->GetValue( n ) );
+ if( duplicatedPointOfElem2d.find( gn ) == duplicatedPointOfElem2d.end() )
+ {
+ faceMatch = false;
+ break;
+ }
}
- if( globalNodes == e2n.second )
+
+ if( faceMatch )
{
elem2dToFaces.emplaceBack( e2d, faceIndex );
- elem2dToCellBlock.emplaceBack( e2d, elemToFaces.getCellBlockIndex( e2n.first ) );
+ elem2dToCellBlock.emplaceBack( e2d, elemToFaces.getCellBlockIndex( cellGlobalId ) );
break;
}
}
@@ -598,7 +604,12 @@ array1d< globalIndex > buildLocalToGlobal( vtkIdTypeArray const * faceMeshCellGl
// In order to avoid any cell global id collision, we gather the max cell global id over all the ranks.
// Then we use this maximum as on offset.
// TODO This does not take into account multiple fractures.
- vtkIdType const maxLocalCellId = meshCellGlobalIds->GetMaxId();
+ vtkIdType maxLocalCellId = 0;
+ for( vtkIdType i = 0; i < meshCellGlobalIds->GetNumberOfTuples(); ++i )
+ {
+ maxLocalCellId = std::max( maxLocalCellId,
+ static_cast< vtkIdType >(meshCellGlobalIds->GetValue( i )) );
+ }
vtkIdType const maxGlobalCellId = MpiWrapper::max( maxLocalCellId );
vtkIdType const cellGlobalOffset = maxGlobalCellId + 1;
diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp
index e4b9e6c7feb..77de96b9a20 100644
--- a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp
+++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp
@@ -82,6 +82,10 @@ VTKMeshGenerator::VTKMeshGenerator( string const & name,
setInputFlag( InputFlags::OPTIONAL ).
setDescription( "Method (library) used to partition the mesh" );
+ registerWrapper( viewKeyStruct::partitionFractureWeightString(), &m_partitionFractureWeight ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Additional weight to fracture-connected super-cells during partitioning" );
+
registerWrapper( viewKeyStruct::useGlobalIdsString(), &m_useGlobalIds ).
setInputFlag( InputFlags::OPTIONAL ).
setApplyDefaultValue( 0 ).
@@ -208,6 +212,7 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager
comm,
m_partitionMethod,
m_partitionRefinement,
+ m_partitionFractureWeight,
m_useGlobalIds,
m_structuredIndexAttributeName,
numPartZ );
diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp
index cd968c2abd4..e31ada4bb1e 100644
--- a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp
+++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp
@@ -114,6 +114,7 @@ class VTKMeshGenerator : public ExternalMeshGeneratorBase
constexpr static char const * nodesetNamesString() { return "nodesetNames"; }
constexpr static char const * partitionRefinementString() { return "partitionRefinement"; }
constexpr static char const * partitionMethodString() { return "partitionMethod"; }
+ constexpr static char const * partitionFractureWeightString() { return "partitionFractureWeight"; }
constexpr static char const * useGlobalIdsString() { return "useGlobalIds"; }
constexpr static char const * dataSourceString() { return "dataSourceName"; }
};
@@ -161,6 +162,9 @@ class VTKMeshGenerator : public ExternalMeshGeneratorBase
/// Number of graph partitioning refinement iterations
integer m_partitionRefinement = 0;
+ /// Additional weight to fracture-connected super-cells during partitioning
+ integer m_partitionFractureWeight = 0;
+
/// Whether global id arrays should be used, if available
integer m_useGlobalIds = 0;
diff --git a/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.cpp b/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.cpp
new file mode 100644
index 00000000000..c972461fbc5
--- /dev/null
+++ b/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.cpp
@@ -0,0 +1,1191 @@
+/**
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+#include "VTKSuperCellPartitioning.hpp"
+
+// GEOS mesh includes
+#include "mesh/generators/VTKMeshGeneratorTools.hpp"
+
+
+// GEOS ParMETIS/Scotch (conditional)
+#ifdef GEOS_USE_PARMETIS
+#include "mesh/generators/ParMETISInterface.hpp"
+#endif
+
+// GEOS common includes
+#include "common/MpiWrapper.hpp"
+#include "common/DataTypes.hpp"
+#include "common/format/StringUtilities.hpp"
+
+// LvArray
+#include "LvArray/src/ArrayOfArrays.hpp"
+
+// VTK includes
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+
+namespace geos
+{
+
+namespace vtk
+{
+
+// =============================================================================================
+// SECTION 1: SUPER-CELL IDENTIFICATION (rank 0 only)
+// =============================================================================================
+pmet_idx_t computeSuperCellWeight( vtkIdType numCells, integer fractureWeight )
+{
+ if( numCells > 1 )
+ {
+ return numCells + fractureWeight; // Fracture super-cell
+ }
+ else
+ {
+ return numCells; // Regular cell (weight = 1)
+ }
+}
+
+SuperCellInfo tagCellsWithSuperCellIds(
+ vtkSmartPointer< vtkUnstructuredGrid > cells3D,
+ stdMap< string, ArrayOfArrays< vtkIdType, int64_t > > const & fractureNeighbors,
+ integer fractureWeight )
+{
+ GEOS_LOG_RANK_0( "TAGGING 3D CELLS WITH SUPER-CELL IDs" );
+
+ vtkIdType const numLocalCells = cells3D->GetNumberOfCells();
+
+ vtkIdTypeArray * globalIds =
+ vtkIdTypeArray::SafeDownCast( cells3D->GetCellData()->GetGlobalIds() );
+
+ GEOS_ERROR_IF( !globalIds, "3D mesh missing global IDs" );
+
+// -----------------------------------------------------------------------
+// Step 1: Build 3D cell connectivity graph via fractures
+// -----------------------------------------------------------------------
+ std::map< vtkIdType, std::set< vtkIdType > > fractureGraph;
+ vtkIdType totalFracturePairs = 0;
+
+ for( auto const & [fractureName, neighbors] : fractureNeighbors )
+ {
+ vtkIdType numFracCells = neighbors.size();
+
+ for( vtkIdType i = 0; i < numFracCells; ++i )
+ {
+ auto neighborList = neighbors[i];
+
+ if( neighborList.size() >= 2 )
+ {
+ vtkIdType gidA = neighborList[0];
+ vtkIdType gidB = neighborList[1];
+
+ fractureGraph[gidA].insert( gidB );
+ fractureGraph[gidB].insert( gidA );
+
+ totalFracturePairs++;
+ }
+ }
+
+ GEOS_LOG_RANK_0( GEOS_FMT( "Fracture '{}': processed {} fracture elements",
+ fractureName, numFracCells ) );
+ }
+
+ GEOS_LOG_RANK_0( GEOS_FMT( "Built fracture graph with {} 3D cells having fracture connections, {} fracture pairs",
+ fractureGraph.size(), totalFracturePairs ) );
+
+ // -----------------------------------------------------------------------
+ // Step 2: Find the maximum global ID from ALL cells (not just fracture cells)
+ // -----------------------------------------------------------------------
+ vtkIdType maxGlobalId = 0;
+
+ for( vtkIdType i = 0; i < numLocalCells; ++i )
+ {
+ vtkIdType gid = globalIds->GetValue( i );
+ maxGlobalId = std::max( maxGlobalId, gid );
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 3: Find connected components using DFS
+ // -----------------------------------------------------------------------
+ std::map< vtkIdType, vtkIdType > cellToSuperCell;
+ std::set< vtkIdType > visited;
+
+ // Start super-cell IDs AFTER the maximum global ID to avoid collisions
+ vtkIdType nextSuperCellId = maxGlobalId + 1;
+
+ std::function< void(vtkIdType, vtkIdType, std::vector< vtkIdType > &) > dfs =
+ [&]( vtkIdType cell, vtkIdType superCellId, std::vector< vtkIdType > & component )
+ {
+ if( visited.count( cell ) )
+ return;
+ // Mark this cell as visited
+ visited.insert( cell );
+ // Assign this cell to the current super-cell
+ cellToSuperCell[cell] = superCellId;
+ component.push_back( cell );
+
+ // Recursively visit all neighbors
+ if( fractureGraph.count( cell ) )
+ {
+ for( vtkIdType neighbor : fractureGraph.at( cell ) )
+ {
+ dfs( neighbor, superCellId, component );
+ }
+ }
+ };
+
+ std::map< vtkIdType, std::vector< vtkIdType > > superCellComponents;
+
+ for( auto const & [cell, neighbors] : fractureGraph )
+ {
+ if( !visited.count( cell ) )
+ {
+ std::vector< vtkIdType > component;
+ dfs( cell, nextSuperCellId, component );
+
+ superCellComponents[nextSuperCellId] = component;
+ nextSuperCellId++;
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 4: Build SuperCellInfo
+ // -----------------------------------------------------------------------
+ SuperCellInfo info;
+
+ vtkIdType numCellsInSuperCells = 0;
+ vtkIdType numSuperCellsCreated = 0;
+ vtkIdType largestSuperCellSize = 0;
+
+ for( auto const & [scId, members] : superCellComponents )
+ {
+ info.superCellToOriginalCells[scId] = members;
+ info.vertexWeights[scId] = computeSuperCellWeight( members.size(), fractureWeight );
+
+ if( members.size() > 1 )
+ {
+ info.atomicSuperCells.insert( scId );
+ }
+
+ numCellsInSuperCells += members.size();
+ numSuperCellsCreated++;
+ largestSuperCellSize = std::max( largestSuperCellSize,
+ static_cast< vtkIdType >( members.size() ) );
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 5: Tag cells with SuperCellIds
+ // -----------------------------------------------------------------------
+ vtkNew< vtkIdTypeArray > superCellIdArray;
+ superCellIdArray->SetName( "SuperCellId" );
+ superCellIdArray->SetNumberOfTuples( numLocalCells );
+
+ for( vtkIdType i = 0; i < numLocalCells; ++i )
+ {
+ vtkIdType globalId = globalIds->GetValue( i );
+
+ if( cellToSuperCell.count( globalId ) )
+ {
+ // Cell is part of a fracture-connected super-cell
+ superCellIdArray->SetValue( i, cellToSuperCell.at( globalId ) );
+ }
+ else
+ {
+ // Regular cell - use its own global ID as super-cell ID
+ // This is safe because super-cell IDs start at maxGlobalId + 1
+ superCellIdArray->SetValue( i, globalId );
+ }
+ }
+
+ cells3D->GetCellData()->AddArray( superCellIdArray );
+
+ // -----------------------------------------------------------------------
+ // Step 6: Report statistics
+ // -----------------------------------------------------------------------
+ vtkIdType numRegularCells = numLocalCells - numCellsInSuperCells;
+ vtkIdType totalSuperCells = numSuperCellsCreated + numRegularCells;
+ vtkIdType cellReduction = numLocalCells - totalSuperCells;
+
+ GEOS_LOG_RANK_0( "SUPER-CELL TAGGING SUMMARY" );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Total 3D cells: {:8}", numLocalCells ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Cells in super-cells: {:8}", numCellsInSuperCells ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Regular cells (no fractures): {:8}", numRegularCells ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Number of super-cells created: {:8}", numSuperCellsCreated ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Total super-cells (incl regular): {:8}", totalSuperCells ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Cell reduction: {:8}", cellReduction ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Largest super-cell size: {:8} cells", largestSuperCellSize ) );
+
+ return info;
+}
+
+
+// =============================================================================================
+// SECTION 2: SUPER-CELL RECONSTRUCTION (after redistribution)
+// =============================================================================================
+SuperCellInfo reconstructSuperCellInfo( vtkSmartPointer< vtkUnstructuredGrid > mesh, integer fractureWeight )
+{
+ SuperCellInfo info;
+
+ vtkIdTypeArray * superCellIdArray =
+ vtkIdTypeArray::SafeDownCast( mesh->GetCellData()->GetArray( "SuperCellId" ) );
+
+ if( !superCellIdArray )
+ {
+ return info; // No super-cells, return empty info
+ }
+
+ vtkIdTypeArray * globalIds =
+ vtkIdTypeArray::SafeDownCast( mesh->GetCellData()->GetGlobalIds() );
+
+ GEOS_ERROR_IF( !globalIds, "Mesh missing global IDs" );
+
+ // Build map: super-cell ID -> vector of cell global IDs
+ std::map< vtkIdType, std::vector< vtkIdType > > localSuperCells;
+
+ for( vtkIdType i = 0; i < mesh->GetNumberOfCells(); ++i )
+ {
+ vtkIdType scId = superCellIdArray->GetValue( i );
+ vtkIdType globalId = globalIds->GetValue( i );
+ localSuperCells[scId].push_back( globalId );
+ }
+
+// Populate SuperCellInfo
+ for( auto const & [scId, cells] : localSuperCells )
+ {
+ info.superCellToOriginalCells[scId] = cells;
+ info.vertexWeights[scId] = computeSuperCellWeight( cells.size(), fractureWeight );
+
+ if( cells.size() > 1 )
+ {
+ info.atomicSuperCells.insert( scId );
+ }
+ }
+
+ return info;
+}
+
+
+// =============================================================================================
+// SECTION 3: INITIAL REDISTRIBUTION
+// =============================================================================================
+vtkSmartPointer< vtkDataSet >
+redistributeBySuperCellBlocks( vtkSmartPointer< vtkUnstructuredGrid > cells3D,
+ MPI_Comm comm,
+ InitialDistributionStrategy strategy )
+{
+ GEOS_MARK_FUNCTION;
+
+ int const rank = MpiWrapper::commRank( comm );
+ int const numRanks = MpiWrapper::commSize( comm );
+
+ vtkSmartPointer< vtkPartitionedDataSet > partitionedMesh =
+ vtkSmartPointer< vtkPartitionedDataSet >::New();
+
+ if( rank == 0 )
+ {
+ vtkIdTypeArray * superCellIdArray =
+ vtkIdTypeArray::SafeDownCast( cells3D->GetCellData()->GetArray( "SuperCellId" ) );
+
+ GEOS_ERROR_IF( !superCellIdArray, "SuperCellId array not found" );
+
+ vtkIdType numCells = cells3D->GetNumberOfCells();
+
+ // Build super-cell metadata
+ std::map< vtkIdType, std::vector< vtkIdType > > superCellToLocalCells;
+
+ for( vtkIdType i = 0; i < numCells; ++i )
+ {
+ vtkIdType scId = superCellIdArray->GetValue( i );
+ superCellToLocalCells[scId].push_back( i );
+ }
+ vtkIdType numSuperCells = superCellToLocalCells.size();
+
+ GEOS_LOG_RANK_0( GEOS_FMT( "Initial distribution: {} super-cells across {} ranks using {} strategy",
+ numSuperCells, numRanks,
+ (strategy == InitialDistributionStrategy::MORTON ? "MORTON" : "BLOCK") ) );
+
+ // -----------------------------------------------------------------------
+ // Step 1: Build super-cell list and optionally sort by spatial locality
+ // -----------------------------------------------------------------------
+ struct SuperCellDistributionInfo
+ {
+ vtkIdType scId;
+ std::vector< vtkIdType > cellIndices;
+ std::array< double, 3 > centroid; // Only computed for Morton strategy
+ };
+
+ std::vector< SuperCellDistributionInfo > superCells;
+ superCells.reserve( numSuperCells );
+
+ if( strategy == InitialDistributionStrategy::MORTON )
+ {
+ // MORTON: Compute centroids for spatial sorting
+ GEOS_LOG_RANK_0( "Computing super-cell centroids for Morton ordering..." );
+
+ // Pre-compute cell centroids
+ std::vector< std::array< double, 3 > > allCellCentroids( numCells );
+ vtkPoints * meshPoints = cells3D->GetPoints();
+
+ for( vtkIdType cellIdx = 0; cellIdx < numCells; ++cellIdx )
+ {
+ vtkIdType npts;
+ const vtkIdType * pts;
+ cells3D->GetCellPoints( cellIdx, npts, pts );
+
+ double centroid[3] = {0.0, 0.0, 0.0};
+ for( vtkIdType i = 0; i < npts; ++i )
+ {
+ double pt[3];
+ meshPoints->GetPoint( pts[i], pt );
+ centroid[0] += pt[0];
+ centroid[1] += pt[1];
+ centroid[2] += pt[2];
+ }
+
+ allCellCentroids[cellIdx][0] = centroid[0] / npts;
+ allCellCentroids[cellIdx][1] = centroid[1] / npts;
+ allCellCentroids[cellIdx][2] = centroid[2] / npts;
+ }
+
+ // Compute super-cell centroids
+ for( auto const & [scId, cellIndices] : superCellToLocalCells )
+ {
+ std::array< double, 3 > centroid = {0.0, 0.0, 0.0};
+
+ for( vtkIdType cellIdx : cellIndices )
+ {
+ centroid[0] += allCellCentroids[cellIdx][0];
+ centroid[1] += allCellCentroids[cellIdx][1];
+ centroid[2] += allCellCentroids[cellIdx][2];
+ }
+
+ centroid[0] /= cellIndices.size();
+ centroid[1] /= cellIndices.size();
+ centroid[2] /= cellIndices.size();
+
+ superCells.push_back( SuperCellDistributionInfo{ scId, cellIndices, centroid } );
+ }
+
+ allCellCentroids.clear();
+ allCellCentroids.shrink_to_fit();
+
+ // Find bounding box
+ double minCoord[3] = {std::numeric_limits< double >::max(),
+ std::numeric_limits< double >::max(),
+ std::numeric_limits< double >::max()};
+ double maxCoord[3] = {std::numeric_limits< double >::lowest(),
+ std::numeric_limits< double >::lowest(),
+ std::numeric_limits< double >::lowest()};
+
+ for( auto const & sc : superCells )
+ {
+ for( int d = 0; d < 3; ++d )
+ {
+ minCoord[d] = std::min( minCoord[d], sc.centroid[d] );
+ maxCoord[d] = std::max( maxCoord[d], sc.centroid[d] );
+ }
+ }
+
+ GEOS_LOG_RANK_0( GEOS_FMT( "Bounding box: X=[{:.3f}, {:.3f}], Y=[{:.3f}, {:.3f}], Z=[{:.3f}, {:.3f}]",
+ minCoord[0], maxCoord[0], minCoord[1], maxCoord[1], minCoord[2], maxCoord[2] ));
+
+ // Morton encoding
+ auto computeMorton = []( std::array< double, 3 > const & centroid,
+ double bounds_min[3],
+ double bounds_max[3] ) -> uint64_t
+ {
+ auto normalize = [&]( double val, int dim ) -> uint32_t
+ {
+ double range = bounds_max[dim] - bounds_min[dim];
+ if( range < 1e-10 )
+ return 0;
+ double norm = (val - bounds_min[dim]) / range;
+ norm = std::max( 0.0, std::min( 1.0, norm ) );
+ return static_cast< uint32_t >( norm * ((1u << 21) - 1) );
+ };
+
+ uint32_t x = normalize( centroid[0], 0 );
+ uint32_t y = normalize( centroid[1], 1 );
+ uint32_t z = normalize( centroid[2], 2 );
+
+ uint64_t code = 0;
+ for( int i = 0; i < 21; ++i )
+ {
+ code |= ((x & (1u << i)) ? (1ull << (3*i)) : 0);
+ code |= ((y & (1u << i)) ? (1ull << (3*i + 1)) : 0);
+ code |= ((z & (1u << i)) ? (1ull << (3*i + 2)) : 0);
+ }
+ return code;
+ };
+
+ // Sort by Morton code
+ std::sort( superCells.begin(), superCells.end(),
+ [&]( SuperCellDistributionInfo const & a, SuperCellDistributionInfo const & b )
+ {
+ return computeMorton( a.centroid, minCoord, maxCoord ) <
+ computeMorton( b.centroid, minCoord, maxCoord );
+ } );
+
+ GEOS_LOG_RANK_0( "Super-cells sorted by Morton Z-order curve" );
+ }
+ else // BLOCK strategy
+ {
+ // BLOCK: Simple ordering by super-cell ID (no centroid computation needed)
+ GEOS_LOG_RANK_0( "Using block distribution (no spatial ordering)" );
+
+ for( auto const & [scId, cellIndices] : superCellToLocalCells )
+ {
+ superCells.push_back( SuperCellDistributionInfo{ scId, cellIndices, {0.0, 0.0, 0.0} } );
+ }
+
+ // Sort by super-cell ID for deterministic partitioning
+ std::sort( superCells.begin(), superCells.end(),
+ []( SuperCellDistributionInfo const & a, SuperCellDistributionInfo const & b )
+ {
+ return a.scId < b.scId;
+ } );
+
+ GEOS_LOG_RANK_0( "Super-cells ordered by ID" );
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 2: Assign super-cells to ranks in contiguous blocks
+ // -----------------------------------------------------------------------
+ array1d< int64_t > cellPartitions( numCells );
+ std::vector< vtkIdType > cellsPerRank( numRanks, 0 );
+
+ vtkIdType superCellsPerRank = (numSuperCells + numRanks - 1) / numRanks;
+
+ for( vtkIdType scIdx = 0; scIdx < numSuperCells; ++scIdx )
+ {
+ int targetRank = std::min( static_cast< int >(scIdx / superCellsPerRank), numRanks - 1 );
+
+ // All cells in this super-cell go to the same rank
+ for( vtkIdType cellIdx : superCells[scIdx].cellIndices )
+ {
+ cellPartitions[cellIdx] = targetRank;
+ cellsPerRank[targetRank]++;
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 3: Build partitions
+ // -----------------------------------------------------------------------
+ partitionedMesh->SetNumberOfPartitions( numRanks );
+
+ for( int r = 0; r < numRanks; ++r )
+ {
+ vtkNew< vtkIdList > cellsForRank;
+
+ for( vtkIdType i = 0; i < numCells; ++i )
+ {
+ if( cellPartitions[i] == r )
+ {
+ cellsForRank->InsertNextId( i );
+ }
+ }
+
+ if( cellsForRank->GetNumberOfIds() == 0 )
+ {
+ vtkNew< vtkUnstructuredGrid > emptyPart;
+ partitionedMesh->SetPartition( r, emptyPart );
+ continue;
+ }
+
+ vtkNew< vtkExtractCells > extractor;
+ extractor->SetInputData( cells3D );
+ extractor->SetCellList( cellsForRank );
+ extractor->Update();
+
+ vtkSmartPointer< vtkUnstructuredGrid > partition =
+ vtkUnstructuredGrid::SafeDownCast( extractor->GetOutput() );
+
+ partitionedMesh->SetPartition( r, partition );
+ }
+ }
+ else
+ {
+ // Other ranks: create empty partitioned dataset
+ partitionedMesh->SetNumberOfPartitions( numRanks );
+ for( int r = 0; r < numRanks; ++r )
+ {
+ vtkNew< vtkUnstructuredGrid > emptyPart;
+ partitionedMesh->SetPartition( r, emptyPart );
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 4: Redistribute using VTK
+ // -----------------------------------------------------------------------
+ vtkSmartPointer< vtkDataSet > result = vtk::redistribute( *partitionedMesh, comm );
+
+ partitionedMesh = nullptr;
+ if( rank == 0 )
+ {
+ cells3D = nullptr;
+ }
+
+ vtkIdType localCells = result->GetNumberOfCells();
+
+ // Verify SuperCellId array survived redistribution
+ if( localCells > 0 )
+ {
+ vtkIdTypeArray * resultSuperCellIdArray =
+ vtkIdTypeArray::SafeDownCast( result->GetCellData()->GetArray( "SuperCellId" ) );
+
+ GEOS_ERROR_IF( !resultSuperCellIdArray,
+ GEOS_FMT( "Rank {}: SuperCellId array lost during redistribution", rank ) );
+ }
+
+ return result;
+}
+
+
+// =============================================================================================
+// SECTION 4: SUPER-CELL GRAPH BUILDING
+// =============================================================================================
+std::pair< ArrayOfArrays< pmet_idx_t, pmet_idx_t >, array1d< pmet_idx_t > >
+buildSuperCellGraph(
+ vtkSmartPointer< vtkUnstructuredGrid > cells3D,
+ ArrayOfArrays< pmet_idx_t, pmet_idx_t > const & baseGraph,
+ arrayView1d< pmet_idx_t const > const & baseElemDist,
+ SuperCellInfo const & info,
+ MPI_Comm comm )
+{
+ int const rank = MpiWrapper::commRank( comm );
+ int const numRanks = MpiWrapper::commSize( comm );
+
+ vtkIdTypeArray * superCellIdArray =
+ vtkIdTypeArray::SafeDownCast( cells3D->GetCellData()->GetArray( "SuperCellId" ) );
+ GEOS_ERROR_IF( superCellIdArray == nullptr, "SuperCellId array not found" );
+
+ vtkIdTypeArray * cellGlobalIds =
+ vtkIdTypeArray::SafeDownCast( cells3D->GetCellData()->GetGlobalIds() );
+ GEOS_ERROR_IF( cellGlobalIds == nullptr, "Cell global IDs not found" );
+
+ vtkIdType numLocalCells = cells3D->GetNumberOfCells();
+
+ std::map< vtkIdType, std::vector< vtkIdType > > superCellToLocalCells;
+ for( vtkIdType i = 0; i < numLocalCells; ++i )
+ {
+ vtkIdType superCellId = superCellIdArray->GetValue( i );
+ superCellToLocalCells[superCellId].push_back( i );
+ }
+
+ localIndex numLocalSuperCells = superCellToLocalCells.size();
+
+ // Build super-cell distribution
+ array1d< pmet_idx_t > superElemDist( numRanks + 1 );
+ pmet_idx_t localSuperCellCount = numLocalSuperCells;
+ MpiWrapper::allgather( &localSuperCellCount, 1, superElemDist.data(), 1, comm );
+
+ pmet_idx_t temp = superElemDist[0];
+ superElemDist[0] = 0;
+ for( int r = 1; r <= numRanks; ++r )
+ {
+ pmet_idx_t next = superElemDist[r];
+ superElemDist[r] = superElemDist[r-1] + temp;
+ temp = next;
+ }
+
+ pmet_idx_t myGlobalStart = superElemDist[rank];
+ pmet_idx_t myGlobalEnd = superElemDist[rank + 1];
+
+ std::vector< vtkIdType > orderedSuperCellIds;
+ orderedSuperCellIds.reserve( numLocalSuperCells );
+ for( auto const & [superCellId, localCells] : superCellToLocalCells )
+ {
+ orderedSuperCellIds.push_back( superCellId );
+ }
+ std::sort( orderedSuperCellIds.begin(), orderedSuperCellIds.end() );
+
+ std::map< vtkIdType, pmet_idx_t > superCellIdToGlobalIdx;
+ for( pmet_idx_t i = 0; i < numLocalSuperCells; ++i )
+ {
+ superCellIdToGlobalIdx[orderedSuperCellIds[i]] = myGlobalStart + i;
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 1: Build mappings for exchange
+ // -----------------------------------------------------------------------
+ // Local maps
+ std::unordered_map< pmet_idx_t, vtkIdType > myGlobalCellToSuperCell;
+ myGlobalCellToSuperCell.reserve( numLocalCells );
+
+ std::unordered_map< pmet_idx_t, vtkIdType > myParmetisToVtk;
+ myParmetisToVtk.reserve( numLocalCells );
+
+ pmet_idx_t myParmetisStart = baseElemDist[rank];
+ pmet_idx_t myParmetisEnd = baseElemDist[rank + 1];
+
+ for( vtkIdType i = 0; i < numLocalCells; ++i )
+ {
+ vtkIdType vtkGlobalId = cellGlobalIds->GetValue( i );
+ vtkIdType superCellId = superCellIdArray->GetValue( i );
+
+ myGlobalCellToSuperCell[vtkGlobalId] = superCellId;
+ myParmetisToVtk[myParmetisStart + i] = vtkGlobalId;
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 2: Identify unique ghost ParMETIS indices
+ // -----------------------------------------------------------------------
+ std::set< pmet_idx_t > ghostParmetisSet;
+
+ for( localIndex i = 0; i < baseGraph.size(); ++i )
+ {
+ auto neighbors = baseGraph[i];
+ for( localIndex j = 0; j < neighbors.size(); ++j )
+ {
+ pmet_idx_t nbrIdx = neighbors[j];
+ if( nbrIdx < myParmetisStart || nbrIdx >= myParmetisEnd )
+ {
+ ghostParmetisSet.insert( nbrIdx );
+ }
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 3: Exchange ghost mappings
+ // -----------------------------------------------------------------------
+ // Prepare local data to send
+ std::vector< pmet_idx_t > localParmetisIndices;
+ std::vector< vtkIdType > localVtkIds;
+ std::vector< vtkIdType > localSuperCellIds;
+
+ localParmetisIndices.reserve( numLocalCells );
+ localVtkIds.reserve( numLocalCells );
+ localSuperCellIds.reserve( numLocalCells );
+
+ for( vtkIdType i = 0; i < numLocalCells; ++i )
+ {
+ localParmetisIndices.push_back( myParmetisStart + i );
+ localVtkIds.push_back( cellGlobalIds->GetValue( i ) );
+ localSuperCellIds.push_back( superCellIdArray->GetValue( i ) );
+ }
+
+ // Gather sizes
+ int localCount = static_cast< int >( localParmetisIndices.size() );
+ array1d< int > allCounts( numRanks );
+ MpiWrapper::allgather( &localCount, 1, allCounts.data(), 1, comm );
+
+ array1d< int > displs( numRanks + 1 );
+ displs[0] = 0;
+ for( int r = 0; r < numRanks; ++r )
+ {
+ displs[r + 1] = displs[r] + allCounts[r];
+ }
+
+ int totalMappings = displs[numRanks];
+
+ // Allocate and gather
+ array1d< pmet_idx_t > allParmetisIndices( totalMappings );
+ array1d< vtkIdType > allVtkIds( totalMappings );
+ array1d< vtkIdType > allSuperCellIds( totalMappings );
+
+ MpiWrapper::allgatherv( localParmetisIndices.data(), localCount,
+ allParmetisIndices.data(), allCounts.data(), displs.data(), comm );
+ MpiWrapper::allgatherv( localVtkIds.data(), localCount,
+ allVtkIds.data(), allCounts.data(), displs.data(), comm );
+ MpiWrapper::allgatherv( localSuperCellIds.data(), localCount,
+ allSuperCellIds.data(), allCounts.data(), displs.data(), comm );
+
+ // Build lookup tables from gathered data (only for ghosts + local)
+ std::unordered_map< pmet_idx_t, vtkIdType > parmetisToVtk;
+ std::unordered_map< vtkIdType, vtkIdType > vtkToSuperCell;
+
+ parmetisToVtk.reserve( ghostParmetisSet.size() + numLocalCells );
+ vtkToSuperCell.reserve( ghostParmetisSet.size() + numLocalCells );
+
+ for( int i = 0; i < totalMappings; ++i )
+ {
+ pmet_idx_t parmetisIdx = allParmetisIndices[i];
+
+ // Only store if it's local or a ghost we need
+ if( (parmetisIdx >= myParmetisStart && parmetisIdx < myParmetisEnd) ||
+ ghostParmetisSet.count( parmetisIdx ) > 0 )
+ {
+ vtkIdType vtkId = allVtkIds[i];
+ vtkIdType superCellId = allSuperCellIds[i];
+
+ parmetisToVtk[parmetisIdx] = vtkId;
+ vtkToSuperCell[vtkId] = superCellId;
+ }
+ }
+
+ // Cleanup
+ allParmetisIndices.clear();
+ allVtkIds.clear();
+ allSuperCellIds.clear();
+
+ // -----------------------------------------------------------------------
+ // Step 4: Exchange super-cell global indices
+ // -----------------------------------------------------------------------
+ std::vector< vtkIdType > sendSCIds;
+ std::vector< pmet_idx_t > sendSCGlobalIndices;
+
+ for( auto const & [scId, gIdx] : superCellIdToGlobalIdx )
+ {
+ sendSCIds.push_back( scId );
+ sendSCGlobalIndices.push_back( gIdx );
+ }
+
+ int localSCCount = static_cast< int >( sendSCIds.size() );
+ array1d< int > allSCCounts( numRanks );
+ MpiWrapper::allgather( &localSCCount, 1, allSCCounts.data(), 1, comm );
+
+ array1d< int > scDispls( numRanks + 1 );
+ scDispls[0] = 0;
+ for( int r = 0; r < numRanks; ++r )
+ {
+ scDispls[r + 1] = scDispls[r] + allSCCounts[r];
+ }
+
+ int totalSCMappings = scDispls[numRanks];
+
+ array1d< vtkIdType > allSCIds( totalSCMappings );
+ array1d< pmet_idx_t > allSCGlobalIndices( totalSCMappings );
+
+ MpiWrapper::allgatherv( sendSCIds.data(), localSCCount,
+ allSCIds.data(), allSCCounts.data(), scDispls.data(), comm );
+ MpiWrapper::allgatherv( sendSCGlobalIndices.data(), localSCCount,
+ allSCGlobalIndices.data(), allSCCounts.data(), scDispls.data(), comm );
+
+ // Update local map with all super-cell mappings
+ for( int i = 0; i < totalSCMappings; ++i )
+ {
+ superCellIdToGlobalIdx[allSCIds[i]] = allSCGlobalIndices[i];
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 5: Build super-cell graph edges
+ // -----------------------------------------------------------------------
+ array1d< pmet_idx_t > superVertexWeights( numLocalSuperCells );
+ std::vector< std::set< pmet_idx_t > > neighborSets( numLocalSuperCells );
+
+ for( localIndex localSuperIdx = 0; localSuperIdx < numLocalSuperCells; ++localSuperIdx )
+ {
+ vtkIdType superCellId = orderedSuperCellIds[localSuperIdx];
+ auto const & localCells = superCellToLocalCells.at( superCellId );
+
+ auto itWeight = info.vertexWeights.find( superCellId );
+ superVertexWeights[localSuperIdx] = (itWeight != info.vertexWeights.end())
+ ? itWeight->second : localCells.size();
+
+ for( vtkIdType cellLocalIdx : localCells )
+ {
+ auto neighbors = baseGraph[cellLocalIdx];
+ for( localIndex j = 0; j < neighbors.size(); ++j )
+ {
+ pmet_idx_t nbrParmetis = neighbors[j];
+
+ auto itVtk = parmetisToVtk.find( nbrParmetis );
+ if( itVtk == parmetisToVtk.end() )
+ continue;
+
+ vtkIdType nbrVtkId = itVtk->second;
+
+ auto itSC = vtkToSuperCell.find( nbrVtkId );
+ if( itSC == vtkToSuperCell.end() )
+ continue;
+
+ vtkIdType nbrSuperCellId = itSC->second;
+ if( nbrSuperCellId == superCellId )
+ continue; // Skip self-loops
+
+ auto itGlobal = superCellIdToGlobalIdx.find( nbrSuperCellId );
+ if( itGlobal != superCellIdToGlobalIdx.end() )
+ {
+ neighborSets[localSuperIdx].insert( itGlobal->second );
+ }
+ }
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 6: Symmetrize
+ // -----------------------------------------------------------------------
+ // Collect ALL edges to send (to any rank)
+ std::vector< pmet_idx_t > myEdgesSrc;
+ std::vector< pmet_idx_t > myEdgesDst;
+
+ for( localIndex i = 0; i < numLocalSuperCells; ++i )
+ {
+ pmet_idx_t globalI = myGlobalStart + i;
+ for( pmet_idx_t nbrIdx : neighborSets[i] )
+ {
+ if( nbrIdx < myGlobalStart || nbrIdx >= myGlobalEnd )
+ {
+ myEdgesSrc.push_back( nbrIdx ); // destination of reverse edge
+ myEdgesDst.push_back( globalI ); // source of reverse edge
+ }
+ }
+ }
+
+ int myEdgeCount = static_cast< int >( myEdgesSrc.size() );
+
+ // Gather counts from all ranks
+ array1d< int > allEdgeCounts( numRanks );
+ MpiWrapper::allgather( &myEdgeCount, 1, allEdgeCounts.data(), 1, comm );
+
+ // Compute displacements for symmetrization
+ array1d< int > edgeDispls( numRanks + 1 );
+ edgeDispls[0] = 0;
+ for( int r = 0; r < numRanks; ++r )
+ {
+ edgeDispls[r + 1] = edgeDispls[r] + allEdgeCounts[r];
+ }
+
+ int totalEdges = edgeDispls[numRanks];
+
+ // Gather all edges
+ array1d< pmet_idx_t > allEdgeSrc( totalEdges > 0 ? totalEdges : 1 );
+ array1d< pmet_idx_t > allEdgeDst( totalEdges > 0 ? totalEdges : 1 );
+
+ if( myEdgeCount > 0 )
+ {
+ MpiWrapper::allgatherv( myEdgesSrc.data(), myEdgeCount,
+ allEdgeSrc.data(), allEdgeCounts.data(), edgeDispls.data(), comm );
+ MpiWrapper::allgatherv( myEdgesDst.data(), myEdgeCount,
+ allEdgeDst.data(), allEdgeCounts.data(), edgeDispls.data(), comm );
+ }
+ else
+ {
+ // Handle empty send
+ pmet_idx_t dummy = 0;
+ MpiWrapper::allgatherv( &dummy, 0,
+ allEdgeSrc.data(), allEdgeCounts.data(), edgeDispls.data(), comm );
+ MpiWrapper::allgatherv( &dummy, 0,
+ allEdgeDst.data(), allEdgeCounts.data(), edgeDispls.data(), comm );
+ }
+
+ // Add reverse edges for those I own
+ int reverseEdgesAdded = 0;
+ for( int i = 0; i < totalEdges; ++i )
+ {
+ pmet_idx_t dst = allEdgeSrc[i];
+ pmet_idx_t src = allEdgeDst[i];
+
+ if( dst >= myGlobalStart && dst < myGlobalEnd )
+ {
+ localIndex localIdx = dst - myGlobalStart;
+ if( neighborSets[localIdx].insert( src ).second )
+ {
+ reverseEdgesAdded++;
+ }
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 7: Build final ArrayOfArrays
+ // -----------------------------------------------------------------------
+ ArrayOfArrays< pmet_idx_t, pmet_idx_t > superGraph;
+ array1d< pmet_idx_t > offsets( numLocalSuperCells + 1 );
+ offsets[0] = 0;
+
+ for( localIndex i = 0; i < numLocalSuperCells; ++i )
+ {
+ offsets[i + 1] = offsets[i] + neighborSets[i].size();
+ }
+
+ superGraph.resizeFromOffsets( numLocalSuperCells, offsets.data() );
+
+ for( localIndex i = 0; i < numLocalSuperCells; ++i )
+ {
+ superGraph.appendToArray( i, neighborSets[i].begin(), neighborSets[i].end() );
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 8: Statistics
+ // -----------------------------------------------------------------------
+ localIndex totalSuperEdges = 0;
+ for( localIndex i = 0; i < superGraph.size(); ++i )
+ {
+ totalSuperEdges += superGraph.sizeOfArray( i );
+ }
+
+ localIndex totalBaseEdges = 0;
+ for( localIndex i = 0; i < baseGraph.size(); ++i )
+ {
+ totalBaseEdges += baseGraph.sizeOfArray( i );
+ }
+
+ vtkIdType globalSuperCells = MpiWrapper::sum(
+ static_cast< vtkIdType >(numLocalSuperCells), comm );
+ vtkIdType globalSuperEdges = MpiWrapper::sum(
+ static_cast< vtkIdType >(totalSuperEdges), comm );
+ vtkIdType globalBaseCells = MpiWrapper::sum(
+ static_cast< vtkIdType >(numLocalCells), comm );
+ vtkIdType globalBaseEdges = MpiWrapper::sum(
+ static_cast< vtkIdType >(totalBaseEdges), comm );
+
+ GEOS_LOG_RANK_0( "SUPER-CELL GRAPH:" );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Super-graph nodes: {:>10L}", globalSuperCells ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Super-graph edges: {:>10L}", globalSuperEdges ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Base graph nodes: {:>10L}", globalBaseCells ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Base graph edges: {:>10L}", globalBaseEdges ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Node reduction: {:>10L} cells ({:.1f}%)",
+ globalBaseCells - globalSuperCells,
+ 100.0 * (globalBaseCells - globalSuperCells) / globalBaseCells ) );
+ return std::make_pair( std::move( superGraph ), std::move( superVertexWeights ) );
+}
+
+void validateSuperCellGraph(
+ ArrayOfArrays< pmet_idx_t, pmet_idx_t > const & superGraph,
+ arrayView1d< pmet_idx_t const > const & superElemDist,
+ arrayView1d< pmet_idx_t const > const & vertexWeights,
+ MPI_Comm comm )
+{
+ int const rank = MpiWrapper::commRank( comm );
+ int const numRanks = MpiWrapper::commSize( comm );
+
+ GEOS_LOG_RANK_0( "Running graph integrity checks..." );
+
+ pmet_idx_t myStart = superElemDist[rank];
+ pmet_idx_t myEnd = superElemDist[rank + 1];
+
+ // -----------------------------------------------------------------------
+ // Check 1: No self-loops
+ // -----------------------------------------------------------------------
+ {
+ for( localIndex i = 0; i < superGraph.size(); ++i )
+ {
+ pmet_idx_t myGlobalId = myStart + i;
+
+ for( localIndex j = 0; j < superGraph.sizeOfArray( i ); ++j )
+ {
+ GEOS_ERROR_IF( superGraph[i][j] == myGlobalId,
+ GEOS_FMT( "Rank {}: Super-cell {} (global {}) has self-loop",
+ rank, i, myGlobalId ) );
+ }
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Check 2: All neighbors in valid global range
+ // -----------------------------------------------------------------------
+ {
+ pmet_idx_t globalMin = 0;
+ pmet_idx_t globalMax = superElemDist[numRanks] - 1;
+
+ for( localIndex i = 0; i < superGraph.size(); ++i )
+ {
+ for( localIndex j = 0; j < superGraph.sizeOfArray( i ); ++j )
+ {
+ pmet_idx_t neighbor = superGraph[i][j];
+
+ GEOS_ERROR_IF( neighbor< globalMin || neighbor > globalMax,
+ GEOS_FMT( "Rank {}: Super-cell {} has out-of-range neighbor {} (valid: [{}, {}])",
+ rank, i, neighbor, globalMin, globalMax ) );
+ }
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Check 3: No duplicate edges
+ // -----------------------------------------------------------------------
+ {
+ for( localIndex i = 0; i < superGraph.size(); ++i )
+ {
+ std::set< pmet_idx_t > uniqueNeighbors;
+
+ for( localIndex j = 0; j < superGraph.sizeOfArray( i ); ++j )
+ {
+ pmet_idx_t neighbor = superGraph[i][j];
+
+ GEOS_ERROR_IF( !uniqueNeighbors.insert( neighbor ).second,
+ GEOS_FMT( "Rank {}: Super-cell {} has duplicate neighbor {}",
+ rank, i, neighbor ) );
+ }
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Check 4: Isolated vertices
+ // -----------------------------------------------------------------------
+ {
+ localIndex isolated = 0;
+
+ for( localIndex i = 0; i < superGraph.size(); ++i )
+ {
+ if( superGraph.sizeOfArray( i ) == 0 )
+ {
+ isolated++;
+ if( isolated <= 5 )
+ {
+ pmet_idx_t globalId = myStart + i;
+ GEOS_LOG_RANK( GEOS_FMT( "WARNING: Super-cell {} (global {}) has no neighbors (isolated)",
+ i, globalId ) );
+ }
+ }
+ }
+ GEOS_WARNING_IF( isolated > 0, GEOS_FMT( "Found {} isolated vertices ", isolated ) );
+ }
+
+ // -----------------------------------------------------------------------
+ // Check 5: Verify weights
+ // -----------------------------------------------------------------------
+ {
+ if( !vertexWeights.empty() )
+ {
+ GEOS_ERROR_IF( vertexWeights.size() != superGraph.size(),
+ GEOS_FMT( "Rank {}: Vertex weights size ({}) != graph size ({})",
+ rank, vertexWeights.size(), superGraph.size() ) );
+
+ for( localIndex i = 0; i < vertexWeights.size(); ++i )
+ {
+ GEOS_ERROR_IF( vertexWeights[i] <= 0,
+ GEOS_FMT( "Rank {}: Super-cell {} has invalid weight {}",
+ rank, i, vertexWeights[i] ) );
+ }
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Check 6: Graph symmetry (local edges only)
+ // -----------------------------------------------------------------------
+ GEOS_LOG_RANK_0( "Checking local graph symmetry..." );
+ {
+ std::unordered_set< uint64_t > localOutgoingNeighbors;
+
+ // Build edge set
+ for( localIndex i = 0; i < superGraph.size(); ++i )
+ {
+ pmet_idx_t globalI = myStart + i;
+ auto neighbors = superGraph[i];
+
+ for( localIndex j = 0; j < neighbors.size(); ++j )
+ {
+ pmet_idx_t globalJ = neighbors[j];
+
+ if( globalJ >= myStart && globalJ < myEnd )
+ {
+ uint64_t edgeKey = (static_cast< uint64_t >( globalI ) << 32) |
+ static_cast< uint64_t >( globalJ );
+ localOutgoingNeighbors.insert( edgeKey );
+ }
+ }
+ }
+
+ // Verify symmetry
+ for( localIndex i = 0; i < superGraph.size(); ++i )
+ {
+ pmet_idx_t globalI = myStart + i;
+ auto neighbors = superGraph[i];
+
+ for( localIndex j = 0; j < neighbors.size(); ++j )
+ {
+ pmet_idx_t globalJ = neighbors[j];
+
+ if( globalJ >= myStart && globalJ < myEnd )
+ {
+ uint64_t reverseKey = (static_cast< uint64_t >( globalJ ) << 32) |
+ static_cast< uint64_t >( globalI );
+
+ GEOS_ERROR_IF( localOutgoingNeighbors.find( reverseKey ) == localOutgoingNeighbors.end(),
+ GEOS_FMT( "Rank {}: Asymmetric edge ({} -> {}) - reverse edge missing",
+ rank, globalI, globalJ ) );
+ }
+ }
+ }
+ }
+
+}
+
+
+// =============================================================================================
+// SECTION 5: PARTITION UNPACKING (after ParMETIS)
+// =============================================================================================
+array1d< int64_t >
+unpackSuperCellPartitioning(
+ vtkSmartPointer< vtkUnstructuredGrid > cells3D,
+ array1d< int64_t > const & superPartitioning,
+ stdMap< vtkIdType, localIndex > const & superCellIdToLocalIdx,
+ MPI_Comm comm )
+{
+ int const rank = MpiWrapper::commRank( comm );
+
+ GEOS_LOG_RANK_0( " UNPACKING SUPER-CELL PARTITIONING" );
+
+ // -----------------------------------------------------------------------
+ // Step 1: Get super-cell ID array
+ // -----------------------------------------------------------------------
+ vtkIdTypeArray * superCellIdArray =
+ vtkIdTypeArray::SafeDownCast( cells3D->GetCellData()->GetArray( "SuperCellId" ) );
+
+ GEOS_ERROR_IF( superCellIdArray == nullptr,
+ GEOS_FMT( "Rank {}: SuperCellId array not found", rank ) );
+
+ GEOS_ERROR_IF( static_cast< size_t >(superPartitioning.size()) != superCellIdToLocalIdx.size(),
+ GEOS_FMT( "Rank {}: Super-cell partitioning size ({}) doesn't match number of local super-cells ({})",
+ rank, superPartitioning.size(), superCellIdToLocalIdx.size() ) );
+
+ // -----------------------------------------------------------------------
+ // Step 2: Assign cells to ranks based on super-cell partitioning
+ // -----------------------------------------------------------------------
+ vtkIdType const numCells = cells3D->GetNumberOfCells();
+ array1d< int64_t > cellPartitioning( numCells );
+
+ for( vtkIdType i = 0; i < numCells; ++i )
+ {
+ vtkIdType superCellId = superCellIdArray->GetValue( i );
+ auto it = superCellIdToLocalIdx.find( superCellId );
+
+ GEOS_ERROR_IF( it == superCellIdToLocalIdx.end(),
+ GEOS_FMT( "Rank {}: Cell {} has unknown super-cell ID {}",
+ rank, i, superCellId ) );
+
+ cellPartitioning[i] = superPartitioning[it->second];
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 3: Validate that super-cells weren't split
+ // -----------------------------------------------------------------------
+ stdMap< vtkIdType, std::set< int64_t > > superCellToRanks;
+
+ for( vtkIdType i = 0; i < numCells; ++i ) // ← Reuse numCells (no redeclaration)
+ {
+ vtkIdType const scId = superCellIdArray->GetValue( i );
+ superCellToRanks.get_inserted( scId ).insert( cellPartitioning[i] );
+ }
+
+ vtkIdType numSplitSuperCells = 0;
+ for( auto const & [superCellId, ranks] : superCellToRanks )
+ {
+ if( ranks.size() > 1 )
+ {
+ ++numSplitSuperCells;
+ }
+ }
+
+ vtkIdType totalSplitSuperCells = MpiWrapper::sum( numSplitSuperCells, comm );
+
+ GEOS_ERROR_IF( totalSplitSuperCells > 0,
+ GEOS_FMT( "Partitioning failed: {:L} super-cells were split across ranks!",
+ totalSplitSuperCells ) );
+
+ GEOS_LOG_RANK_0( "Super-cell partitioning validated successfully" );
+
+ return cellPartitioning;
+}
+
+} // namespace vtk
+} // namespace geos
diff --git a/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.hpp b/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.hpp
new file mode 100644
index 00000000000..565140c9651
--- /dev/null
+++ b/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.hpp
@@ -0,0 +1,191 @@
+/**
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+#ifndef GEOS_MESH_GENERATORS_VTKSUPERCELLPARTITIONING_HPP
+#define GEOS_MESH_GENERATORS_VTKSUPERCELLPARTITIONING_HPP
+
+
+#include "common/MpiWrapper.hpp"
+#include "common/TimingMacros.hpp"
+#include "mesh/generators/VTKMeshGeneratorTools.hpp" // for PartitionMethod
+
+
+#include
+
+// GEOS types
+#include "common/DataTypes.hpp"
+#include "LvArray/src/ArrayOfArrays.hpp"
+
+// VTK forward declarations
+class vtkUnstructuredGrid;
+class vtkDataSet;
+
+// Get pmet_idx_t definition from ParMETIS interface
+#include "mesh/generators/ParMETISInterface.hpp" // Defines pmet_idx_t as int64_t
+
+
+
+namespace geos
+{
+
+namespace vtk
+{
+
+/**
+ * @brief Strategy for initial super-cell scatter across ranks (before ParMETIS refinement)
+ */
+enum class InitialDistributionStrategy
+{
+ MORTON, /// Distribute by Morton Z-curve for spatial locality
+ BLOCK /// Distribute by simple contiguous blocks (faster)
+};
+
+/**
+ * @brief Metadata about super-cells for partitioning
+ *
+ * A super-cell is a group of 3D cells that must stay together during partitioning.
+ * Typically these are cells connected by fractures.
+ */
+struct SuperCellInfo
+{
+ /// Map: SuperCellId -> vector of global cell IDs in that super-cell
+ std::map< vtkIdType, std::vector< vtkIdType > > superCellToOriginalCells;
+
+ /// Map: SuperCellId -> weight (number of cells in super-cell)
+ std::map< vtkIdType, vtkIdType > vertexWeights;
+
+ /// Set of SuperCellIds that contain multiple cells (atomic units)
+ std::set< vtkIdType > atomicSuperCells;
+};
+
+/**
+ * @brief Tag 3D cells with super-cell IDs based on fracture connectivity
+ *
+ * Creates a "SuperCellId" cell data array where:
+ * - Cells connected by fractures share the same super-cell ID
+ * - Regular cells have their own unique super-cell ID (= their global ID)
+ *
+ * This runs on rank 0 only, using the fracture neighbor information.
+ *
+ * @param cells3D The 3D volumetric cells (modified in-place to add SuperCellId array)
+ * @param fractureNeighbors Map of fracture name to neighbor mapping (fracture element -> 3D cell neighbors)
+ * @param fractureWeight Additional weight applied to fracture super-cells for load balancing
+ * @return Super-cell metadata for partitioning
+ */
+SuperCellInfo tagCellsWithSuperCellIds(
+ vtkSmartPointer< vtkUnstructuredGrid > cells3D,
+ stdMap< string, ArrayOfArrays< vtkIdType, int64_t > > const & fractureNeighbors,
+ integer fractureWeight );
+
+/**
+ * @brief Reconstruct super-cell info from the SuperCellId array
+ *
+ * After mesh redistribution, each rank needs to rebuild its local super-cell metadata
+ * from the SuperCellId cell data array. Fracture super-cells are weighted to improve
+ * load balancing by accounting for their higher computational cost.
+ *
+ * @param mesh The distributed mesh with SuperCellId array
+ * @param fractureWeight Additional weight added to fracture super-cells (cells.size() > 1)
+ * to account for contact mechanics computation cost.
+ * @return Local super-cell metadata with weighted super-cells
+ */
+SuperCellInfo reconstructSuperCellInfo( vtkSmartPointer< vtkUnstructuredGrid > mesh,
+ integer fractureWeight );
+
+/**
+ * @brief Initial redistribution preserving super-cell integrity
+ * @param cells3D Input mesh (only non-empty on rank 0)
+ * @param comm MPI communicator
+ * @param strategy Initial distribution strategy:
+ * - MORTON: Sort super-cells by Morton Z-curve for spatial locality
+ * - BLOCK: Simple contiguous block assignment
+ * @return Redistributed mesh with SuperCellId array preserved
+ *
+ * Uses simple round-robin assignment of super-cells to ranks.
+ * This is faster than graph partitioning and suitable for initial distribution.
+ */
+vtkSmartPointer< vtkDataSet >
+redistributeBySuperCellBlocks( vtkSmartPointer< vtkUnstructuredGrid > cells3D,
+ MPI_Comm comm,
+ InitialDistributionStrategy strategy = InitialDistributionStrategy::MORTON );
+
+
+/**
+ * @brief Build a graph where nodes are super-cells (not individual cells)
+ *
+ * Each super-cell becomes a single node in the graph. Edges connect super-cells
+ * whose constituent cells are neighbors in the original mesh.
+ *
+ * @param cells3D The tagged 3D mesh with SuperCellId array
+ * @param baseGraph The original cell-to-cell adjacency graph
+ * @param baseElemDist Element distribution for base graph (numRanks+1 array)
+ * @param info Super-cell metadata
+ * @param comm MPI communicator
+ * @return Pair of (super-cell graph, super-cell vertex weights)
+ */
+std::pair< ArrayOfArrays< pmet_idx_t, pmet_idx_t >, array1d< pmet_idx_t > >
+buildSuperCellGraph(
+ vtkSmartPointer< vtkUnstructuredGrid > cells3D,
+ ArrayOfArrays< pmet_idx_t, pmet_idx_t > const & baseGraph,
+ arrayView1d< pmet_idx_t const > const & baseElemDist,
+ SuperCellInfo const & info,
+ MPI_Comm comm );
+
+/**
+ * @brief Validate super-cell graph integrity before partitioning
+ *
+ * Checks for:
+ * - Self-loops
+ * - Out-of-range neighbor indices
+ * - Duplicate edges
+ * - Isolated vertices
+ * - Invalid vertex weights
+ *
+ * @param superGraph The super-cell adjacency graph
+ * @param superElemDist Element distribution array for super-cells
+ * @param vertexWeights Vertex weights for load balancing
+ * @param comm MPI communicator
+ * @throws If graph has integrity errors
+ */
+void validateSuperCellGraph(
+ ArrayOfArrays< pmet_idx_t, pmet_idx_t > const & superGraph,
+ arrayView1d< pmet_idx_t const > const & superElemDist,
+ arrayView1d< pmet_idx_t const > const & vertexWeights,
+ MPI_Comm comm );
+
+/**
+ * @brief Unpack super-cell partitioning to individual cells
+ *
+ * Maps the super-cell → rank assignments from ParMETIS back to individual cell assignments.
+ * All cells in the same super-cell get assigned to the same rank.
+ *
+ * @param cells3D The 3D mesh with SuperCellId array
+ * @param superPartitioning Super-cell → rank assignments from ParMETIS
+ * @param superCellIdToLocalIdx Mapping from SuperCellId to local super-cell index
+ * @param comm MPI communicator
+ * @return Cell-level partitioning array (cell → rank)
+ */
+array1d< int64_t >
+unpackSuperCellPartitioning(
+ vtkSmartPointer< vtkUnstructuredGrid > cells3D,
+ array1d< int64_t > const & superPartitioning,
+ stdMap< vtkIdType, localIndex > const & superCellIdToLocalIdx,
+ MPI_Comm comm );
+
+} // namespace vtk
+
+} // namespace geos
+
+#endif /* GEOS_MESH_GENERATORS_VTKSUPERCELLPARTITIONING_HPP */
diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp
index e1da6189ac3..19b1c8b6052 100644
--- a/src/coreComponents/mesh/generators/VTKUtilities.cpp
+++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp
@@ -22,6 +22,8 @@
#include "mesh/generators/VTKMeshGeneratorTools.hpp"
#include "mesh/generators/VTKUtilities.hpp"
#include "mesh/MeshFields.hpp"
+#include "mesh/generators/VTKSuperCellPartitioning.hpp"
+
#ifdef GEOS_USE_PARMETIS
#include "mesh/generators/ParMETISInterface.hpp"
@@ -66,6 +68,7 @@
#include
#include
#include
+#include
#ifdef GEOS_USE_MPI
#include
@@ -236,59 +239,38 @@ vtkSmartPointer< vtkCellArray > getCellArray( vtkSmartPointer< vtkDataSet > mesh
/**
- * @brief Build the element to nodes mappings for all the @p meshes.
- * @tparam INDEX_TYPE The indexing type that will be used by the toolbox that will perfomrn the parallel split.
+ * @brief Build the element to nodes mappings implementation
+ * @tparam INDEX_TYPE The indexing type
* @tparam POLICY The computational policy (parallel/serial)
- * @param meshes All the meshes involved (volumic and surfacic (for fractures))
- * @param cells The vtk cell array.
- * @return The mapping.
+ * @param mesh The 3D mesh
+ * @param cells The vtk cell array
+ * @return The mapping for 3D cells only
*/
template< typename INDEX_TYPE, typename POLICY >
ArrayOfArrays< INDEX_TYPE, INDEX_TYPE >
-buildElemToNodesImpl( AllMeshes & meshes,
+buildElemToNodesImpl( vtkSmartPointer< vtkDataSet > mesh,
vtkSmartPointer< vtkCellArray > const & cells )
{
- localIndex const num3dCells = LvArray::integerConversion< localIndex >( meshes.getMainMesh()->GetNumberOfCells() );
+ GEOS_MARK_FUNCTION;
+
+ localIndex const numCells = LvArray::integerConversion< localIndex >( mesh->GetNumberOfCells() );
- localIndex num2dCells = 0;
- stdMap< string, CollocatedNodes > collocatedNodesMap;
- for( auto & [fractureName, fractureMesh]: meshes.getFaceBlocks() )
- {
- num2dCells += fractureMesh->GetNumberOfCells();
- collocatedNodesMap.insert( { fractureName, CollocatedNodes( fractureName, fractureMesh ) } );
- }
- localIndex const numCells = num3dCells + num2dCells;
array1d< INDEX_TYPE > nodeCounts( numCells );
// GetCellSize() is always thread-safe, can run in parallel
- forAll< parallelHostPolicy >( num3dCells, [nodeCounts = nodeCounts.toView(), &cells] ( localIndex const cellIdx )
+ forAll< parallelHostPolicy >( numCells, [nodeCounts = nodeCounts.toView(), &cells] ( localIndex const cellIdx )
{
nodeCounts[cellIdx] = LvArray::integerConversion< INDEX_TYPE >( cells->GetCellSize( cellIdx ) );
} );
- localIndex offset = num3dCells;
- for( auto & [fractureName, fractureMesh]: meshes.getFaceBlocks() )
- {
- CollocatedNodes const & collocatedNodes = collocatedNodesMap.at( fractureName );
- forAll< parallelHostPolicy >( fractureMesh->GetNumberOfCells(), [&, nodeCounts = nodeCounts.toView(), fracture = fractureMesh.Get()] ( localIndex const cellIdx )
- {
- nodeCounts[cellIdx + offset] = 0;
- // We are doing a very strict allocation because some TPLs rely on not having any over allocation.
- for( vtkIdType const pointId: *fracture->GetCell( cellIdx )->GetPointIds() )
- {
- nodeCounts[cellIdx + offset] += collocatedNodes[pointId].size();
- }
- } );
- offset += fractureMesh->GetNumberOfCells();
- }
-
+ // Create strictly allocated ArrayOfArrays using resizeFromCapacities
ArrayOfArrays< INDEX_TYPE, INDEX_TYPE > elemToNodes;
elemToNodes.template resizeFromCapacities< parallelHostPolicy >( numCells, nodeCounts.data() );
- vtkIdTypeArray const & globalPointId = *vtkIdTypeArray::FastDownCast( meshes.getMainMesh()->GetPointData()->GetGlobalIds() );
+ vtkIdTypeArray const & globalPointId = *vtkIdTypeArray::FastDownCast( mesh->GetPointData()->GetGlobalIds() );
// GetCellAtId() is conditionally thread-safe, use POLICY argument
- forAll< POLICY >( num3dCells, [&cells, &globalPointId, elemToNodes = elemToNodes.toView()] ( localIndex const cellIdx )
+ forAll< POLICY >( numCells, [&cells, &globalPointId, elemToNodes = elemToNodes.toView()] ( localIndex const cellIdx )
{
vtkIdType numPts;
vtkIdType const * points;
@@ -300,44 +282,27 @@ buildElemToNodesImpl( AllMeshes & meshes,
}
} );
- offset = num3dCells; // Restarting the loop from the beginning.
- for( auto & [fractureName, fractureMesh]: meshes.getFaceBlocks() )
- {
- CollocatedNodes const & collocatedNodes = collocatedNodesMap.at( fractureName );
- for( vtkIdType i = 0; i < fractureMesh->GetNumberOfCells(); ++i )
- {
- for( vtkIdType const pointId: *fractureMesh->GetCell( i )->GetPointIds() )
- {
- for( vtkIdType const & tmp: collocatedNodes[pointId] )
- {
- elemToNodes.emplaceBack( offset + i, tmp );
- }
- }
- }
- offset += fractureMesh->GetNumberOfCells();
- }
-
return elemToNodes;
}
-
/**
- * @brief Build the element to nodes mappings for all the @p meshes.
- * @tparam INDEX_TYPE The indexing type that will be used by the toolbox that will perfomrn the parallel split.
- * @param meshes All the meshes involved (volumic and surfacic (for fractures))l
- * @return The mapping.
+ * @brief Build the element to nodes mappings for the 3D mesh only
+ * @tparam INDEX_TYPE The indexing type
+ * @param mesh The 3D mesh
+ * @return The mapping for 3D cells only (fractures excluded)
+ * @note Fractures are partitioned separately based on 3D neighbors, not included in ParMETIS graph
*/
template< typename INDEX_TYPE >
ArrayOfArrays< INDEX_TYPE, INDEX_TYPE >
-buildElemToNodes( AllMeshes & meshes )
+buildElemToNodes( vtkSmartPointer< vtkDataSet > mesh )
{
- vtkSmartPointer< vtkCellArray > const & cells = vtk::getCellArray( meshes.getMainMesh() );
+ vtkSmartPointer< vtkCellArray > const & cells = vtk::getCellArray( mesh );
// According to VTK docs, IsStorageShareable() indicates whether pointers extracted via
// vtkCellArray::GetCellAtId() are pointers into internal storage rather than temp buffer
// and thus results can be used in a thread-safe way.
return cells->IsStorageShareable()
- ? buildElemToNodesImpl< INDEX_TYPE, parallelHostPolicy >( meshes, cells )
- : buildElemToNodesImpl< INDEX_TYPE, serialPolicy >( meshes, cells );
+ ? buildElemToNodesImpl< INDEX_TYPE, parallelHostPolicy >( mesh, cells )
+ : buildElemToNodesImpl< INDEX_TYPE, serialPolicy >( mesh, cells );
}
/**
@@ -559,23 +524,23 @@ AllMeshes loadAllMeshes( Path const & filePath,
string const & mainBlockName,
string_array const & faceBlockNames )
{
- int const lastRank = MpiWrapper::commSize() - 1;
vtkSmartPointer< vtkDataSet > main = loadMesh( filePath, mainBlockName );
stdMap< string, vtkSmartPointer< vtkDataSet > > faces;
+ // Load fractures on rank 0 (same as main mesh for 2D cells)
+ // This allows building fracture-to-3D connectivity before redistribution
for( string const & faceBlockName: faceBlockNames )
{
- faces.insert( { faceBlockName, loadMesh( filePath, faceBlockName, lastRank ) } );
+ faces.insert( { faceBlockName, loadMesh( filePath, faceBlockName, 0 ) } );
}
return AllMeshes( main, faces );
}
-
/**
- * @brief Partition the mesh using cell graph methods (ParMETIS or PTScotch)
+ * @brief Partition the 3D mesh using cell graph methods (ParMETIS or PTScotch)
*
- * @param[in] input a collection of vtk main (3D) and fracture mesh
+ * @param[in] mesh3D the 3D main mesh to partition
* @param[in] method the partitioning method
* @param[in] comm the MPI communicator
* @param[in] numParts the number of partitions
@@ -584,7 +549,7 @@ AllMeshes loadAllMeshes( Path const & filePath,
* @return the cell partitioning array
*/
array1d< int64_t >
-partitionByCellGraph( AllMeshes & input,
+partitionByCellGraph( vtkSmartPointer< vtkDataSet > mesh3D,
PartitionMethod const method,
MPI_Comm const comm,
int const numParts,
@@ -593,16 +558,18 @@ partitionByCellGraph( AllMeshes & input,
{
GEOS_MARK_FUNCTION;
- pmet_idx_t const numElems = input.getMainMesh()->GetNumberOfCells();
+ pmet_idx_t const numElems = mesh3D->GetNumberOfCells();
pmet_idx_t const numRanks = MpiWrapper::commSize( comm );
- int const rank = MpiWrapper::commRank( comm );
- int const lastRank = numRanks - 1;
// Value at each index (i.e. MPI rank) of `elemDist` gives the first element index of the MPI rank.
- // It's assumed that MPI ranks spans continuous numbers of elements.
+ // It's assumed that MPI ranks span continuous numbers of elements.
// Thus, the number of elements of each rank can be deduced by subtracting
// the values between two consecutive ranks. To be able to do this even for the last rank,
// a last additional value is appended, and the size of the array is then the comm size plus 1.
+ //
+ // Note: elemDist contains the distribution of 3D cells only.
+ // Fractures are NOT included in the ParMETIS graph - they will be assigned separately
+ // in merge2D3DCellsAndRedistribute() based on their 3D neighbor connectivity.
array1d< pmet_idx_t > const elemDist( numRanks + 1 );
{
array1d< pmet_idx_t > elemCounts;
@@ -610,21 +577,9 @@ partitionByCellGraph( AllMeshes & input,
std::partial_sum( elemCounts.begin(), elemCounts.end(), elemDist.begin() + 1 );
}
- vtkIdType localNumFracCells = 0;
- if( rank == lastRank ) // Let's add artificially the fracture to the last rank (for numbering reasons).
- {
- // Adding one fracture element
- for( auto const & [fractureName, fracture]: input.getFaceBlocks() )
- {
- localNumFracCells += fracture->GetNumberOfCells();
- }
- }
- vtkIdType globalNumFracCells = localNumFracCells;
- MpiWrapper::broadcast( globalNumFracCells, lastRank, comm );
- elemDist[lastRank + 1] += globalNumFracCells;
+ // Build element-to-node connectivity for the 3D mesh only
+ ArrayOfArrays< pmet_idx_t, pmet_idx_t > const elemToNodes = buildElemToNodes< pmet_idx_t >( mesh3D );
- // The `elemToNodes` mapping binds element indices (local to the rank) to the global indices of their support nodes.
- ArrayOfArrays< pmet_idx_t, pmet_idx_t > const elemToNodes = buildElemToNodes< pmet_idx_t >( input );
ArrayOfArrays< pmet_idx_t, pmet_idx_t > graph;
#ifdef GEOS_USE_PARMETIS
graph = parmetis::meshToDual( elemToNodes.toViewConst(), elemDist, comm, minCommonNodes );
@@ -664,57 +619,275 @@ partitionByCellGraph( AllMeshes & input,
}
/**
- * @brief Redistribute the mesh using cell graph methods (ParMETIS or PTScotch)
- * @param[in] mesh a vtk grid
- * @param[in] method the partitioning method
- * @param[in] comm the MPI communicator
- * @param[in] numRefinements the number of refinements for PTScotch
- * @return
+ * @brief Redistribute mesh preserving super-cell integrity using weighted graph partitioning
+ *
+ * This function partitions a mesh while ensuring that groups of cells marked with the same
+ * SuperCellId remain together on the same MPI rank. It:
+ * 1. Builds a base cell-to-cell adjacency graph
+ * 2. Collapses cells into super-cell vertices (weighted by cell count)
+ * 3. Partitions the super-cell graph using ParMETIS/PTScotch
+ * 4. Unpacks assignments to individual cells
+ * 5. Redistributes the mesh
+ *
+ * @param mesh Input mesh with "SuperCellId" cell data array
+ * @param method Partitioning algorithm (parmetis or ptscotch)
+ * @param comm MPI communicator
+ * @param numRefinementIterations Number of ParMETIS refinement passes
+ * @return Redistributed mesh with super-cells kept intact
+ *
+ * @pre mesh must be vtkUnstructuredGrid
+ * @pre mesh must have "SuperCellId" integer cell data array
+ * @pre All cells in a super-cell must be topologically connected
+ * @post No super-cell is split across multiple ranks
*/
-AllMeshes
-redistributeByCellGraph( AllMeshes & input,
- PartitionMethod const method,
- MPI_Comm const comm,
- int const numRefinements )
+vtkSmartPointer< vtkDataSet >
+redistributeBySuperCellGraph(
+ vtkSmartPointer< vtkDataSet > mesh,
+ PartitionMethod const method,
+ MPI_Comm comm,
+ int const numRefinementIterations,
+ int const fractureWeight )
{
GEOS_MARK_FUNCTION;
int const rank = MpiWrapper::commRank( comm );
int const numRanks = MpiWrapper::commSize( comm );
- array1d< int64_t > newPartitions = partitionByCellGraph( input, method, comm, numRanks, 3, numRefinements );
- // Extract the partition information related to the fracture mesh.
- stdMap< string, array1d< pmet_idx_t > > newFracturePartitions;
- vtkIdType fracOffset = input.getMainMesh()->GetNumberOfCells();
- vtkIdType localNumFracCells = 0;
- for( auto const & [fractureName, fracture]: input.getFaceBlocks() )
+ // -----------------------------------------------------------------------
+ // Step 1: Build base cell graph (standard adjacency)
+ // -----------------------------------------------------------------------
+ vtkSmartPointer< vtkUnstructuredGrid > ugrid =
+ vtkUnstructuredGrid::SafeDownCast( mesh );
+
+ GEOS_ERROR_IF( !ugrid, "Rank " << rank << ": Mesh is not vtkUnstructuredGrid" );
+
+ ArrayOfArrays< pmet_idx_t, pmet_idx_t > baseCellGraph;
+ array1d< pmet_idx_t > baseElemDist( numRanks + 1 );
{
- localIndex const numFracCells = fracture->GetNumberOfCells();
- localNumFracCells += (rank == numRanks - 1) ? fracture->GetNumberOfCells() : 0;
- array1d< pmet_idx_t > tmp( numFracCells );
- std::copy( newPartitions.begin() + fracOffset, newPartitions.begin() + fracOffset + numFracCells, tmp.begin() );
- newFracturePartitions.insert( { fractureName, tmp } );
- fracOffset += numFracCells;
+ // Build element distribution based on local cell counts
+ pmet_idx_t const numLocalCells = ugrid->GetNumberOfCells();
+
+ array1d< pmet_idx_t > cellCounts;
+ MpiWrapper::allGather( numLocalCells, cellCounts, comm );
+
+ baseElemDist[0] = 0;
+ std::partial_sum( cellCounts.begin(), cellCounts.end(), baseElemDist.begin() + 1 );
+
+ // Build element-to-node connectivity local to each rank
+ ArrayOfArrays< pmet_idx_t, pmet_idx_t > const elemToNodes =
+ buildElemToNodes< pmet_idx_t >( ugrid );
+
+ // Build dual graph (cell-to-cell via shared nodes)
+#ifdef GEOS_USE_PARMETIS
+ int const minCommonNodes = 3; // Minimum shared nodes for edge
+ baseCellGraph = parmetis::meshToDual( elemToNodes.toViewConst(), baseElemDist, comm, minCommonNodes );
+#else
+ GEOS_THROW( "GEOS must be built with ParMETIS support (ENABLE_PARMETIS=ON) "
+ "to build cell graphs for partitioning", InputError );
+#endif
}
- // Now do the same for the 3d mesh, simply by trimming the fracture information.
- newPartitions.resize( newPartitions.size() - localNumFracCells );
- // Now, perform the final steps: first, a new split following the new partitions.
- // Then those newly split meshes will be redistributed across the ranks.
+ // -----------------------------------------------------------------------
+ // Step 2: Reconstruct super-cell info on all ranks (from SuperCellId array)
+ // -----------------------------------------------------------------------
+ SuperCellInfo localSuperCellInfo = reconstructSuperCellInfo( ugrid, fractureWeight );
+
+ // -----------------------------------------------------------------------
+ // Step 3: Build super-cell graph
+ // -----------------------------------------------------------------------
+ auto [superCellGraph, superVertexWeights] = buildSuperCellGraph(
+ ugrid,
+ baseCellGraph,
+ baseElemDist,
+ localSuperCellInfo,
+ comm
+ );
- // First for the main 3d mesh...
- vtkSmartPointer< vtkPartitionedDataSet > const splitMesh = splitMeshByPartition( input.getMainMesh(), numRanks, newPartitions.toViewConst() );
- vtkSmartPointer< vtkUnstructuredGrid > finalMesh = vtk::redistribute( *splitMesh, MPI_COMM_GEOS );
- // ... and then for the fractures.
- stdMap< string, vtkSmartPointer< vtkDataSet > > finalFractures;
- for( auto const & [fractureName, fracture]: input.getFaceBlocks() )
+ // -----------------------------------------------------------------------
+ // Step 4: Compute super-cell element distribution
+ // -----------------------------------------------------------------------
+ array1d< pmet_idx_t > superElemDist( numRanks + 1 );
{
- vtkSmartPointer< vtkPartitionedDataSet > const splitFracMesh = splitMeshByPartition( fracture, numRanks, newFracturePartitions[fractureName].toViewConst() );
- vtkSmartPointer< vtkUnstructuredGrid > const finalFracMesh = vtk::redistribute( *splitFracMesh, MPI_COMM_GEOS );
- finalFractures.insert( {fractureName, finalFracMesh} );
+ pmet_idx_t const localSuperCellCount = LvArray::integerConversion< pmet_idx_t >( superCellGraph.size() );
+
+ array1d< pmet_idx_t > superCellCounts;
+ MpiWrapper::allGather( localSuperCellCount, superCellCounts, comm );
+
+ superElemDist[0] = 0;
+ std::partial_sum( superCellCounts.begin(), superCellCounts.end(),
+ superElemDist.begin() + 1 );
}
- return AllMeshes( finalMesh, finalFractures );
+ // -----------------------------------------------------------------------
+ // Step 5: Validate graph before partitioning
+ // -----------------------------------------------------------------------
+ validateSuperCellGraph(
+ superCellGraph,
+ superElemDist.toViewConst(),
+ superVertexWeights.toViewConst(),
+ comm
+ );
+
+ // -----------------------------------------------------------------------
+ // Step 6: Partition super-cell graph using ParMETIS/PTScotch
+ // -----------------------------------------------------------------------
+ array1d< int64_t > superCellPartitioning;
+
+ if( method == PartitionMethod::parmetis )
+ {
+ GEOS_LOG_RANK_0( "Partitioning super-cell graph with ParMETIS..." );
+
+ // Basic validation
+ pmet_idx_t const minGraphSize = MpiWrapper::min( static_cast< pmet_idx_t >( superCellGraph.size() ), comm );
+ GEOS_ERROR_IF( minGraphSize == 0, "At least one rank has an empty super-cell graph!" );
+
+ superCellPartitioning = parmetis::partitionWeighted(
+ superCellGraph.toViewConst(),
+ superVertexWeights.toViewConst(),
+ superElemDist.toViewConst(),
+ numRanks,
+ comm,
+ numRefinementIterations
+ );
+
+ GEOS_LOG_RANK_0( "ParMETIS completed successfully!" );
+ }
+ else
+ {
+ GEOS_ERROR( GEOS_FMT( "Unsupported partition method: {}",
+ EnumStrings< PartitionMethod >::toString( method ) ) );
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 7: Build mapping from SuperCellId to local super-cell index (FIXED)
+ // -----------------------------------------------------------------------
+
+ vtkIdTypeArray * superCellIdArray =
+ vtkIdTypeArray::SafeDownCast( ugrid->GetCellData()->GetArray( "SuperCellId" ) );
+
+ GEOS_ERROR_IF( !superCellIdArray,
+ GEOS_FMT( "SuperCellId array not found on rank {}", rank ) );
+
+ stdMap< vtkIdType, localIndex > superCellIdToLocalIdx;
+
+// Build ordered list of local super-cell IDs
+ stdVector< vtkIdType > orderedSuperCellIds;
+ orderedSuperCellIds.reserve( localSuperCellInfo.superCellToOriginalCells.size() );
+
+ for( auto const & [scId, cells] : localSuperCellInfo.superCellToOriginalCells )
+ {
+ orderedSuperCellIds.push_back( scId );
+ }
+
+ // Sort to ensure consistent ordering
+ std::sort( orderedSuperCellIds.begin(), orderedSuperCellIds.end() );
+
+ localIndex const numLocalSuperCells = LvArray::integerConversion< localIndex >( orderedSuperCellIds.size() );
+ for( localIndex i = 0; i < numLocalSuperCells; ++i )
+ {
+ superCellIdToLocalIdx.insert( { orderedSuperCellIds[i], i } );
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 8: Unpack super-cell partitioning to individual cells
+ // -----------------------------------------------------------------------
+ array1d< int64_t > cellPartitioning = unpackSuperCellPartitioning(
+ ugrid,
+ superCellPartitioning,
+ superCellIdToLocalIdx,
+ comm
+ );
+
+ // -----------------------------------------------------------------------
+ // Step 9: Verify no super-cells were split across ranks
+ // -----------------------------------------------------------------------
+ stdMap< vtkIdType, int64_t > superCellToRank;
+
+ vtkIdType const numCells = ugrid->GetNumberOfCells();
+ for( vtkIdType i = 0; i < numCells; ++i )
+ {
+ vtkIdType const scId = superCellIdArray->GetValue( i );
+ int64_t const targetRank = cellPartitioning[i];
+
+ auto [it, inserted] = superCellToRank.insert( {scId, targetRank} );
+
+ GEOS_ERROR_IF( !inserted && it->second != targetRank,
+ GEOS_FMT( "Super-cell {} is split: assigned to both rank {} and rank {}",
+ scId, it->second, targetRank ) );
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 10: Redistribute mesh according to cell partitioning
+ // -----------------------------------------------------------------------
+ // Split mesh according to partitioning
+ vtkSmartPointer< vtkPartitionedDataSet > splitMesh =
+ splitMeshByPartition( ugrid, numRanks, cellPartitioning.toViewConst() );
+
+ // Redistribute using VTK
+ vtkSmartPointer< vtkDataSet > redistributed = vtk::redistribute( *splitMesh, comm );
+
+ // -----------------------------------------------------------------------
+ // Step 11: Report final distribution statistics
+ // -----------------------------------------------------------------------
+ array1d< vtkIdType > cellsPerRank;
+ vtkIdType const localCells = redistributed->GetNumberOfCells();
+ MpiWrapper::allGather( localCells, cellsPerRank, comm );
+
+ if( rank == 0 )
+ {
+ vtkIdType totalCells = 0;
+ vtkIdType minCells = std::numeric_limits< vtkIdType >::max();
+ vtkIdType maxCells = 0;
+
+ for( int r = 0; r < numRanks; ++r )
+ {
+ totalCells += cellsPerRank[r];
+ minCells = std::min( minCells, cellsPerRank[r] );
+ maxCells = std::max( maxCells, cellsPerRank[r] );
+ }
+
+ real64 const avgCells = static_cast< real64 >( totalCells ) / numRanks;
+ real64 const imbalance = (numRanks > 1) ? (maxCells - avgCells) / avgCells * 100.0 : 0.0;
+
+ GEOS_LOG_RANK_0( "\nFinal 3D cell distribution:" );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Total: {:10} cells", totalCells ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Average: {:10} cells", static_cast< vtkIdType >( avgCells )) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Range: [{}, {}]", minCells, maxCells ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Imbalance: {:.2f}%", imbalance ) );
+ }
+
+ return redistributed;
+}
+
+
+/**
+ * @brief Redistribute the 3D mesh using cell graph methods (ParMETIS or PTScotch)
+ * @param[in] mesh3D The 3D main mesh to partition
+ * @param[in] method The partitioning method
+ * @param[in] comm The MPI communicator
+ * @param[in] numRefinements The number of refinement iterations
+ * @return Redistributed 3D mesh
+ * @note Fractures are NOT partitioned here - they will be assigned later
+ * in merge2D3DCellsAndRedistribute() based on 3D neighbor connectivity
+ */
+vtkSmartPointer< vtkDataSet >
+redistributeByCellGraph( vtkSmartPointer< vtkDataSet > mesh3D,
+ PartitionMethod const method,
+ MPI_Comm const comm,
+ int const numRefinements )
+{
+ GEOS_MARK_FUNCTION;
+
+ int const numRanks = MpiWrapper::commSize( comm );
+
+ // Partition the 3D main mesh only
+ array1d< int64_t > newPartitions = partitionByCellGraph( mesh3D, method, comm, numRanks, 3, numRefinements );
+
+ // Split and redistribute the 3D mesh
+ vtkSmartPointer< vtkPartitionedDataSet > const splitMesh =
+ splitMeshByPartition( mesh3D, numRanks, newPartitions.toViewConst() );
+
+ return vtk::redistribute( *splitMesh, comm );
}
@@ -820,6 +993,656 @@ scatterByBlock( vtkDataSet & mesh )
return result;
}
+/**
+ * @brief Separate 2D and 3D cells from a mesh
+ *
+ * Extracts cells by topological dimension into separate unstructured grids while
+ * maintaining mappings to original cell indices. 1D and 0D cells are ignored.
+ *
+ * @param[in] mesh Input mesh containing cells of mixed dimensions
+ * @param[out] cells3D Unstructured grid containing only 3D volumetric cells
+ * @param[out] cells2D Unstructured grid containing only 2D surface cells
+ * @param[out] cells3DToOriginal Index mapping from extracted 3D cells to original mesh
+ * @param[out] cells2DToOriginal Index mapping from extracted 2D cells to original mesh
+ *
+ * @note The output grids are deep copies and independent of the input mesh
+ */
+static void separateCellsByDimension( vtkDataSet & mesh,
+ vtkSmartPointer< vtkUnstructuredGrid > & cells3D,
+ vtkSmartPointer< vtkUnstructuredGrid > & cells2D,
+ array1d< vtkIdType > & cells3DToOriginal,
+ array1d< vtkIdType > & cells2DToOriginal )
+{
+ GEOS_MARK_FUNCTION;
+
+ vtkIdType const numCells = mesh.GetNumberOfCells();
+
+ vtkNew< vtkIdList > indices3D;
+ vtkNew< vtkIdList > indices2D;
+
+ // Classify cells by dimension
+ for( vtkIdType i = 0; i < numCells; ++i )
+ {
+ unsigned char const cellType = mesh.GetCellType( i );
+ int const cellDim = vtkCellTypes::GetDimension( cellType );
+ if( cellDim == 3 )
+ {
+ indices3D->InsertNextId( i );
+ }
+ else if( cellDim == 2 )
+ {
+ indices2D->InsertNextId( i );
+ }
+ // Note: 1D edges and 0D vertices are intentionally ignored
+ }
+
+ GEOS_LOG_RANK_0( GEOS_FMT( "Separating mesh: {} 3D cells, {} 2D cells (from {} total)",
+ indices3D->GetNumberOfIds(),
+ indices2D->GetNumberOfIds(),
+ numCells ) );
+
+ // Extract and deep-copy 3D cells
+ vtkNew< vtkExtractCells > extractor3D;
+ extractor3D->SetInputDataObject( &mesh );
+ extractor3D->SetExtractAllCells( false );
+ extractor3D->SetCellList( indices3D );
+ extractor3D->Update();
+
+ cells3D = vtkSmartPointer< vtkUnstructuredGrid >::New();
+ cells3D->DeepCopy( extractor3D->GetOutput() );
+
+ // Store index mapping for 3D cells
+ cells3DToOriginal.resize( indices3D->GetNumberOfIds() );
+ for( vtkIdType i = 0; i < indices3D->GetNumberOfIds(); ++i )
+ {
+ cells3DToOriginal[i] = indices3D->GetId( i );
+ }
+
+ // Extract and deep-copy 2D cells
+ vtkNew< vtkExtractCells > extractor2D;
+ extractor2D->SetInputDataObject( &mesh );
+ extractor2D->SetExtractAllCells( false );
+ extractor2D->SetCellList( indices2D );
+ extractor2D->Update();
+
+ cells2D = vtkSmartPointer< vtkUnstructuredGrid >::New();
+ cells2D->DeepCopy( extractor2D->GetOutput() );
+
+ // Store index mapping for 2D cells
+ cells2DToOriginal.resize( indices2D->GetNumberOfIds() );
+ for( vtkIdType i = 0; i < indices2D->GetNumberOfIds(); ++i )
+ {
+ cells2DToOriginal[i] = indices2D->GetId( i );
+ }
+}
+
+/**
+ * @brief Build mapping from 2D cells to their neighboring 3D cells using global cell IDs
+ * @param[in] mesh Original mesh containing both 2D and 3D cells with global IDs
+ * @param[in] cells2DToOriginal Mapping from local 2D indices to original mesh indices
+ * @param[in] cells3DToOriginal Mapping from local 3D indices to original mesh indices
+ * @return ArrayOfArrays mapping each 2D cell index to global IDs of neighboring 3D cells
+ */
+static ArrayOfArrays< vtkIdType, int64_t >
+build2DTo3DNeighbors( vtkDataSet & mesh,
+ arrayView1d< vtkIdType const > cells2DToOriginal,
+ arrayView1d< vtkIdType const > cells3DToOriginal )
+{
+ GEOS_MARK_FUNCTION;
+
+ // Retrieve global cell ID array
+ vtkDataArray * globalCellIds = mesh.GetCellData()->GetGlobalIds();
+ GEOS_ERROR_IF( globalCellIds == nullptr,
+ "Global cell IDs must be present in mesh for 2D-3D neighbor mapping" );
+
+
+// Build reverse lookup: original mesh index -> global cell ID (3D cells only)
+ stdUnorderedMap< vtkIdType, int64_t > original3DToGlobalId;
+ original3DToGlobalId.reserve( cells3DToOriginal.size() );
+
+ for( localIndex i = 0; i < cells3DToOriginal.size(); ++i )
+ {
+ vtkIdType const origIdx = cells3DToOriginal[i];
+ int64_t const globalId = static_cast< int64_t >( globalCellIds->GetTuple1( origIdx ) );
+ original3DToGlobalId.emplace( origIdx, globalId );
+ }
+
+ ArrayOfArrays< vtkIdType, int64_t > neighbors2Dto3D;
+ neighbors2Dto3D.reserve( cells2DToOriginal.size() );
+
+ // Topology statistics
+ localIndex numStandalone = 0;
+ localIndex numBoundary = 0; // 1 neighbor
+ localIndex numInternal = 0; // 2 neighbors
+ localIndex numJunction = 0; // >2 neighbors
+
+ vtkNew< vtkIdList > neighborCells;
+ vtkNew< vtkIdList > pointIds2D;
+
+ // Build neighbor list for each 2D cell
+ for( localIndex i = 0; i < cells2DToOriginal.size(); ++i )
+ {
+ vtkIdType const origIdx2D = cells2DToOriginal[i];
+ mesh.GetCellPoints( origIdx2D, pointIds2D );
+
+ // Find all cells sharing ALL nodes with this 2D cell (exact face match)
+ neighborCells->Reset();
+ mesh.GetCellNeighbors( origIdx2D, pointIds2D, neighborCells );
+
+ // Filter for 3D neighbors and retrieve their global IDs
+ array1d< int64_t > neighbor3DGlobalIds;
+ neighbor3DGlobalIds.reserve( neighborCells->GetNumberOfIds() );
+
+ for( vtkIdType n = 0; n < neighborCells->GetNumberOfIds(); ++n )
+ {
+ vtkIdType const neighborIdx = neighborCells->GetId( n );
+
+ // Only consider 3D cells
+ if( mesh.GetCell( neighborIdx )->GetCellDimension() == 3 )
+ {
+ auto it = original3DToGlobalId.find( neighborIdx );
+ GEOS_ERROR_IF( it == original3DToGlobalId.end(),
+ GEOS_FMT( "3D neighbor at index {} not found in 3D cell mapping", neighborIdx ) );
+
+ neighbor3DGlobalIds.emplace_back( it->second );
+ }
+ }
+
+ // Update topology statistics
+ localIndex const numNeighbors = neighbor3DGlobalIds.size();
+ switch( numNeighbors )
+ {
+ case 0: numStandalone++; break;
+ case 1: numBoundary++; break;
+ case 2: numInternal++; break;
+ default: numJunction++; break;
+ }
+
+ neighbors2Dto3D.appendArray( neighbor3DGlobalIds.begin(), neighbor3DGlobalIds.end() );
+ }
+
+ // Print diagnostic summary
+ GEOS_LOG_RANK_0( "\n2D-to-3D Neighbor Topology " );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Total 2D cells: {}", cells2DToOriginal.size() ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Standalone (0 neighbors): {}", numStandalone ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Boundary (1 neighbor): {}", numBoundary ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Internal (2 neighbors): {}", numInternal ) );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Junction (>2 neighbors):{}", numJunction ) );
+
+ // Standalone 2D cells indicate mesh topology errors
+ GEOS_ERROR_IF( numStandalone > 0,
+ GEOS_FMT( "{} orphaned 2D cells detected with no 3D neighbors. "
+ "These may be artifacts or detached surfaces. "
+ "Please clean the mesh or verify the geometry.", numStandalone ) );
+
+ return neighbors2Dto3D;
+}
+
+/**
+ * @brief Assign partition ownership for 2D cells based on their 3D neighbor locations
+ *
+ * Uses a deterministic tie-breaking rule (minimum global ID) to assign 2D cells that
+ * have multiple 3D neighbors (e.g., internal fractures).
+ *
+ * @param[in] neighbors2Dto3D Neighbor map from build2DTo3DNeighbors()
+ * @param[in] globalIdTo3DPartition Complete map of 3D cell global ID -> rank assignment
+ * @param[in] comm MPI communicator (unused, kept for API consistency)
+ * @return Partition (rank) assignment for each 2D cell
+ */
+static array1d< int64_t >
+assign2DCellsTo3DPartitions( ArrayOfArrays< vtkIdType, int64_t > const & neighbors2Dto3D,
+ stdUnorderedMap< int64_t, int64_t > const & globalIdTo3DPartition )
+{
+ GEOS_MARK_FUNCTION;
+
+ array1d< int64_t > partitions2D( neighbors2Dto3D.size() );
+
+ for( localIndex i = 0; i < neighbors2Dto3D.size(); ++i )
+ {
+ auto neighbors = neighbors2Dto3D[i];
+
+ GEOS_ASSERT_MSG( neighbors.size() > 0,
+ "2D cell should have at least one neighbor (enforced by build2DTo3DNeighbors)" );
+
+ // Deterministic tie-breaking: choose neighbor with minimum global ID
+ int64_t minGlobalId = neighbors[0];
+ for( localIndex n = 1; n < neighbors.size(); ++n )
+ {
+ if( neighbors[n] < minGlobalId )
+ {
+ minGlobalId = neighbors[n];
+ }
+ }
+
+ // Look up partition assignment from the complete 3D partition map
+ auto it = globalIdTo3DPartition.find( minGlobalId );
+ GEOS_ERROR_IF( it == globalIdTo3DPartition.end(),
+ GEOS_FMT( "3D neighbor with global ID {} not found in partition map", minGlobalId ) );
+
+ partitions2D[i] = it->second;
+ }
+
+ return partitions2D;
+}
+
+AllMeshes
+merge2D3DCellsAndRedistribute( vtkSmartPointer< vtkDataSet > redistributed3D,
+ vtkSmartPointer< vtkUnstructuredGrid > cells2D,
+ ArrayOfArrays< vtkIdType, int64_t > const & neighbors2Dto3D,
+ stdMap< string, vtkSmartPointer< vtkDataSet > > const & unpartitionedFractures,
+ stdMap< string, ArrayOfArrays< vtkIdType, int64_t > > const & fractureNeighbors,
+ stdVector< string > const & fractureNames,
+ MPI_Comm const comm )
+{
+ GEOS_MARK_FUNCTION;
+
+ int const rank = MpiWrapper::commRank( comm );
+ int const numRanks = MpiWrapper::commSize( comm );
+
+ // -----------------------------------------------------------------------
+ // Step 1: Build complete 3D partition map (all ranks participate)
+ // -----------------------------------------------------------------------
+ vtkDataArray * redistributed3DGlobalIds = redistributed3D->GetCellData()->GetGlobalIds();
+ GEOS_ERROR_IF( redistributed3DGlobalIds == nullptr,
+ "Global cell IDs required in redistributed 3D mesh" );
+
+ vtkIdType const local3DCells = redistributed3D->GetNumberOfCells();
+
+ // Gather counts from all ranks
+ array1d< int > cellCounts;
+ MpiWrapper::allGather( static_cast< int >( local3DCells ), cellCounts, comm );
+
+ // Compute offsets
+ array1d< int > offsets( numRanks );
+ int totalCells = 0;
+ for( int r = 0; r < numRanks; ++r )
+ {
+ offsets[r] = totalCells;
+ totalCells += cellCounts[r];
+ }
+
+ // Gather local global IDs
+ array1d< int64_t > localGlobalIds( local3DCells );
+ for( vtkIdType i = 0; i < local3DCells; ++i )
+ {
+ localGlobalIds[i] = static_cast< int64_t >( redistributed3DGlobalIds->GetTuple1( i ) );
+ }
+
+ // All-gather global IDs to all ranks
+ array1d< int64_t > allGlobalIds( totalCells );
+ MpiWrapper::allgatherv( localGlobalIds.data(), static_cast< int >( local3DCells ),
+ allGlobalIds.data(), cellCounts.data(), offsets.data(),
+ comm );
+
+ // Build complete partition map: global cell ID -> owning rank
+ stdUnorderedMap< int64_t, int64_t > complete3DPartitionMap;
+ complete3DPartitionMap.reserve( totalCells );
+
+ for( int r = 0; r < numRanks; ++r )
+ {
+ for( int i = 0; i < cellCounts[r]; ++i )
+ {
+ int64_t const globalId = allGlobalIds[offsets[r] + i];
+ complete3DPartitionMap.emplace( globalId, r );
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 2: Rank 0 assigns and redistributes 2D cells
+ // -----------------------------------------------------------------------
+ vtkSmartPointer< vtkPartitionedDataSet > split2DCells = vtkSmartPointer< vtkPartitionedDataSet >::New();
+ vtkIdType expected2DCells = 0;
+
+ if( rank == 0 )
+ {
+ vtkIdType const numCells2D = cells2D->GetNumberOfCells();
+ expected2DCells = numCells2D;
+
+ array1d< int64_t > partitions2D = assign2DCellsTo3DPartitions(
+ neighbors2Dto3D,
+ complete3DPartitionMap );
+
+ split2DCells = splitMeshByPartition( cells2D, numRanks, partitions2D.toViewConst() );
+ }
+ else
+ {
+ split2DCells->SetNumberOfPartitions( numRanks );
+ for( int r = 0; r < numRanks; ++r )
+ {
+ vtkNew< vtkUnstructuredGrid > emptyPart;
+ split2DCells->SetPartition( r, emptyPart );
+ }
+ }
+
+ vtkSmartPointer< vtkUnstructuredGrid > local2DCells = vtk::redistribute( *split2DCells, comm );
+
+ // Conservation check
+ vtkIdType const total2DCells = MpiWrapper::sum( local2DCells->GetNumberOfCells(), comm );
+ MpiWrapper::broadcast( expected2DCells, 0, comm );
+
+ GEOS_ERROR_IF( total2DCells != expected2DCells,
+ GEOS_FMT( "2D cell redistribution failed: expected {} cells, got {} cells",
+ expected2DCells, total2DCells ) );
+
+ // -----------------------------------------------------------------------
+ // Step 3: Redistribute fracture elements
+ // -----------------------------------------------------------------------
+ stdMap< string, vtkSmartPointer< vtkDataSet > > redistributedFractures;
+
+ for( string const & fractureName : fractureNames )
+ {
+ vtkSmartPointer< vtkPartitionedDataSet > splitFracture =
+ vtkSmartPointer< vtkPartitionedDataSet >::New();
+
+ vtkIdType expectedFractureCells = 0;
+
+ if( rank == 0 )
+ {
+ auto fracIt = unpartitionedFractures.find( fractureName );
+ vtkSmartPointer< vtkDataSet > unpartitionedFracture =
+ (fracIt != unpartitionedFractures.end()) ? fracIt->second : nullptr;
+
+ if( !unpartitionedFracture || unpartitionedFracture->GetNumberOfCells() == 0 )
+ {
+ // Empty fracture - create empty partitions
+ splitFracture->SetNumberOfPartitions( numRanks );
+ for( int r = 0; r < numRanks; ++r )
+ {
+ vtkNew< vtkUnstructuredGrid > emptyPart;
+ splitFracture->SetPartition( r, emptyPart );
+ }
+ }
+ else
+ {
+ vtkIdType const numFractureCells = unpartitionedFracture->GetNumberOfCells();
+ expectedFractureCells = numFractureCells;
+
+ auto fracNeighborIt = fractureNeighbors.find( fractureName );
+ GEOS_ERROR_IF( fracNeighborIt == fractureNeighbors.end(),
+ "Fracture '" << fractureName << "' not found in fractureNeighbors map" );
+
+ ArrayOfArrays< vtkIdType, int64_t > const & neighbors = fracNeighborIt->second;
+
+ GEOS_ERROR_IF( neighbors.size() != numFractureCells,
+ "Fracture '" << fractureName << "' has " << numFractureCells
+ << " cells but " << neighbors.size() << " neighbor entries" );
+
+ // Ensure GlobalIds exist on fracture mesh
+ vtkDataArray * fracGlobalIds = unpartitionedFracture->GetCellData()->GetGlobalIds();
+ GEOS_ERROR_IF( !fracGlobalIds,
+ "Fracture '" << fractureName << "' missing GlobalIds - should have been set by manageGlobalIds()" );
+
+ // Assign fractures to ranks using DETERMINISTIC min-neighbor strategy
+ array1d< int64_t > partitionsFracture( numFractureCells );
+
+ for( localIndex i = 0; i < numFractureCells; ++i )
+ {
+ auto neighbors3D = neighbors[i];
+
+ if( neighbors3D.size() > 0 )
+ {
+ // Find neighbor with MINIMUM global ID (deterministic)
+ int64_t minNeighborId = neighbors3D[0];
+ for( localIndex n = 1; n < neighbors3D.size(); ++n )
+ {
+ if( neighbors3D[n] < minNeighborId )
+ minNeighborId = neighbors3D[n];
+ }
+
+ // Assign fracture to rank owning that neighbor
+ auto it = complete3DPartitionMap.find( minNeighborId );
+ GEOS_ERROR_IF( it == complete3DPartitionMap.end(),
+ "Fracture '" << fractureName << "' element " << i
+ << " neighbor " << minNeighborId << " not in partition map" );
+
+ partitionsFracture[i] = it->second;
+ }
+ else
+ {
+ GEOS_ERROR( "Fracture '" << fractureName << "' element " << i << " has no 3D neighbors" );
+ }
+ }
+
+ splitFracture = splitMeshByPartition( unpartitionedFracture, numRanks,
+ partitionsFracture.toViewConst() );
+ }
+ }
+ else
+ {
+ // Non-root ranks: create empty partitions
+ splitFracture->SetNumberOfPartitions( numRanks );
+ for( int r = 0; r < numRanks; ++r )
+ {
+ vtkNew< vtkUnstructuredGrid > emptyPart;
+ splitFracture->SetPartition( r, emptyPart );
+ }
+ }
+
+ // Redistribute fracture across ranks
+ vtkSmartPointer< vtkUnstructuredGrid > localFracture =
+ vtk::redistribute( *splitFracture, comm );
+
+ // Conservation check
+ vtkIdType const totalFractureCells = MpiWrapper::sum( localFracture->GetNumberOfCells(), comm );
+ MpiWrapper::broadcast( expectedFractureCells, 0, comm );
+ GEOS_ERROR_IF( totalFractureCells != expectedFractureCells,
+ "Fracture '" << fractureName << "' redistribution lost cells: expected "
+ << expectedFractureCells << ", got " << totalFractureCells );
+
+ redistributedFractures.insert( { fractureName, localFracture } );
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 4: Merge local 3D and 2D cells on each rank
+ // -----------------------------------------------------------------------
+ vtkSmartPointer< vtkUnstructuredGrid > mergedMesh = vtkSmartPointer< vtkUnstructuredGrid >::New();
+
+ if( local2DCells->GetNumberOfCells() > 0 )
+ {
+ vtkNew< vtkAppendFilter > appendFilter;
+ appendFilter->AddInputData( redistributed3D );
+ appendFilter->AddInputData( local2DCells );
+ appendFilter->MergePointsOn(); // Ensures shared nodes are not duplicated
+ appendFilter->Update();
+
+ mergedMesh->DeepCopy( appendFilter->GetOutput() );
+ }
+ else
+ {
+ mergedMesh->DeepCopy( redistributed3D );
+ }
+
+ return AllMeshes( mergedMesh, redistributedFractures );
+}
+
+/**
+ * @brief Find 3D cells whose faces exactly match a fracture element
+ *
+ * A 3D cell matches if it has a face that shares all nodes with the fracture element,
+ * accounting for collocated nodes at split interfaces.
+ *
+ * @param fractureNodeIds Local node IDs of the fracture element
+ * @param collocatedNodes Mapping from local node ID to all collocated global IDs
+ * @param nodesToCells Reverse map from global node ID to cells containing that node
+ * @return Global IDs of matching 3D cells (typically 0-2 neighbors for fractures)
+ */
+stdVector< vtkIdType > findMatchingCellsForFractureElement(
+ vtkIdList * fractureNodeIds,
+ CollocatedNodes const & collocatedNodes,
+ stdMap< vtkIdType, std::set< vtkIdType > > const & nodesToCells )
+{
+ vtkIdType const numFractureNodes = fractureNodeIds->GetNumberOfIds();
+
+ // Build set of ALL collocated nodes for this fracture element
+ std::unordered_set< vtkIdType > fractureCollocatedNodes;
+ fractureCollocatedNodes.reserve( numFractureNodes * 2 ); // Typical case: 2 versions per node
+
+ for( vtkIdType j = 0; j < numFractureNodes; ++j )
+ {
+ vtkIdType const localNodeIdx = fractureNodeIds->GetId( j );
+ stdVector< vtkIdType > const & ns = collocatedNodes[ localNodeIdx ];
+ fractureCollocatedNodes.insert( ns.begin(), ns.end() );
+ }
+
+ // Build map: candidate cellId -> set of its nodes that match fracture's collocated nodes
+ std::unordered_map< vtkIdType, std::unordered_set< vtkIdType > > cellToMatchedNodes;
+ cellToMatchedNodes.reserve( fractureCollocatedNodes.size() );
+
+ for( vtkIdType const & collocatedNode : fractureCollocatedNodes )
+ {
+ auto it = nodesToCells.find( collocatedNode );
+ if( it != nodesToCells.cend() )
+ {
+ for( vtkIdType const & cellId : it->second )
+ {
+ cellToMatchedNodes[cellId].insert( collocatedNode );
+ }
+ }
+ }
+
+ // Filter to cells that form a valid matching face
+ stdVector< vtkIdType > matchingCells;
+ matchingCells.reserve( 2 ); // Most fractures have 0-2 neighbors
+
+ for( auto const & [cellId, matchedNodes] : cellToMatchedNodes )
+ {
+ // Must match exactly the number of fracture nodes
+ if( matchedNodes.size() != static_cast< std::size_t >( numFractureNodes ) )
+ {
+ continue;
+ }
+
+ // Verify each fracture node has at least one collocated version in the matched set
+ // (ensures we matched a true face, not just any numFractureNodes nodes)
+ bool allFractureNodesRepresented = true;
+
+ for( vtkIdType j = 0; j < numFractureNodes; ++j )
+ {
+ vtkIdType const localNodeIdx = fractureNodeIds->GetId( j );
+ stdVector< vtkIdType > const & nodeCollocated = collocatedNodes[ localNodeIdx ];
+
+ // Check if ANY collocated version of this node is in the matched set
+ bool nodeRepresented = std::any_of(
+ nodeCollocated.begin(),
+ nodeCollocated.end(),
+ [&matchedNodes]( vtkIdType collocNode ) {
+ return matchedNodes.count( collocNode ) > 0;
+ }
+ );
+
+ if( !nodeRepresented )
+ {
+ allFractureNodesRepresented = false;
+ break;
+ }
+ }
+
+ if( allFractureNodesRepresented )
+ {
+ matchingCells.push_back( cellId );
+
+ // Early exit - most fractures have ≤2 neighbors (boundary or internal)
+ if( matchingCells.size() >= 2 )
+ {
+ return matchingCells;
+ }
+ }
+ }
+
+ return matchingCells;
+}
+
+/**
+ * @brief Build fracture-to-3D neighbor connectivity using collocated nodes
+ *
+ * @param fractureMesh The fracture mesh (2D elements)
+ * @param originalMesh The original 3D mesh containing both geometries
+ * @param cells3DToOriginal Mapping from local 3D cell index to original mesh index
+ * @return ArrayOfArrays mapping each fracture element to global IDs of neighboring 3D cells
+ */
+ArrayOfArrays< vtkIdType, int64_t >
+buildFractureTo3DNeighbors( vtkSmartPointer< vtkDataSet > fractureMesh,
+ vtkSmartPointer< vtkDataSet > originalMesh,
+ arrayView1d< vtkIdType const > cells3DToOriginal )
+{
+ GEOS_MARK_FUNCTION;
+
+ vtkIdType const numFractureElems = fractureMesh->GetNumberOfCells();
+ vtkDataArray * globalCellIds = originalMesh->GetCellData()->GetGlobalIds();
+ vtkDataArray * globalNodeIds = originalMesh->GetPointData()->GetGlobalIds();
+
+ GEOS_ERROR_IF( globalCellIds == nullptr, "Original mesh must have GlobalIds for cells" );
+ GEOS_ERROR_IF( globalNodeIds == nullptr, "Original mesh must have GlobalIds for points" );
+
+ // Build mapping: original mesh index -> global cell ID (3D cells only)
+ std::unordered_map< vtkIdType, int64_t > original3DToGlobalId;
+ original3DToGlobalId.reserve( cells3DToOriginal.size() );
+
+ for( localIndex i = 0; i < cells3DToOriginal.size(); ++i )
+ {
+ vtkIdType const origIdx = cells3DToOriginal[i];
+ int64_t const globalId = static_cast< int64_t >( globalCellIds->GetTuple1( origIdx ) );
+ original3DToGlobalId.emplace( origIdx, globalId );
+ }
+
+ // Build node-to-cell connectivity for 3D cells
+ stdMap< vtkIdType, std::set< vtkIdType > > nodeToOriginalCells;
+
+ for( localIndex i = 0; i < cells3DToOriginal.size(); ++i )
+ {
+ vtkIdType const origCellIdx = cells3DToOriginal[i];
+ vtkCell * cell = originalMesh->GetCell( origCellIdx );
+ vtkIdList * pointIds = cell->GetPointIds();
+
+ for( vtkIdType p = 0; p < pointIds->GetNumberOfIds(); ++p )
+ {
+ vtkIdType const nodeLocalId = pointIds->GetId( p );
+ vtkIdType const nodeGlobalId = static_cast< vtkIdType >( globalNodeIds->GetTuple1( nodeLocalId ) );
+ nodeToOriginalCells.get_inserted( nodeGlobalId ).insert( origCellIdx );
+ }
+ }
+
+ // Build collocated nodes wrapper for fracture mesh
+ CollocatedNodes collocatedNodesObj( "fracture", fractureMesh, false ); // Skip MPI check for local operation
+
+ // Build fracture-to-3D neighbor mapping
+ ArrayOfArrays< vtkIdType, int64_t > result;
+ result.reserve( numFractureElems );
+
+ for( vtkIdType fracElemId = 0; fracElemId < numFractureElems; ++fracElemId )
+ {
+ vtkCell * fracCell = fractureMesh->GetCell( fracElemId );
+ vtkIdList * fracPointIds = fracCell->GetPointIds();
+
+ // Find matching 3D cells using exact node matching
+ stdVector< vtkIdType > matchingOriginalCells = findMatchingCellsForFractureElement(
+ fracPointIds,
+ collocatedNodesObj,
+ nodeToOriginalCells );
+
+ // Convert to global IDs
+ array1d< int64_t > neighbor3DGlobalIds;
+ neighbor3DGlobalIds.reserve( matchingOriginalCells.size() );
+
+ for( vtkIdType origCellIdx : matchingOriginalCells )
+ {
+ auto it = original3DToGlobalId.find( origCellIdx );
+ if( it != original3DToGlobalId.end() )
+ {
+ neighbor3DGlobalIds.emplace_back( it->second );
+ }
+ }
+
+ // Error on orphaned fractures
+ GEOS_ERROR_IF( neighbor3DGlobalIds.empty(),
+ GEOS_FMT( "Fracture element {} has no 3D neighbors", fracElemId ) );
+
+ result.appendArray( neighbor3DGlobalIds.begin(), neighbor3DGlobalIds.end() );
+ }
+
+ return result;
+}
+
vtkSmartPointer< vtkUnstructuredGrid >
@@ -919,7 +1742,7 @@ redistributeByAreaGraphAndLayer( AllMeshes & input,
if( haveLayer0 )
{
AllMeshes layer0input( layer0, {} ); // fracture mesh not supported yet
- layer0Parts = partitionByCellGraph( layer0input, method, subComm, numPartA, 3, numRefinements );
+ layer0Parts = partitionByCellGraph( layer0input.getMainMesh(), method, subComm, numPartA, 3, numRefinements );
MpiWrapper::commFree( subComm );
}
@@ -1195,7 +2018,33 @@ ensureNoEmptyRank( vtkSmartPointer< vtkDataSet > mesh,
return vtk::redistribute( *splitMesh, MPI_COMM_GEOS );
}
-
+/**
+ * @brief Redistribute mesh across MPI ranks ensuring 2D-3D co-location
+ *
+ * This function implements a multi-stage redistribution workflow:
+ * 1. Separates 2D (surfaces/fractures) and 3D (volumes) cells
+ * 2. Builds neighbor connectivity (2D->3D, fracture->3D)
+ * 3. Tags super-cells (groups that must stay together) if fractures present
+ * 4. Redistributes 3D mesh (with super-cell constraints if present, else standard/structured)
+ * 5. Assigns 2D/fractures to ranks based on 3D neighbor ownership
+ * 6. Merges all components on each rank
+ *
+ * @param logLevel Logging verbosity level
+ * @param loadedMesh Input mesh (may contain mixed 2D/3D cells)
+ * @param namesToFractures Map of fracture name -> fracture mesh (must be on rank 0)
+ * @param comm MPI communicator
+ * @param method Partitioning algorithm (parmetis/ptscotch)
+ * @param partitionRefinement Number of refinement iterations
+ * @param useGlobalIds Whether to generate/use global IDs
+ * @param structuredIndexAttributeName Attribute for structured mesh layers (optional, mutually exclusive with super-cells)
+ * @param numPartZ Number of partitions in Z direction for structured meshes
+ * @return Redistributed AllMeshes with main mesh + fractures
+ *
+ * @pre Fractures must be on rank 0
+ * @post 2D elements co-located with their 3D neighbors
+ * @post Super-cells (if present) remain intact on single ranks
+ * @post Structured meshes partitioned by layers (if structuredIndexAttributeName provided and no super-cells)
+ */
AllMeshes
redistributeMeshes( integer const logLevel,
vtkSmartPointer< vtkDataSet > loadedMesh,
@@ -1203,12 +2052,16 @@ redistributeMeshes( integer const logLevel,
MPI_Comm const comm,
PartitionMethod const method,
int const partitionRefinement,
+ int const partitionFractureWeight,
int const useGlobalIds,
string const & structuredIndexAttributeName,
int const numPartZ )
{
GEOS_MARK_FUNCTION;
+ int const rank = MpiWrapper::commRank( comm );
+ int const numRanks = MpiWrapper::commSize( comm );
+
stdVector< vtkSmartPointer< vtkDataSet > > fractures;
for( auto & nameToFracture: namesToFractures )
{
@@ -1218,69 +2071,294 @@ redistributeMeshes( integer const logLevel,
// Generate global IDs for vertices and cells, if needed
vtkSmartPointer< vtkDataSet > mesh = manageGlobalIds( loadedMesh, useGlobalIds, !std::empty( fractures ) );
- if( MpiWrapper::commRank( comm ) != ( MpiWrapper::commSize( comm ) - 1 ) )
+ // Verify fractures are on rank 0
+ if( rank != 0 )
{
- for( auto nameToFracture: namesToFractures )
+ for( auto const & [name, fracture]: namesToFractures )
{
- GEOS_ASSERT_EQ( nameToFracture.second->GetNumberOfCells(), 0 );
+ GEOS_UNUSED_VAR( name );
+ GEOS_ASSERT_EQ( fracture->GetNumberOfCells(), 0 );
}
}
- // Determine if redistribution is required
- vtkIdType const minCellsOnAnyRank = MpiWrapper::min( mesh->GetNumberOfCells(), comm );
- if( minCellsOnAnyRank == 0 )
+ // -----------------------------------------------------------------------
+ // Step 1: Separate 2D and 3D cells from main mesh
+ // -----------------------------------------------------------------------
+ vtkSmartPointer< vtkUnstructuredGrid > cells3D, cells2D;
+ array1d< vtkIdType > cells3DToOriginal, cells2DToOriginal;
+ separateCellsByDimension( *mesh, cells3D, cells2D, cells3DToOriginal, cells2DToOriginal );
+
+ // Build 2D-to-3D neighbor mapping using global IDs
+ ArrayOfArrays< vtkIdType, int64_t > neighbors2Dto3D;
+ if( cells2D->GetNumberOfCells() > 0 )
+ {
+ neighbors2Dto3D = build2DTo3DNeighbors( *mesh,
+ cells2DToOriginal.toViewConst(),
+ cells3DToOriginal.toViewConst() );
+ }
+ else
+ {
+ neighbors2Dto3D.resize( 0, 0 );
+ }
+
+ GEOS_LOG_RANK_0( GEOS_FMT( "Separated main mesh into {} 3D cells and {} 2D cells",
+ cells3D->GetNumberOfCells(), cells2D->GetNumberOfCells() ));
+
+ // -----------------------------------------------------------------------
+ // Step 1b: Build fracture-to-3D neighbor mapping (rank 0 only)
+ // -----------------------------------------------------------------------
+ stdMap< string, ArrayOfArrays< vtkIdType, int64_t > > fractureNeighbors;
+ stdVector< string > fractureNames;
+
+ if( rank == 0 )
+ {
+ for( auto const & [fractureName, fractureMesh]: namesToFractures )
+ {
+ fractureNames.push_back( fractureName );
+ }
+ GEOS_LOG_RANK_0( GEOS_FMT( "Building fracture neighbors for {} fractures", fractureNames.size() ));
+ }
+
+ // Broadcast fracture names to all ranks
+ {
+ int numFractures = LvArray::integerConversion< int >( fractureNames.size() );
+ MpiWrapper::broadcast( numFractures, 0, comm );
+
+ if( rank != 0 )
+ {
+ fractureNames.resize( numFractures );
+ }
+
+ for( string & name : fractureNames )
+ {
+ MpiWrapper::broadcast( name, 0, comm );
+ }
+ }
+
+ // Initialize empty neighbor arrays on ALL ranks
+ for( auto const & fractureName : fractureNames )
+ {
+ fractureNeighbors.get_inserted( fractureName ).resize( 0, 0 );
+ }
+
+ // Build actual neighbors only on rank 0
+ if( rank == 0 )
{
- // Redistribute the mesh over all ranks using simple octree partitions
- mesh = redistributeByKdTree( *mesh );
+ for( auto & [fractureName, fractureMesh]: namesToFractures )
+ {
+ if( fractureMesh && fractureMesh->GetNumberOfCells() > 0 )
+ {
+ GEOS_LOG_RANK_0( GEOS_FMT( "Building neighbors for fracture '{}' with {} elements",
+ fractureName, fractureMesh->GetNumberOfCells() ));
+
+ fractureNeighbors[fractureName] = buildFractureTo3DNeighbors(
+ fractureMesh,
+ mesh,
+ cells3DToOriginal.toViewConst()
+ );
+
+ GEOS_LOG_RANK_0( GEOS_FMT( "Finished building {} neighbors for '{}'",
+ fractureNeighbors[fractureName].size(), fractureName ));
+ }
+ }
+ mesh = nullptr; // Free the original mesh
}
- // Check if a rank does not have a cell after the redistribution
- // If this is the case, we need a fix otherwise the next redistribution will fail
- // We expect this function to only be called in some pathological cases
- if( MpiWrapper::min( mesh->GetNumberOfCells(), comm ) == 0 )
+ MpiWrapper::barrier( comm );
+
+ // -----------------------------------------------------------------------
+ // Step 2: Tag 3D cells with super-cell IDs (rank 0 only, if fractures exist)
+ // -----------------------------------------------------------------------
+ SuperCellInfo superCellInfo;
+ bool hasSuperCells = false;
+
+ if( rank == 0 && !fractureNames.empty() )
{
- mesh = ensureNoEmptyRank( mesh, comm );
+ superCellInfo = tagCellsWithSuperCellIds( cells3D, fractureNeighbors, partitionFractureWeight );
+
+ vtkIdTypeArray * scArray =
+ vtkIdTypeArray::SafeDownCast( cells3D->GetCellData()->GetArray( "SuperCellId" ) );
+ hasSuperCells = (scArray != nullptr);
+
+ if( hasSuperCells )
+ {
+ GEOS_LOG_RANK_0( GEOS_FMT( "Tagged {} super-cells from fracture connectivity",
+ superCellInfo.atomicSuperCells.size() ));
+ }
}
- AllMeshes result;
- // Redistribute the mesh again using higher-quality graph partitioner
- if( !structuredIndexAttributeName.empty() )
+ // Broadcast whether we have super-cells
{
- AllMeshes input( mesh, namesToFractures );
- result = redistributeByAreaGraphAndLayer( input,
- method,
- structuredIndexAttributeName,
- comm,
- numPartZ,
- partitionRefinement - 1 );
+ int hasSuperCellsInt = hasSuperCells ? 1 : 0;
+ MpiWrapper::broadcast( hasSuperCellsInt, 0, comm );
+ hasSuperCells = (hasSuperCellsInt > 0);
}
- else if( partitionRefinement > 0 )
+
+ // -----------------------------------------------------------------------
+ // Step 3: Redistribute the 3D cells
+ // -----------------------------------------------------------------------
+ vtkIdType const minCellsOnAnyRank = MpiWrapper::min( cells3D->GetNumberOfCells(), comm );
+
+ vtkSmartPointer< vtkDataSet > redistributed3D;
+
+ // Initial redistribution if some ranks are empty
+ if( minCellsOnAnyRank == 0 )
{
- AllMeshes input( mesh, namesToFractures );
- result = redistributeByCellGraph( input, method, comm, partitionRefinement - 1 );
+ if( hasSuperCells )
+ {
+ GEOS_LOG_RANK_0( "Initial redistribution preserving super-cells..." );
+ redistributed3D = redistributeBySuperCellBlocks( cells3D, comm );
+ }
+ else
+ {
+ GEOS_LOG_RANK_0( "Initial redistribution using KD-tree..." );
+ redistributed3D = redistributeByKdTree( *cells3D );
+
+ // Safety check for pathological cases
+ if( MpiWrapper::min( redistributed3D->GetNumberOfCells(), comm ) == 0 )
+ {
+ redistributed3D = ensureNoEmptyRank( redistributed3D, comm );
+ }
+ }
}
else
{
- result.setMainMesh( mesh );
- result.setFaceBlocks( namesToFractures );
+ redistributed3D = cells3D;
}
+ // Free the original 3D mesh.
+ cells3D = nullptr;
+
+ // Check all ranks have cells after initial redistribution
+ GEOS_ERROR_IF( redistributed3D->GetNumberOfCells() == 0, "Rank has no cells after initial redistribution." );
- // Logging some information about the redistribution.
+
+ // Refine partitioning based on mesh characteristics
+ if( partitionRefinement > 0 )
{
- string const pattern = "{}: {}";
- stdVector< string > messages;
- messages.push_back( GEOS_FMT( pattern, "Local mesh size", result.getMainMesh()->GetNumberOfCells() ) );
- for( auto const & [faceName, faceMesh]: result.getFaceBlocks() )
+ // Check if we still have super-cells (they might have been added during redistribution)
+ vtkIdTypeArray * scArray =
+ vtkIdTypeArray::SafeDownCast( redistributed3D->GetCellData()->GetArray( "SuperCellId" ) );
+
+ int localHasSuperCells = (scArray != nullptr) ? 1 : 0;
+ int globalHasSuperCells = MpiWrapper::max( localHasSuperCells, comm );
+
+ if( globalHasSuperCells > 0 )
{
- messages.push_back( GEOS_FMT( pattern, faceName, faceMesh->GetNumberOfCells() ) );
+ // Use super-cell aware partitioning (fractures present - keep super-cells intact)
+ GEOS_LOG_RANK_0( "Refining partition with super-cell constraints..." );
+
+ redistributed3D = redistributeBySuperCellGraph(
+ redistributed3D,
+ method,
+ comm,
+ partitionRefinement - 1,
+ partitionFractureWeight
+ );
}
- if( logLevel >= 5 )
+ else if( !structuredIndexAttributeName.empty() )
+ {
+ // Use structured mesh layered partitioning (no fractures/super-cells)
+ GEOS_LOG_RANK_0( "Refining partition with structured mesh layers..." );
+
+ // Wrap in AllMeshes for compatibility with redistributeByAreaGraphAndLayer
+ AllMeshes tempWrapper;
+ tempWrapper.setMainMesh( redistributed3D );
+ tempWrapper.setFaceBlocks( stdMap< string, vtkSmartPointer< vtkDataSet > >() ); // Empty fractures
+
+ tempWrapper = redistributeByAreaGraphAndLayer(
+ tempWrapper,
+ method,
+ structuredIndexAttributeName,
+ comm,
+ numPartZ,
+ partitionRefinement - 1 );
+
+ redistributed3D = tempWrapper.getMainMesh();
+ }
+ else
{
- GEOS_LOG_RANK( stringutilities::join( messages, ", " ) );
+ // Use standard graph partitioning (no fractures, no structured mesh)
+ GEOS_LOG_RANK_0( "Refining partition with standard cell graph..." );
+
+ redistributed3D = redistributeByCellGraph(
+ redistributed3D,
+ method,
+ comm,
+ partitionRefinement - 1
+ );
}
+
+ // Check all ranks have cells after refined redistribution
+ GEOS_ERROR_IF( redistributed3D->GetNumberOfCells() == 0, "Rank has no cells after refined redistribution." );
+
}
- return result;
+ // -----------------------------------------------------------------------
+ // Step 4: Merge 2D cells AND fractures back with redistributed 3D cells
+ // -----------------------------------------------------------------------
+ AllMeshes finalResult = merge2D3DCellsAndRedistribute(
+ redistributed3D,
+ cells2D,
+ neighbors2Dto3D,
+ namesToFractures,
+ fractureNeighbors,
+ fractureNames,
+ comm );
+
+ // -----------------------------------------------------------------------
+ // Step 5: Final logging
+ // -----------------------------------------------------------------------
+ if( logLevel >=5 )
+ {
+ vtkIdType local2DCells = 0;
+ vtkIdType local3DCells = 0;
+ vtkIdType localFractureCells = 0;
+
+ vtkSmartPointer< vtkDataSet > finalMesh = finalResult.getMainMesh();
+ for( vtkIdType i = 0; i < finalMesh->GetNumberOfCells(); ++i )
+ {
+ int const dim = finalMesh->GetCell( i )->GetCellDimension();
+ if( dim == 2 )
+ local2DCells++;
+ else if( dim == 3 )
+ local3DCells++;
+ }
+
+ for( auto const & [faceName, faceMesh] : finalResult.getFaceBlocks() )
+ {
+ localFractureCells += faceMesh->GetNumberOfCells();
+ }
+
+ array1d< vtkIdType > all2D, all3D, allFracture;
+ MpiWrapper::allGather( local2DCells, all2D, comm );
+ MpiWrapper::allGather( local3DCells, all3D, comm );
+ MpiWrapper::allGather( localFractureCells, allFracture, comm );
+
+ if( rank == 0 )
+ {
+ GEOS_LOG_RANK_0( "\n----------------------------------------------------------" );
+ GEOS_LOG_RANK_0( "| Rk | 3D Cells | 2D Cells | Fractures | Total |" );
+ GEOS_LOG_RANK_0( "----------------------------------------------------------" );
+ vtkIdType sum2D = 0, sum3D = 0, sumFracture = 0;
+ for( int r = 0; r < numRanks; ++r )
+ {
+ sum2D += all2D[r];
+ sum3D += all3D[r];
+ sumFracture += allFracture[r];
+
+ GEOS_LOG_RANK_0( GEOS_FMT( "| {:>2} | {:9} | {:9} | {:9} | {:13} |",
+ r, all3D[r], all2D[r], allFracture[r],
+ all3D[r] + all2D[r] + allFracture[r] ) );
+ }
+ GEOS_LOG_RANK_0( "----------------------------------------------------------" );
+ GEOS_LOG_RANK_0( GEOS_FMT( "|Tot | {:9} | {:9} | {:9} | {:13} |",
+ sum3D, sum2D, sumFracture,
+ sum3D + sum2D + sumFracture ) );
+ GEOS_LOG_RANK_0( "----------------------------------------------------------" );
+ }
+ }
+
+ return finalResult;
}
/**
diff --git a/src/coreComponents/mesh/generators/VTKUtilities.hpp b/src/coreComponents/mesh/generators/VTKUtilities.hpp
index ef660d7641c..23ebb5dcb23 100644
--- a/src/coreComponents/mesh/generators/VTKUtilities.hpp
+++ b/src/coreComponents/mesh/generators/VTKUtilities.hpp
@@ -31,11 +31,14 @@
#include
#include
+
namespace geos
{
namespace vtk
{
+class CollocatedNodes;
+
/**
* @brief Choice of advanced mesh partitioner
*/
@@ -164,6 +167,7 @@ redistributeMeshes( integer const logLevel,
MPI_Comm const comm,
PartitionMethod const method,
int const partitionRefinement,
+ int const partitionFractureWeight,
int const useGlobalIds,
string const & structuredIndexAttributeName,
int const numPartZ );
@@ -243,6 +247,22 @@ void importRegularField( stdVector< vtkIdType > const & cellIds,
void importRegularField( vtkDataArray * vtkArray,
dataRepository::WrapperBase & wrapper );
+/**
+ * @brief Find 3D cells whose faces exactly match a fracture element.
+ *
+ * A 3D cell matches if it has a face that shares all nodes with the fracture element,
+ * accounting for collocated nodes at split interfaces.
+ *
+ * @param fractureNodeIds Local node IDs of the fracture element
+ * @param collocatedNodes Mapping from local node ID to all collocated global IDs
+ * @param nodesToCells Reverse map from global node ID to cells containing that node
+ * @return Global IDs of matching 3D cells (typically 0-2 neighbors for fractures)
+ */
+stdVector< vtkIdType > findMatchingCellsForFractureElement(
+ vtkIdList * fractureNodeIds,
+ CollocatedNodes const & collocatedNodes,
+ stdMap< vtkIdType, std::set< vtkIdType > > const & nodesToCells );
+
} // namespace vtk