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 72b8e218102..3c6c403b2e6 100644
--- a/src/coreComponents/mesh/generators/CollocatedNodes.cpp
+++ b/src/coreComponents/mesh/generators/CollocatedNodes.cpp
@@ -24,29 +24,51 @@ 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( GEOS_FMT( "Available point data fields in '{}':", faceBlockName ) );
+ for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )
+ {
+ GEOS_LOG_RANK_0( GEOS_FMT( " - {} of type '{}'",
+ faceMesh->GetPointData()->GetArrayName( i ),
+ faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() ) );
+ }
+ GEOS_ERROR( GEOS_FMT( "Could not find valid field \"{}\" for fracture \"{}\".",
+ COLLOCATED_NODES,
+ faceBlockName ) );
+ }
+ }
+ else
+ {
+ if( !collocatedNodes )
+ {
+ GEOS_LOG_RANK_0( GEOS_FMT( "Available point data fields in '{}':", faceBlockName ) );
+ for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )
+ {
+ GEOS_LOG_RANK_0( GEOS_FMT( " - {} of type '{}'",
+ faceMesh->GetPointData()->GetArrayName( i ),
+ faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() ) );
+ }
+ GEOS_ERROR( GEOS_FMT( "Could not find valid field \"{}\" for fracture \"{}\" on this rank.",
+ COLLOCATED_NODES,
+ faceBlockName ) );
}
- GEOS_ERROR( GEOS_FMT( "Could not find valid field \"{}\" for fracture \"{}\".",
- COLLOCATED_NODES,
- 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 d177313f44c..e68f7eca1ea 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 8e948f762e1..b4cc93b5087 100644
--- a/src/coreComponents/mesh/generators/ParMETISInterface.hpp
+++ b/src/coreComponents/mesh/generators/ParMETISInterface.hpp
@@ -65,6 +65,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..21f0d17fafe 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 b90bd1651e7..4a881414782 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..4302819e3e9
--- /dev/null
+++ b/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.cpp
@@ -0,0 +1,1052 @@
+/**
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 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
+ }
+}
+
+SuperCellInfo tagCellsWithSuperCellIds(
+ vtkSmartPointer< vtkUnstructuredGrid > cells3D,
+ stdMap< string, ArrayOfArrays< localIndex, int64_t > > const & fractureNeighbors,
+
+ integer fractureWeight )
+{
+ 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 totalFractureElements = 0;
+
+ for( auto const & [fractureName, neighbors] : fractureNeighbors )
+ {
+ totalFractureElements += static_cast< vtkIdType >( neighbors.size() );
+
+ for( vtkIdType i = 0; i < neighbors.size(); ++i )
+ {
+ auto neighborList = neighbors[i];
+
+ GEOS_ERROR_IF( neighborList.size() != 2,
+ GEOS_FMT( "Fracture '{}' element {} has {} 3D neighbors (expected exactly 2). ",
+ fractureName, i, neighborList.size() ) );
+
+ vtkIdType gidA = neighborList[0];
+ vtkIdType gidB = neighborList[1];
+
+ fractureGraph[gidA].insert( gidB );
+ fractureGraph[gidB].insert( gidA );
+ }
+ }
+
+ GEOS_LOG_RANK_0( GEOS_FMT( "Built fracture graph: {} 3D cells connected via {} fracture elements",
+ fractureGraph.size(), totalFractureElements ) );
+
+ // -----------------------------------------------------------------------
+ // 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;
+ visited.insert( 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 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 )
+ {
+ // Mark as atomic super-cell (cannot be split during partitioning)
+ info.atomicSuperCells.insert( scId );
+ }
+
+ 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
+ // -----------------------------------------------------------------------
+ GEOS_LOG_RANK_0( GEOS_FMT( "Super-cells: {} created (largest: {} cells)",
+ numSuperCellsCreated, 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 );
+ }
+
+ 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 const 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 ({} ordering)",
+ 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 )
+ {
+ vtkPoints * meshPoints = cells3D->GetPoints();
+
+ // Compute super-cell centroids directly
+ for( auto const & [scId, cellIndices] : superCellToLocalCells )
+ {
+ std::array< double, 3 > centroid = {0.0, 0.0, 0.0};
+ vtkIdType totalPoints = 0;
+
+ // Accumulate all points from all cells in this super-cell
+ for( vtkIdType cellIdx : cellIndices )
+ {
+ vtkIdType npts;
+ const vtkIdType * pts;
+ cells3D->GetCellPoints( cellIdx, npts, pts );
+
+ 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];
+ }
+ totalPoints += npts;
+ }
+
+ // Average over all points
+ centroid[0] /= totalPoints;
+ centroid[1] /= totalPoints;
+ centroid[2] /= totalPoints;
+
+ superCells.push_back( SuperCellDistributionInfo{ scId, cellIndices, centroid } );
+ }
+
+ // 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] );
+ }
+ }
+
+ // 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 );
+ } );
+
+ }
+ else
+ {
+ // BLOCK: Simple ordering by super-cell ID (no centroid computation needed)
+ 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;
+ } );
+ }
+
+ // -----------------------------------------------------------------------
+ // 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 global index distribution across ranks
+ 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
+ array1d< pmet_idx_t > localParmetisIndices;
+ array1d< vtkIdType > localVtkIds;
+ array1d< vtkIdType > localSuperCellIds;
+
+ localParmetisIndices.reserve( numLocalCells );
+ localVtkIds.reserve( numLocalCells );
+ localSuperCellIds.reserve( numLocalCells );
+
+ for( vtkIdType i = 0; i < numLocalCells; ++i )
+ {
+ localParmetisIndices.emplace_back( myParmetisStart + i );
+ localVtkIds.emplace_back( cellGlobalIds->GetValue( i ) );
+ localSuperCellIds.emplace_back( superCellIdArray->GetValue( i ) );
+ }
+
+ // Gather all mappings
+ array1d< pmet_idx_t > allParmetisIndices;
+ array1d< vtkIdType > allVtkIds;
+ array1d< vtkIdType > allSuperCellIds;
+
+ MpiWrapper::allGatherv( localParmetisIndices.toViewConst(), allParmetisIndices, comm );
+ MpiWrapper::allGatherv( localVtkIds.toViewConst(), allVtkIds, comm );
+ MpiWrapper::allGatherv( localSuperCellIds.toViewConst(), allSuperCellIds, 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( localIndex i = 0; i < allParmetisIndices.size(); ++i )
+ {
+ pmet_idx_t const parmetisIdx = allParmetisIndices[i];
+
+ // Only store if it's local or a ghost we need
+ if( (parmetisIdx >= myParmetisStart && parmetisIdx < myParmetisEnd) ||
+ ghostParmetisSet.count( parmetisIdx ) > 0 )
+ {
+ vtkIdType const vtkId = allVtkIds[i];
+ vtkIdType const superCellId = allSuperCellIds[i];
+
+ parmetisToVtk[parmetisIdx] = vtkId;
+ vtkToSuperCell[vtkId] = superCellId;
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 4: Exchange super-cell global indices
+ // -----------------------------------------------------------------------
+ array1d< vtkIdType > sendSCIds;
+ array1d< pmet_idx_t > sendSCGlobalIndices;
+
+ sendSCIds.reserve( superCellIdToGlobalIdx.size() );
+ sendSCGlobalIndices.reserve( superCellIdToGlobalIdx.size() );
+
+ for( auto const & [scId, gIdx] : superCellIdToGlobalIdx )
+ {
+ sendSCIds.emplace_back( scId );
+ sendSCGlobalIndices.emplace_back( gIdx );
+ }
+
+ array1d< vtkIdType > allSCIds;
+ array1d< pmet_idx_t > allSCGlobalIndices;
+
+ MpiWrapper::allGatherv( sendSCIds.toViewConst(), allSCIds, comm );
+ MpiWrapper::allGatherv( sendSCGlobalIndices.toViewConst(), allSCGlobalIndices, comm );
+
+ // Update local map with all super-cell mappings
+ for( localIndex i = 0; i < allSCIds.size(); ++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 );
+ }
+
+ // -----------------------------------------------------------------------
+ // 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() );
+ }
+
+ 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 );
+
+ 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 for local edges
+ // -----------------------------------------------------------------------
+ {
+ 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 );
+
+ // -----------------------------------------------------------------------
+ // 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 )
+ {
+ 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 ) );
+
+ 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..65d895fec6b
--- /dev/null
+++ b/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.hpp
@@ -0,0 +1,167 @@
+/**
+ * ------------------------------------------------------------------------------------------------------------
+ * 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/DataTypes.hpp"
+#include "common/MpiWrapper.hpp"
+#include "common/TimingMacros.hpp"
+#include "mesh/generators/VTKMeshGeneratorTools.hpp"
+#include "mesh/generators/ParMETISInterface.hpp"
+#include "LvArray/src/ArrayOfArrays.hpp"
+
+#include
+
+class vtkUnstructuredGrid;
+class vtkDataSet;
+
+namespace geos
+{
+namespace vtk
+{
+
+/**
+ * @brief Initial distribution strategy for super-cells
+ */
+enum class InitialDistributionStrategy
+{
+ MORTON, ///< Morton Z-curve ordering for spatial locality
+ BLOCK ///< Contiguous block distribution
+};
+
+/**
+ * @brief Super-cell metadata for constrained partitioning
+ *
+ * Groups cells that must remain co-located (e.g., fracture-connected cells).
+ */
+struct SuperCellInfo
+{
+ /// SuperCellId to global cell IDs
+ std::map< vtkIdType, std::vector< vtkIdType > > superCellToOriginalCells;
+
+ /// SuperCellId to vertex weight
+ std::map< vtkIdType, vtkIdType > vertexWeights;
+
+ /// SuperCellIds containing multiple cells
+ std::set< vtkIdType > atomicSuperCells;
+};
+
+/**
+ * @brief Assign super-cell IDs based on fracture connectivity
+ *
+ * Cells sharing a fracture get the same ID; isolated cells use their global ID.
+ * Adds a "SuperCellId" array to the mesh. Runs on rank 0.
+ *
+ * @param cells3D Volumetric mesh (modified in-place)
+ * @param fractureNeighbors Fracture element to 3D neighbor mappings
+ * @param fractureWeight Load balancing weight boost for fracture super-cells
+ * @return Super-cell metadata
+ */
+SuperCellInfo tagCellsWithSuperCellIds(
+ vtkSmartPointer< vtkUnstructuredGrid > cells3D,
+ stdMap< string, ArrayOfArrays< localIndex, int64_t > > const & fractureNeighbors,
+ integer fractureWeight );
+
+/**
+ * @brief Rebuild super-cell metadata from SuperCellId array
+ *
+ * After redistribution, each rank reconstructs local super-cell info.
+ * Applies fracture weights to account for contact computation cost.
+ *
+ * @param mesh Distributed mesh with SuperCellId array
+ * @param fractureWeight Weight boost for multi-cell super-cells
+ * @return Local super-cell metadata
+ */
+SuperCellInfo reconstructSuperCellInfo( vtkSmartPointer< vtkUnstructuredGrid > mesh,
+ integer fractureWeight );
+
+/**
+ * @brief Initial redistribution preserving super-cell integrity
+ *
+ * Distributes super-cells across ranks without graph partitioning.
+ * Faster than ParMETIS for initial scatter.
+ *
+ * @param cells3D Input mesh (non-empty on rank 0 only)
+ * @param comm MPI communicator
+ * @param strategy Distribution strategy (MORTON or BLOCK)
+ * @return Redistributed mesh with preserved SuperCellId array
+ */
+vtkSmartPointer< vtkDataSet >
+redistributeBySuperCellBlocks( vtkSmartPointer< vtkUnstructuredGrid > cells3D,
+ MPI_Comm comm,
+ InitialDistributionStrategy strategy = InitialDistributionStrategy::MORTON );
+
+/**
+ * @brief Build super-cell adjacency graph
+ *
+ * Collapses the cell graph into a super-cell graph where each vertex
+ * represents a super-cell and edges connect neighboring super-cells.
+ *
+ * @param cells3D Mesh with SuperCellId array
+ * @param baseGraph Original cell-to-cell adjacency
+ * @param baseElemDist Element distribution for base graph
+ * @param info Super-cell metadata
+ * @param comm MPI communicator
+ * @return (super-cell graph, 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 before partitioning
+ *
+ * Checks for self-loops, invalid indices, duplicate edges,
+ * isolated vertices, and invalid weights.
+ *
+ * @param superGraph Super-cell adjacency graph
+ * @param superElemDist Element distribution array
+ * @param vertexWeights Vertex weights
+ * @param comm MPI communicator
+ */
+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 Map super-cell partitions back to individual cells
+ *
+ * Converts super-cell -> rank assignments into cell -> rank assignments.
+ * All cells in a super-cell go to the same rank.
+ *
+ * @param cells3D Mesh with SuperCellId array
+ * @param superPartitioning Super-cell to rank assignments from ParMETIS
+ * @param superCellIdToLocalIdx SuperCellId to local index mapping
+ * @param comm MPI communicator
+ * @return Cell-level partition assignments
+ */
+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 360bc74ad4d..463e2ca05fd 100644
--- a/src/coreComponents/mesh/generators/VTKUtilities.cpp
+++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp
@@ -22,6 +22,7 @@
#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"
@@ -237,59 +238,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;
@@ -301,46 +281,30 @@ 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 );
}
+
/**
* @brief Split a mesh by partitioning it
*
@@ -593,13 +557,14 @@ 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 );
@@ -618,7 +583,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,
@@ -627,16 +592,10 @@ 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.
- // 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 3D cells only. Fractures are assigned later based on 3D neighbors.
array1d< pmet_idx_t > const elemDist( numRanks + 1 );
{
array1d< pmet_idx_t > elemCounts;
@@ -644,21 +603,8 @@ 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;
-
- // 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 );
+ // Build element-to-node connectivity for the 3D mesh only
+ ArrayOfArrays< pmet_idx_t, pmet_idx_t > const elemToNodes = buildElemToNodes< pmet_idx_t >( mesh3D );
ArrayOfArrays< pmet_idx_t, pmet_idx_t > graph;
#ifdef GEOS_USE_PARMETIS
graph = parmetis::meshToDual( elemToNodes.toViewConst(), elemDist, comm, minCommonNodes );
@@ -698,57 +644,260 @@ 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 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
*/
-AllMeshes
-redistributeByCellGraph( AllMeshes & input,
+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 );
+}
+
+
+/**
+ * @brief Partition mesh while keeping super-cells (cell groups) intact
+ *
+ * Builds a weighted graph where vertices represent super-cells, then partitions
+ * using ParMETIS/PTScotch. Ensures cells with the same SuperCellId stay together.
+ *
+ * @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
+ */
+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, "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
+ );
+
+ // -----------------------------------------------------------------------
+ // Step 4: Compute super-cell element distribution
+ // -----------------------------------------------------------------------
+ array1d< pmet_idx_t > superElemDist( numRanks + 1 );
+ {
+ 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 );
+ }
+
+ // -----------------------------------------------------------------------
+ // 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
+ );
+ }
+ else
+ {
+ GEOS_ERROR( GEOS_FMT( "Unsupported partition method: {}",
+ EnumStrings< PartitionMethod >::toString( method ) ) );
+ }
- // 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 7: Map SuperCellId to local super-cell indices for partition unpacking
+ // -----------------------------------------------------------------------
+
+ 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 )
{
- vtkSmartPointer< vtkPartitionedDataSet > const splitFracMesh = splitMeshByPartition( fracture, numRanks, newFracturePartitions[fractureName].toViewConst() );
- vtkSmartPointer< vtkUnstructuredGrid > const finalFracMesh = vtk::redistribute( *splitFracMesh, MPI_COMM_GEOS );
- finalFractures.insert( {fractureName, finalFracMesh} );
+ orderedSuperCellIds.push_back( scId );
}
- return AllMeshes( finalMesh, finalFractures );
+ // 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( GEOS_FMT( "Final 3D cells redistribution: {} total, {} avg/rank, [{}-{}] range, {:.1f}% imbalance",
+ totalCells, static_cast< vtkIdType >(avgCells),
+ minCells, maxCells, imbalance ) );
+ }
+
+ return redistributed;
}
@@ -980,12 +1129,22 @@ build2DTo3DNeighbors( vtkDataSet & mesh,
}
// Print diagnostic summary
- GEOS_LOG_RANK_0( "\n2D-to-3D Neighbor Topology" );
- GEOS_LOG_RANK_0( GEOS_FMT( " Total 2D cells: {}", cells2DIndices.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 ) );
+ if( numStandalone > 0 || numJunction > 0 )
+ {
+ // Show detailed breakdown when anomalies exist
+ GEOS_LOG_RANK_0( "\n2D-to-3D Neighbor Topology" );
+ GEOS_LOG_RANK_0( GEOS_FMT( " Total 2D cells: {}", cells2DIndices.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):{}\n", numJunction ) );
+ }
+ else
+ {
+ // Condensed output for normal cases
+ GEOS_LOG_RANK_0( GEOS_FMT( "2D cells: {} (1-neighbor/boundary: {}, 2-neighbor/internal: {})",
+ cells2DIndices.size(), numBoundary, numInternal ) );
+ }
// Standalone 2D cells indicate mesh topology errors
GEOS_ERROR_IF( numStandalone > 0,
@@ -997,20 +1156,23 @@ build2DTo3DNeighbors( vtkDataSet & mesh,
}
/**
- * @brief Assign 2D cells to partitions based on their 3D neighbor locations
+ * @brief Assign cells to partitions based on their 3D neighbor locations
*
+ * Generic function for assigning 2D cells or fracture elements to partitions.
+ * Uses deterministic tie-breaking (minimum neighbor global ID) to ensure
+ * each cell is co-located with at least one of its 3D neighbors.
*
- * @param[in] neighbors2Dto3D Neighbor mapping (2D cell index to 3D global IDs)
+ * @param[in] neighborsTo3D Neighbor mapping (cell index -> 3D global IDs)
* @param[in] local3DGlobalIds Global IDs of 3D cells on this rank
* @param[in] local3DPartitions Partition assignments for local 3D cells
* @param[in] comm MPI communicator
- * @return Partition assignments for each 2D cell
+ * @return Partition assignments for each input cell
*/
static array1d< int >
-assign2DCellsTo3DPartitions( ArrayOfArrays< localIndex, int64_t > const & neighbors2Dto3D,
- arrayView1d< int64_t const > local3DGlobalIds,
- arrayView1d< int const > local3DPartitions,
- MPI_Comm const comm )
+assignCellsBasedOn3DNeighbors( ArrayOfArrays< localIndex, int64_t > const & neighbors2Dto3D,
+ arrayView1d< int64_t const > local3DGlobalIds,
+ arrayView1d< int const > local3DPartitions,
+ MPI_Comm const comm )
{
GEOS_MARK_FUNCTION;
@@ -1159,28 +1321,40 @@ extractGlobalIds( vtkDataSet & mesh )
return result;
}
+
/**
- * @brief Redistribute 2D cells and merge with already-redistributed 3D cells
+ * @brief Redistribute 2D cells and fractures, then merge with already-redistributed 3D cells
*
- * This function completes the 2D/3D redistribution workflow:
+ * This function completes the 2D/3D/fracture redistribution workflow:
* 1. Assigns 2D cells to partitions based on their 3D neighbor locations
* 2. Redistributes 2D cells to appropriate ranks
- * 3. Merges local 2D and 3D cells on each rank into a unified mesh
+ * 3. Assigns fracture cells to partitions based on their 3D neighbor locations
+ * 4. Redistributes fracture cells to appropriate ranks
+ * 5. Merges local 2D and 3D cells on each rank into a unified mesh
+ *
+ * Both 2D cells and fractures are assigned to ranks using the same strategy:
+ * - Find all 3D neighbors of each 2D/fracture element
+ * - Assign to the rank owning the 3D neighbor with minimum global ID (deterministic)
+ * - Uses efficient MPI communication (only gathers needed 3D partition info)
*
* @param[in] redistributed3D Already partitioned 3D cells with global IDs
* @param[in] originalMesh Original mesh containing 2D cells (on rank 0)
* @param[in] cells2DIndices Indices of 2D cells in original mesh
- * @param[in] neighbors2Dto3D Pre-computed 2D-to-3D neighbor mapping
- * @param[in] redistributedFractures Already partitioned fracture meshes (pass-through)
+ * @param[in] neighbors2Dto3D Pre-computed 2D-to-3D neighbor mapping (fracture element -> 3D cell global IDs)
+ * @param[in] unpartitionedFractures Unpartitioned fracture meshes (on rank 0, empty on others)
+ * @param[in] fractureNeighbors Pre-computed fracture-to-3D neighbor mappings (on rank 0, empty on others)
+ * @param[in] fractureNames Ordered list of fracture names (consistent across all ranks)
* @param[in] comm MPI communicator
- * @return Complete AllMeshes object with merged main mesh and fractures
+ * @return Complete AllMeshes object with merged main mesh (3D+2D) and redistributed fractures
*/
static AllMeshes
redistribute2DAndMergeWith3D( vtkSmartPointer< vtkDataSet > redistributed3D,
vtkSmartPointer< vtkDataSet > originalMesh,
arrayView1d< vtkIdType const > cells2DIndices,
ArrayOfArrays< localIndex, int64_t > const & neighbors2Dto3D,
- stdMap< string, vtkSmartPointer< vtkDataSet > > const & redistributedFractures,
+ stdMap< string, vtkSmartPointer< vtkDataSet > > const & unpartitionedFractures,
+ stdMap< string, ArrayOfArrays< localIndex, int64_t > > const & fractureNeighbors,
+ stdVector< string > const & fractureNames,
MPI_Comm const comm )
{
GEOS_MARK_FUNCTION;
@@ -1188,35 +1362,32 @@ redistribute2DAndMergeWith3D( vtkSmartPointer< vtkDataSet > redistributed3D,
int const rank = MpiWrapper::commRank( comm );
int const numRanks = MpiWrapper::commSize( comm );
-
- // Step 1: Assign 2D cells to partitions based on 3D neighbor locations
+ // -----------------------------------------------------------------------
+ // Step 1: Assign and redistribute 2D cells
+ // -----------------------------------------------------------------------
array1d< int64_t > local3DGlobalIds = extractGlobalIds( *redistributed3D );
array1d< int > partitions3D( redistributed3D->GetNumberOfCells() );
partitions3D.setValues< parallelHostPolicy >( rank );
bool const hasLocal2DCells = (rank == 0 && !cells2DIndices.empty());
- // Collective 2D partition assignment (pass empty array if no local 2D cells)
- array1d< int > partitions2D = assign2DCellsTo3DPartitions(
+ array1d< int > partitions2D = assignCellsBasedOn3DNeighbors(
hasLocal2DCells ? neighbors2Dto3D : ArrayOfArrays< localIndex, int64_t >{},
local3DGlobalIds.toViewConst(),
partitions3D.toViewConst(),
comm );
- // Extract 2D cells where applicable
vtkSmartPointer< vtkUnstructuredGrid > cells2D = hasLocal2DCells
? extractCellsByIndices( *originalMesh, cells2DIndices )
: vtkSmartPointer< vtkUnstructuredGrid >::New();
- // Step 2: Split and redistribute 2D cells
- // All ranks participate (cells2D and partitions2D are empty on non-root ranks)
vtkSmartPointer< vtkPartitionedDataSet > split2D =
splitMeshByPartition( cells2D, numRanks, partitions2D.toViewConst() );
vtkSmartPointer< vtkUnstructuredGrid > redistributed2D =
vtk::redistribute( *split2D, comm );
- // Conservation check - verify no 2D cells were lost
+ // Conservation check
vtkIdType const total2DCells = MpiWrapper::sum( redistributed2D->GetNumberOfCells(), comm );
vtkIdType expected2DCells = cells2DIndices.size();
MpiWrapper::broadcast( expected2DCells, 0, comm );
@@ -1225,29 +1396,294 @@ redistribute2DAndMergeWith3D( vtkSmartPointer< vtkDataSet > redistributed3D,
GEOS_FMT( "2D cell redistribution failed: expected {} cells, got {} cells",
expected2DCells, total2DCells ) );
- // Step 3: Merge local 2D and 3D cells on each rank
+ // -----------------------------------------------------------------------
+ // Step 2: Redistribute fractures
+ // -----------------------------------------------------------------------
+ stdMap< string, vtkSmartPointer< vtkDataSet > > redistributedFractures;
+
+ for( string const & fractureName : fractureNames )
+ {
+ vtkIdType expectedFractureCells = 0;
+
+ ArrayOfArrays< localIndex, int64_t > localFractureNeighbors;
+ vtkSmartPointer< vtkDataSet > unpartitionedFracture;
+
+ if( rank == 0 )
+ {
+ auto fracIt = unpartitionedFractures.find( fractureName );
+ unpartitionedFracture = (fracIt != unpartitionedFractures.end()) ? fracIt->second : nullptr;
+
+ if( unpartitionedFracture && unpartitionedFracture->GetNumberOfCells() > 0 )
+ {
+ expectedFractureCells = unpartitionedFracture->GetNumberOfCells();
+
+ auto fracNeighborIt = fractureNeighbors.find( fractureName );
+ GEOS_ERROR_IF( fracNeighborIt == fractureNeighbors.end(),
+ GEOS_FMT( "Fracture '{}' not found in fractureNeighbors map", fractureName ) );
+
+ localFractureNeighbors = fracNeighborIt->second;
+ }
+ }
+
+ array1d< int > partitionsFracture = assignCellsBasedOn3DNeighbors(
+ localFractureNeighbors, // Empty on non-root ranks
+ local3DGlobalIds.toViewConst(),
+ partitions3D.toViewConst(),
+ comm );
+
+ // Split and redistribute
+ vtkSmartPointer< vtkPartitionedDataSet > splitFracture;
+
+ if( rank == 0 && unpartitionedFracture && unpartitionedFracture->GetNumberOfCells() > 0 )
+ {
+ splitFracture = splitMeshByPartition( unpartitionedFracture, numRanks,
+ partitionsFracture.toViewConst() );
+ }
+ else
+ {
+ splitFracture = vtkSmartPointer< vtkPartitionedDataSet >::New();
+ splitFracture->SetNumberOfPartitions( numRanks );
+ for( int r = 0; r < numRanks; ++r )
+ {
+ splitFracture->SetPartition( r, vtkSmartPointer< vtkUnstructuredGrid >::New() );
+ }
+ }
+
+ vtkSmartPointer< vtkUnstructuredGrid > localFracture =
+ vtk::redistribute( *splitFracture, comm );
+
+ vtkIdType const totalFractureCells = MpiWrapper::sum( localFracture->GetNumberOfCells(), comm );
+ MpiWrapper::broadcast( expectedFractureCells, 0, comm );
+
+ GEOS_ERROR_IF( totalFractureCells != expectedFractureCells,
+ GEOS_FMT( "Fracture '{}' redistribution lost cells: expected {}, got {}",
+ fractureName, expectedFractureCells, totalFractureCells ) );
+
+ redistributedFractures.insert( { fractureName, localFracture } );
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 3: Merge local 2D and 3D cells
+ // -----------------------------------------------------------------------
vtkSmartPointer< vtkUnstructuredGrid > mergedMesh = vtkSmartPointer< vtkUnstructuredGrid >::New();
if( redistributed2D->GetNumberOfCells() > 0 )
{
- // Merge 3D and 2D cells using VTK append filter
vtkNew< vtkAppendFilter > appendFilter;
appendFilter->AddInputData( redistributed3D );
appendFilter->AddInputData( redistributed2D );
- appendFilter->MergePointsOn(); // Ensures shared nodes are not duplicated
+ appendFilter->MergePointsOn();
appendFilter->Update();
mergedMesh->ShallowCopy( appendFilter->GetOutput() );
}
else
{
- // Only 3D cells on this rank
mergedMesh->ShallowCopy( 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 & cellToNodes : cellToMatchedNodes )
+ {
+ vtkIdType const cellId = cellToNodes.first;
+ auto const & matchedNodes = cellToNodes.second;
+
+ // 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
+ *
+ * @param originalMesh The original mesh containing both 3D cells and fractures
+ * @param fractureMesh The fracture mesh (separate mesh on rank 0)
+ * @param cells3DIndices Indices of 3D cells in the original mesh
+ * @return Mapping: fracture element index -> global IDs of neighboring 3D cells
+ */
+static ArrayOfArrays< localIndex, int64_t >
+buildFractureTo3DNeighbors( vtkDataSet & originalMesh,
+ vtkSmartPointer< vtkDataSet > fractureMesh,
+ arrayView1d< vtkIdType const > cells3DIndices )
+{
+ GEOS_MARK_FUNCTION;
+
+ vtkIdType const numFractureElems = fractureMesh->GetNumberOfCells();
+
+ vtkDataArray * meshGlobalCellIds = originalMesh.GetCellData()->GetGlobalIds();
+ vtkDataArray * meshGlobalNodeIds = originalMesh.GetPointData()->GetGlobalIds();
+
+ GEOS_ERROR_IF( meshGlobalCellIds == nullptr, "Original mesh must have cell GlobalIds" );
+ GEOS_ERROR_IF( meshGlobalNodeIds == nullptr, "Original mesh must have node GlobalIds" );
+
+ // -----------------------------------------------------------------------
+ // Step 1: Build reverse lookup: original mesh index -> global cell ID (3D cells only)
+ // -----------------------------------------------------------------------
+ stdUnorderedMap< vtkIdType, int64_t > mesh3DIdxToGlobalId;
+ mesh3DIdxToGlobalId.reserve( cells3DIndices.size() );
+
+ for( localIndex i = 0; i < cells3DIndices.size(); ++i )
+ {
+ vtkIdType const origIdx = cells3DIndices[i];
+ int64_t const globalId = static_cast< int64_t >( meshGlobalCellIds->GetTuple1( origIdx ) );
+ mesh3DIdxToGlobalId.emplace( origIdx, globalId );
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 2: Build node-to-3D-cell connectivity for the original mesh
+ // -----------------------------------------------------------------------
+ stdMap< vtkIdType, std::set< vtkIdType > > nodeGlobalIdToCells3D;
+
+ for( localIndex i = 0; i < cells3DIndices.size(); ++i )
+ {
+ vtkIdType const origCellIdx = cells3DIndices[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 >( meshGlobalNodeIds->GetTuple1( nodeLocalId ) );
+
+ nodeGlobalIdToCells3D.get_inserted( nodeGlobalId ).insert( origCellIdx );
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 3: Build collocated nodes for fracture mesh
+ // -----------------------------------------------------------------------
+ CollocatedNodes collocatedNodes( "fracture", fractureMesh, false );
+
+ // -----------------------------------------------------------------------
+ // Step 4: Build fracture-to-3D neighbor mapping
+ // -----------------------------------------------------------------------
+ ArrayOfArrays< localIndex, int64_t > result;
+ result.reserve( numFractureElems );
+
+ for( vtkIdType fracElemIdx = 0; fracElemIdx < numFractureElems; ++fracElemIdx )
+ {
+ vtkCell * fracCell = fractureMesh->GetCell( fracElemIdx );
+ vtkIdList * fracPointIds = fracCell->GetPointIds();
+
+ // Find matching 3D cells using exact node matching
+ stdVector< vtkIdType > matchingOriginalCells = findMatchingCellsForFractureElement(
+ fracPointIds,
+ collocatedNodes,
+ nodeGlobalIdToCells3D );
+
+ // Convert original mesh indices to global cell IDs
+ array1d< int64_t > neighbor3DGlobalIds;
+ neighbor3DGlobalIds.reserve( matchingOriginalCells.size() );
+
+ for( vtkIdType origCellIdx : matchingOriginalCells )
+ {
+ auto it = mesh3DIdxToGlobalId.find( origCellIdx );
+ if( it != mesh3DIdxToGlobalId.end() )
+ {
+ neighbor3DGlobalIds.emplace_back( it->second );
+ }
+ }
+
+ // Error on orphaned fractures
+ GEOS_ERROR_IF( neighbor3DGlobalIds.empty(),
+ GEOS_FMT( "Fracture element {} has no 3D neighbors", fracElemIdx ) );
+
+ result.appendArray( neighbor3DGlobalIds.begin(), neighbor3DGlobalIds.end() );
+ }
+
+ return result;
+}
+
vtkSmartPointer< vtkUnstructuredGrid >
threshold( vtkDataSet & mesh,
@@ -1347,7 +1783,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 );
}
@@ -1632,6 +2068,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 )
@@ -1655,11 +2092,25 @@ redistributeMeshes( integer const logLevel,
mesh = vtkSmartPointer< vtkUnstructuredGrid >::New();
}
+ // Verify fractures are on rank 0
+ if( rank != 0 )
+ {
+ for( auto const & [name, fracture]: namesToFractures )
+ {
+ GEOS_UNUSED_VAR( name );
+ GEOS_ASSERT_EQ( fracture->GetNumberOfCells(), 0 );
+ }
+ }
+
+ // -----------------------------------------------------------------------
// Step 1: Classify cells by dimension
+ // -----------------------------------------------------------------------
array1d< vtkIdType > cells3DIndices, cells2DIndices;
classifyCellsByDimension( *mesh, cells3DIndices, cells2DIndices );
+ // -----------------------------------------------------------------------
// Step 2: Build 2D-to-3D neighbor mapping
+ // -----------------------------------------------------------------------
ArrayOfArrays< localIndex, int64_t > neighbors2Dto3D;
if( !cells2DIndices.empty() )
{
@@ -1668,60 +2119,188 @@ redistributeMeshes( integer const logLevel,
cells3DIndices.toViewConst() );
}
- // Step 3: Extract and redistribute 3D cells (+ fractures)
- vtkIdType const minCellsOnAnyRank = MpiWrapper::min( static_cast< vtkIdType >( cells3DIndices.size() ), comm );
+ // -----------------------------------------------------------------------
+ // Step 3: Build fracture-to-3D neighbor mappings
+ // -----------------------------------------------------------------------
+ stdMap< string, ArrayOfArrays< localIndex, int64_t > > fractureNeighbors;
+ stdVector< string > fractureNames;
- AllMeshes result3DAndFractures;
+ // Collect fracture names on rank 0
+ if( rank == 0 )
+ {
+ for( auto const & [fractureName, fractureMesh]: namesToFractures )
+ {
+ fractureNames.push_back( fractureName );
+ }
+ }
- if( minCellsOnAnyRank == 0 )
+ // Broadcast fracture names to all ranks
{
- // Extract 3D cells for initial redistribution
- vtkSmartPointer< vtkUnstructuredGrid > cells3D = extractCellsByIndices( *mesh, cells3DIndices.toViewConst() );
+ int numFractures = LvArray::integerConversion< int >( fractureNames.size() );
+ MpiWrapper::broadcast( numFractures, 0, comm );
- // Redistribute the 3D mesh over all ranks using simple octree partitions
- vtkSmartPointer< vtkDataSet > redistributed3D = redistributeByKdTree( *cells3D );
+ if( rank != 0 )
+ {
+ fractureNames.resize( numFractures );
+ }
- // Check if a rank does not have a cell after the redistribution
- if( MpiWrapper::min( redistributed3D->GetNumberOfCells(), comm ) == 0 )
+ for( string & name : fractureNames )
{
- redistributed3D = ensureNoEmptyRank( redistributed3D, comm );
+ MpiWrapper::broadcast( name, 0, comm );
}
+ }
- result3DAndFractures.setMainMesh( redistributed3D );
- result3DAndFractures.setFaceBlocks( namesToFractures );
+ // Initialize empty neighbor arrays
+ for( auto const & fractureName : fractureNames )
+ {
+ fractureNeighbors.get_inserted( fractureName ).resize( 0, 0 );
}
- else
+
+ // Build actual neighbors only on rank 0
+ if( rank == 0 )
{
- // Extract 3D cells for partitioning
- vtkSmartPointer< vtkUnstructuredGrid > cells3D = extractCellsByIndices( *mesh, cells3DIndices.toViewConst() );
- result3DAndFractures.setMainMesh( cells3D );
- result3DAndFractures.setFaceBlocks( namesToFractures );
+ for( auto const & [fractureName, fractureMesh]: namesToFractures )
+ {
+ if( fractureMesh && fractureMesh->GetNumberOfCells() > 0 )
+ {
+ fractureNeighbors[fractureName] = buildFractureTo3DNeighbors(
+ *mesh,
+ fractureMesh,
+ cells3DIndices.toViewConst() );
+ }
+ }
}
- // Fine-tune redistribution using higher-quality graph partitioner
- if( !structuredIndexAttributeName.empty() )
+ // -----------------------------------------------------------------------
+ // Step 4: Extract cells and tag super-cells
+ // -----------------------------------------------------------------------
+ SuperCellInfo superCellInfo;
+ bool hasSuperCells = false;
+
+ vtkSmartPointer< vtkUnstructuredGrid > cells3D = extractCellsByIndices( *mesh, cells3DIndices );
+
+ if( rank == 0 && !fractureNames.empty() )
{
- result3DAndFractures = redistributeByAreaGraphAndLayer( result3DAndFractures,
- method,
- structuredIndexAttributeName,
- comm,
- numPartZ,
- partitionRefinement - 1 );
+ superCellInfo = tagCellsWithSuperCellIds( cells3D, fractureNeighbors, partitionFractureWeight );
+
+ // Verify SuperCellId array was actually created
+ 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() ));
+ }
}
- else if( partitionRefinement > 0 )
+
+ // Broadcast whether we have super-cells
{
- result3DAndFractures = redistributeByCellGraph( result3DAndFractures, method, comm, partitionRefinement - 1 );
+ int hasSuperCellsInt = hasSuperCells ? 1 : 0;
+ MpiWrapper::broadcast( hasSuperCellsInt, 0, comm );
+ hasSuperCells = (hasSuperCellsInt > 0);
}
- // Step 4: Redistribute 2D cells and merge with 3D
- AllMeshes finalResult = redistribute2DAndMergeWith3D( result3DAndFractures.getMainMesh(),
- mesh,
- cells2DIndices.toViewConst(),
- neighbors2Dto3D,
- result3DAndFractures.getFaceBlocks(),
- comm );
+ // -----------------------------------------------------------------------
+ // Step 5: Initial redistribution of 3D cells
+ // -----------------------------------------------------------------------
+ vtkIdType const minCellsOnAnyRank = MpiWrapper::min( cells3D->GetNumberOfCells(), comm );
- // Step 5: Final logging
+ vtkSmartPointer< vtkDataSet > redistributed3D;
+
+ if( minCellsOnAnyRank == 0 )
+ {
+ if( hasSuperCells )
+ {
+ redistributed3D = redistributeBySuperCellBlocks( cells3D, comm );
+ }
+ else
+ {
+ GEOS_LOG_RANK_0( "Initial redistribution using KD-tree..." );
+ redistributed3D = redistributeByKdTree( *cells3D );
+
+ if( MpiWrapper::min( redistributed3D->GetNumberOfCells(), comm ) == 0 )
+ {
+ redistributed3D = ensureNoEmptyRank( redistributed3D, comm );
+ }
+ }
+ }
+ else
+ {
+ // All ranks already have cells - no redistribution needed
+ redistributed3D = cells3D;
+ }
+
+ // Check all ranks have cells after redistribution
+ GEOS_ERROR_IF( redistributed3D->GetNumberOfCells() == 0, "Rank has no cells after initial redistribution." );
+
+
+ // -----------------------------------------------------------------------
+ // Step 6: Refined partitioning
+ // -----------------------------------------------------------------------
+ if( partitionRefinement > 0 )
+ {
+ if( hasSuperCells )
+ {
+ GEOS_LOG_RANK_0( "Refining partition with super-cell constraints..." );
+
+ redistributed3D = redistributeBySuperCellGraph(
+ redistributed3D,
+ method,
+ comm,
+ partitionRefinement - 1,
+ partitionFractureWeight );
+ }
+ 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
+ {
+ // Use standard cell graph partitioning
+ GEOS_LOG_RANK_0( "Refining partition with standard cell graph..." );
+
+ redistributed3D = redistributeByCellGraph(
+ redistributed3D,
+ method,
+ comm,
+ partitionRefinement - 1 );
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Step 7: Merge 2D cells and fractures back with redistributed 3D cells
+ // -----------------------------------------------------------------------
+ AllMeshes finalResult = redistribute2DAndMergeWith3D(
+ redistributed3D,
+ mesh,
+ cells2DIndices.toViewConst(),
+ neighbors2Dto3D,
+ namesToFractures,
+ fractureNeighbors,
+ fractureNames,
+ comm );
+
+ // -----------------------------------------------------------------------
+ // Step 8: Final logging
+ // -----------------------------------------------------------------------
if( logLevel >= 5 )
{
vtkIdType local2DCells = 0;
@@ -1775,6 +2354,7 @@ redistributeMeshes( integer const logLevel,
return finalResult;
}
+
/**
* @brief Identify the GEOSX type of the polyhedron
*
diff --git a/src/coreComponents/mesh/generators/VTKUtilities.hpp b/src/coreComponents/mesh/generators/VTKUtilities.hpp
index ef660d7641c..a0181fbfd8b 100644
--- a/src/coreComponents/mesh/generators/VTKUtilities.hpp
+++ b/src/coreComponents/mesh/generators/VTKUtilities.hpp
@@ -35,6 +35,7 @@ namespace geos
{
namespace vtk
{
+class CollocatedNodes;
/**
* @brief Choice of advanced mesh partitioner
@@ -152,6 +153,7 @@ findNeighborRanks( stdVector< vtkBoundingBox > boundingBoxes );
* @param[in] comm the MPI communicator
* @param[in] method the partitioning method
* @param[in] partitionRefinement number of graph partitioning refinement cycles
+ * @param[in] partitionFractureWeight additional weight to fracture-connected super-cells during partitioning
* @param[in] useGlobalIds controls whether global id arrays from the vtk input should be used
* @param[in] structuredIndexAttributeName VTK array name for structured index attribute, if present
* @param[in] numPartZ number of MPI partitions in Z direction (only if @p structuredIndexAttributeName is used)
@@ -164,6 +166,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 +246,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