fix: Fracture/3D cell co-location in parallel mesh redistribution

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(
     <Mesh>
       <VTKMesh
         name="mesh"
         file="{}"
         partitionRefinement="0"
-        useGlobalIds="0"
+        useGlobalIds="{}"
         {} />
     </Mesh>
   )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,14 +167,26 @@ 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 },
         { 0, 0, 1 },
         { 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,15 +195,17 @@ class TestFractureImport : public ::testing::Test
       }
       main->SetPoints( points );
 
-      int constexpr numHexs = 2;
-      vtkIdType const cubes[numHexs][8] = {
-        { 0, 1, 2, 3, 4, 5, 6, 7 },
-        { 8, 9, 10, 11, 12, 13, 14, 15 }
+      int constexpr numHexs = 3;
+      int constexpr pointsPerHex = 8;
+      vtkIdType const cubes[numHexs][pointsPerHex] = {
+        { 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 )
       {
-        main->InsertNextCell( VTK_HEXAHEDRON, 8, cube );
+        main->InsertNextCell( VTK_HEXAHEDRON, pointsPerHex, cube );
       }
 
       vtkNew< vtkIdTypeArray > cellGlobalIds;
@@ -207,15 +227,17 @@ 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;
       double const pointsCoords[numPoints][3] = {
         { 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 )
@@ -225,11 +247,12 @@ class TestFractureImport : public ::testing::Test
       fracture->SetPoints( points );
 
       int constexpr numQuads = 1;
-      vtkIdType const quad[numQuads][4] = { { 0, 1, 2, 3 } };
+      int constexpr pointsPerQuad = 4;
+      vtkIdType const quad[numQuads][pointsPerQuad] = { { 0, 1, 2, 3 } };
       fracture->Allocate( numQuads );
       for( vtkIdType const * q: quad )
       {
-        fracture->InsertNextCell( VTK_QUAD, numPoints, q );
+        fracture->InsertNextCell( VTK_QUAD, pointsPerQuad, q );
       }
 
       vtkNew< vtkIdTypeArray > cellGlobalIds;
@@ -250,15 +273,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 +312,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 +357,6 @@ TEST_F( TestFractureImport, fracture )
   TestMeshImport( m_vtkFile, validate, "fracture" );
 }
 
-
 TEST( VTKImport, cube )
 {
   auto validate = []( CellBlockManagerABC const & cellBlockManager ) -> void

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 )

src/coreComponents/mesh/generators/CollocatedNodes.cpp

@@ -24,29 +24,35 @@ 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 )
+  // With performGlobalCheck, a NULL on the local rank is OK if any other rank found the field
+  // (vtk does not provide fields for empty meshes). Otherwise we require it on every rank.
+  // After MPI max a non-zero result means the field exists on at
+  // least one rank, so the field is not globally missing.
+  bool const isMissing = performGlobalCheck
+                         ? ( MpiWrapper::max( reinterpret_cast< std::uintptr_t >( collocatedNodes ) ) == 0 )
+                         : ( collocatedNodes == nullptr );
+
+  if( isMissing )
   {
-    GEOS_LOG_RANK_0( "Available point data fields in '" << faceBlockName << "':" );
+    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( " - " << faceMesh->GetPointData()->GetArrayName( i ) << " of type '" << faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() << "'" );
+      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 \"{}\".",
+    GEOS_ERROR( GEOS_FMT( "Could not find valid field \"{}\" for fracture \"{}\"{}.",
                           COLLOCATED_NODES,
-                          faceBlockName ) );
+                          faceBlockName,
+                          performGlobalCheck ? "" : " on this rank" ) );
   }
 
   if( collocatedNodes )

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.

src/coreComponents/mesh/generators/ParMETISInterface.cpp

@@ -97,6 +97,8 @@ partition( ArrayOfArraysView< idx_t const, idx_t > const & graph,
            MPI_Comm comm,
            int const numRefinements )
 {
+  GEOS_ERROR_IF( numParts <= 0, "Number of partitions must be strictly positive" );
+
   array1d< idx_t > part( graph.size() ); // all 0 by default
   if( numParts == 1 )
   {
@@ -147,5 +149,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 )
+{
+  GEOS_ERROR_IF( numParts <= 0, "Number of partitions must be strictly positive" );
+
+  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]
+  idx_t options[4] = { 1, 0, 2022, PARMETIS_PSR_UNCOUPLED };
+
+  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

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