feat: turn the split function into a meshDoctor refineMesh action

.github/workflows/python-package.yml

@@ -140,40 +140,66 @@ jobs:
   check_integration_label:
     runs-on: ubuntu-latest
     needs: [build]
+    permissions:
+      pull-requests: read
     outputs:
       has_geos_integration_label: ${{ steps.set-label.outputs.has_label }}
     steps:
       - name: Check if PR has '${{ env.LABEL_TEST_GEOS_INTEGRATION }}' label
         id: set-label
+        # Fetch labels live from the GitHub REST API rather than reading
+        # github.event.pull_request.labels. The event payload is a snapshot
+        # frozen at the time the workflow was first triggered, so labels
+        # added after that first run are invisible to re-runs.
+        env:
+          GITHUB_TOKEN: ${{ github.token }}
+          REQUIRED_LABEL: ${{ env.LABEL_TEST_GEOS_INTEGRATION }}
         run: |
           echo "Checking for label..."
           LABEL_FOUND=false
-          LABELS='${{ toJson(github.event.pull_request.labels.*.name) }}'
+          PR_JSON=$(curl --fail --silent --show-error \
+            -H "Accept: application/vnd.github+json" \
+            -H "Authorization: Bearer ${GITHUB_TOKEN}" \
+            "https://api.github.com/repos/${{ github.repository }}/pulls/${{ github.event.pull_request.number }}")
+          LABELS=$(echo "${PR_JSON}" | jq -crM '[.labels[].name]')
           echo "PR Labels: $LABELS"
-          if echo "$LABELS" | grep -q "${{ env.LABEL_TEST_GEOS_INTEGRATION }}"; then
+          if echo "$LABELS" | jq -e --arg label "${REQUIRED_LABEL}" 'index($label) != null' > /dev/null; then
             LABEL_FOUND=true
-            echo "Label '${{ env.LABEL_TEST_GEOS_INTEGRATION }}' found"
+            echo "Label '${REQUIRED_LABEL}' found"
           fi
-          echo "has_label=$LABEL_FOUND" >> $GITHUB_OUTPUT
+          echo "has_label=$LABEL_FOUND" >> "$GITHUB_OUTPUT"
 
   check_force_integration_label:
     runs-on: ubuntu-latest
     # needs: [build]
+    permissions:
+      pull-requests: read
     outputs:
       has_geos_integration_force_label: ${{ steps.set-label.outputs.has_label }}
     steps:
       - name: Check if PR has '${{ env.LABEL_FORCE_GEOS_INTEGRATION }}' label
         id: set-label
+        # Fetch labels live from the GitHub REST API rather than reading
+        # github.event.pull_request.labels. The event payload is a snapshot
+        # frozen at the time the workflow was first triggered, so labels
+        # added after that first run are invisible to re-runs.
+        env:
+          GITHUB_TOKEN: ${{ github.token }}
+          REQUIRED_LABEL: ${{ env.LABEL_FORCE_GEOS_INTEGRATION }}
         run: |
           echo "Checking for label..."
           LABEL_FOUND=false
-          LABELS='${{ toJson(github.event.pull_request.labels.*.name) }}'
+          PR_JSON=$(curl --fail --silent --show-error \
+            -H "Accept: application/vnd.github+json" \
+            -H "Authorization: Bearer ${GITHUB_TOKEN}" \
+            "https://api.github.com/repos/${{ github.repository }}/pulls/${{ github.event.pull_request.number }}")
+          LABELS=$(echo "${PR_JSON}" | jq -crM '[.labels[].name]')
           echo "PR Labels: $LABELS"
-          if echo "$LABELS" | grep -q "${{ env.LABEL_FORCE_GEOS_INTEGRATION }}"; then
+          if echo "$LABELS" | jq -e --arg label "${REQUIRED_LABEL}" 'index($label) != null' > /dev/null; then
             LABEL_FOUND=true
-            echo "Label '${{ env.LABEL_FORCE_GEOS_INTEGRATION }}' found"
+            echo "Label '${REQUIRED_LABEL}' found"
           fi
-          echo "has_label=$LABEL_FOUND" >> $GITHUB_OUTPUT
+          echo "has_label=$LABEL_FOUND" >> "$GITHUB_OUTPUT"
 
   # Step 3: Check if GEOS integration is required based on changed files
   check_geos_integration_required:

docs/geos_mesh_docs/stats.rst

@@ -11,3 +11,12 @@ geos.mesh.stats.meshQualityMetricHelpers module
     :members:
     :undoc-members:
     :show-inheritance:
+
+
+geos.mesh.stats.CellTypeCounterEnhanced module
+----------------------------------------------
+
+.. automodule:: geos.mesh.stats.CellTypeCounterEnhanced
+    :members:
+    :undoc-members:
+    :show-inheritance:

docs/geos_mesh_docs/utils.rst

@@ -57,3 +57,12 @@ geos.mesh.utils.pyvistaTools module
    :members:
    :undoc-members:
    :show-inheritance:
+
+
+geos.mesh.utils.SplitMesh module
+-----------------------------------------------
+
+.. automodule:: geos.mesh.utils.SplitMesh
+   :members:
+   :undoc-members:
+   :show-inheritance:

docs/geos_processing_docs/generic_processing_tools.rst

@@ -55,12 +55,3 @@ MergeBlockEnhanced
     :members:
     :undoc-members:
     :show-inheritance:
-
-
-SplitMesh
--------------------------
-
-.. automodule:: geos.processing.generic_processing_tools.SplitMesh
-    :members:
-    :undoc-members:
-    :show-inheritance:

docs/geos_processing_docs/pre_processing.rst

@@ -3,14 +3,6 @@ Pre-processing filters
 
 The `pre_processing` module of `geos-processing` package contains filters to pre-process meshes for GEOS simulation.
 
-CellTypeCounterEnhanced filter
--------------------------------------------------------------
-
-.. automodule:: geos.processing.pre_processing.CellTypeCounterEnhanced
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
 MeshQualityEnhanced filter
 ---------------------------------------------------------
 

docs/mesh-doctor.rst

@@ -134,6 +134,17 @@ The ``generateGlobalIds`` can generate `global ids` for the imported ``vtk`` mes
 .. command-output:: mesh-doctor generateGlobalIds --help
    :shell:
 
+``refineMesh``
+""""""""""""""
+
+The ``refineMesh`` module refines a mesh by splitting each cell into smaller cells of the same type using edge midpoints
+(hexahedron -> 8, tetrahedron -> 8, pyramid -> 6 pyramids + 4 tetrahedra, triangle -> 4, quad -> 4).
+Shared edges reuse the same midpoint, so the refined mesh stays conformal, including across 2D faces coincident with 3D cell faces.
+The refinement can be applied several times with the ``--iterations`` option.
+
+.. command-output:: mesh-doctor refineMesh --help
+   :shell:
+
 ``nonConformal``
 """"""""""""""""
 

geos-processing/src/geos/processing/pre_processing/CellTypeCounterEnhanced.py → geos-mesh/src/geos/mesh/stats/CellTypeCounterEnhanced.py

@@ -20,7 +20,7 @@
 
 .. code-block:: python
 
-    from geos.processing.pre_processing.CellTypeCounterEnhanced import CellTypeCounterEnhanced
+    from geos.mesh.stats.CellTypeCounterEnhanced import CellTypeCounterEnhanced
 
     # Filter inputs
     inputMesh: vtkUnstructuredGrid

geos-processing/src/geos/processing/generic_processing_tools/SplitMesh.py → geos-mesh/src/geos/mesh/utils/SplitMesh.py

@@ -8,26 +8,27 @@
 
 from vtkmodules.vtkCommonCore import vtkPoints, vtkIdTypeArray, vtkDataArray
 from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkCellArray, vtkCellData, vtkCell, vtkCellTypes,
-                                            VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_HEXAHEDRON, VTK_PYRAMID, VTK_WEDGE,
-                                            VTK_POLYHEDRON, VTK_POLYGON )
+                                            VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_HEXAHEDRON, VTK_PYRAMID, VTK_WEDGE )
 from vtkmodules.util.numpy_support import numpy_to_vtk, vtk_to_numpy
 
 from geos.utils.Logger import ( getLogger, Logger, CountVerbosityHandler, isHandlerInLogger, getLoggerHandlerType )
-from geos.processing.pre_processing.CellTypeCounterEnhanced import CellTypeCounterEnhanced
+from geos.mesh.stats.CellTypeCounterEnhanced import CellTypeCounterEnhanced
 from geos.mesh.model.CellTypeCounts import CellTypeCounts
 
 __doc__ = """
 SplitMesh module is a vtk filter that splits cells of a mesh composed of tetrahedra, pyramids, hexahedra, triangles, and quads.
 
-.. Warning:: Current implementation only supports meshes composed of either polygons or polyhedra, not both together.
+Polygons (2D) and polyhedra (3D) may be mixed in the same mesh: edge midpoints are
+shared through an internal cache, so a 2D cell coincident with a 3D cell face splits
+conformally with that face. Wedges are not supported.
 
 Filter input and output types are vtkUnstructuredGrid.
 
 To use the filter:
 
 .. code-block:: python
 
-    from geos.processing.generic_processing_tools.SplitMesh import SplitMesh
+    from geos.mesh.utils.SplitMesh import SplitMesh
 
     # Filter inputs
     inputMesh: vtkUnstructuredGrid
@@ -74,6 +75,9 @@ def __init__( self, inputMesh: vtkUnstructuredGrid, speHandler: bool = False ) -
         self.cellTypes: list[ int ]
         self.speHandler: bool = speHandler
         self.handler: None | logging.Handler = None
+        # Cache of edge midpoints keyed on (min_pt_id, max_pt_id) so that
+        # shared edges across adjacent cells reuse the same point ID.
+        self.m_edgeMidpointCache: dict[ tuple[ int, int ], int ] = {}
 
         # Logger
         self.logger: Logger
@@ -133,20 +137,14 @@ def applyFilter( self: Self ) -> None:
         if counts.getTypeCount( VTK_WEDGE ) != 0:
             raise TypeError( "Input mesh contains wedges that are not currently supported." )
 
-        nbPolygon: int = counts.getTypeCount( VTK_POLYGON )
-        nbPolyhedra: int = counts.getTypeCount( VTK_POLYHEDRON )
-        # Current implementation only supports meshes composed of either polygons or polyhedra
-        if nbPolyhedra * nbPolygon != 0:
-            raise TypeError(
-                "Input mesh is composed of both polygons and polyhedra, but it must contains only one of the two." )
-
         nbTet: int = counts.getTypeCount( VTK_TETRA )  # will divide into 8 tets
         nbPyr: int = counts.getTypeCount( VTK_PYRAMID )  # will divide into 6 pyramids and 4 tets so 10 new cells
         nbHex: int = counts.getTypeCount( VTK_HEXAHEDRON )  # will divide into 8 hexes
         nbTriangles: int = counts.getTypeCount( VTK_TRIANGLE )  # will divide into 4 triangles
         nbQuad: int = counts.getTypeCount( VTK_QUAD )  # will divide into 4 quads
-        nbNewPoints: int = 0
-        nbNewPoints = nbHex * 19 + nbTet * 6 + nbPyr * 9 if nbPolyhedra > 0 else nbTriangles * 3 + nbQuad * 5
+        # Upper bound: shared edges (e.g. between a 3D cell and its 2D face) produce one midpoint
+        # not two, thanks to the edge midpoint cache. The actual count will be <= this.
+        nbNewPoints: int = nbHex * 19 + nbTet * 6 + nbPyr * 9 + nbTriangles * 3 + nbQuad * 5
         nbNewCells: int = nbHex * 8 + nbTet * 8 + nbPyr * 10 + nbTriangles * 4 + nbQuad * 4
 
         self.points = vtkPoints()
@@ -159,6 +157,7 @@ def applyFilter( self: Self ) -> None:
         self.originalId.SetName( "OriginalID" )
         self.originalId.Allocate( nbNewCells )
         self.cellTypes = []
+        self.m_edgeMidpointCache = {}
 
         # Define cell type to splitting method mapping
         splitMethods = {
@@ -233,10 +232,16 @@ def _addMidPoint( self: Self, ptA: int, ptB: int ) -> int:
         Returns:
             int: inserted point Id
         """
+        key: tuple[ int, int ] = ( min( ptA, ptB ), max( ptA, ptB ) )
+        cached = self.m_edgeMidpointCache.get( key )
+        if cached is not None:
+            return cached
         ptACoor: npt.NDArray[ np.float64 ] = np.array( self.points.GetPoint( ptA ) )
         ptBCoor: npt.NDArray[ np.float64 ] = np.array( self.points.GetPoint( ptB ) )
         center: npt.NDArray[ np.float64 ] = ( ptACoor + ptBCoor ) / 2.
-        return self.points.InsertNextPoint( center[ 0 ], center[ 1 ], center[ 2 ] )
+        newId: int = self.points.InsertNextPoint( center[ 0 ], center[ 1 ], center[ 2 ] )
+        self.m_edgeMidpointCache[ key ] = newId
+        return newId
 
     def _splitTetrahedron( self: Self, cell: vtkCell, index: int ) -> None:
         r"""Split a tetrahedron.

geos-mesh/tests/conftest.py

@@ -162,6 +162,10 @@ def _get_dataset( datasetType: str ) -> vtkMultiBlockDataSet | vtkPolyData | vtk
             vtkFilename = "data/singlePhasePoromechanics_FaultModel_well_seq/extractAndMergeWell1.vtu"
         elif datasetType == "extractAndMergeVolumeWell1":
             vtkFilename = "data/singlePhasePoromechanics_FaultModel_well_seq/extractAndMergeVolumeWell1.vtm"
+        elif datasetType == "quads2_tris4":
+            vtkFilename = "data/quads2_tris4.vtu"
+        elif datasetType == "hexs3_tets36_pyrs18":
+            vtkFilename = "data/hexs3_tets36_pyrs18.vtu"
 
         datapath: str = str( Path( __file__ ).parent / vtkFilename )
         reader.SetFileName( datapath )

geos-processing/tests/test_CellTypeCounterEnhanced.py → geos-mesh/tests/test_CellTypeCounterEnhanced.py

@@ -6,7 +6,7 @@
 import pytest
 from typing import Any
 
-from geos.processing.pre_processing.CellTypeCounterEnhanced import CellTypeCounterEnhanced
+from geos.mesh.stats.CellTypeCounterEnhanced import CellTypeCounterEnhanced
 from geos.mesh.model.CellTypeCounts import CellTypeCounts
 
 from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkCellTypes, vtkCell, VTK_TRIANGLE, VTK_QUAD,

geos-processing/tests/test_SplitMesh.py → geos-mesh/tests/test_SplitMesh.py

@@ -12,8 +12,8 @@
                                             VTK_QUAD, VTK_TETRA, VTK_HEXAHEDRON, VTK_PYRAMID )
 from vtkmodules.vtkCommonCore import vtkPoints, vtkIdList, vtkDataArray
 
-from geos.mesh.utils.genericHelpers import createSingleCellMesh
-from geos.processing.generic_processing_tools.SplitMesh import SplitMesh
+from geos.mesh.utils.genericHelpers import createMultiCellMesh, createSingleCellMesh
+from geos.mesh.utils.SplitMesh import SplitMesh
 
 ###############################################################
 #                  create single tetra mesh                   #
@@ -296,3 +296,86 @@ def test_splitMesh(
     nbArrayInput: int = cellDataInput.GetNumberOfArrays()
     nbArraySplitted: int = cellData.GetNumberOfArrays()
     assert nbArraySplitted == nbArrayInput + 1, f"Number of arrays should be { nbArrayInput + 1 }."
+
+
+@pytest.mark.parametrize( "cellType,coords", [
+    ( VTK_HEXAHEDRON, [
+        np.array( [ [ 0., 0., 0. ], [ 1., 0., 0. ], [ 1., 1., 0. ], [ 0., 1., 0. ], [ 0., 0., 1. ], [ 1., 0., 1. ],
+                    [ 1., 1., 1. ], [ 0., 1., 1. ] ] ),
+        np.array( [ [ 1., 0., 0. ], [ 2., 0., 0. ], [ 2., 1., 0. ], [ 1., 1., 0. ], [ 1., 0., 1. ], [ 2., 0., 1. ],
+                    [ 2., 1., 1. ], [ 1., 1., 1. ] ] ),
+    ] ),
+    ( VTK_TETRA, [
+        np.array( [ [ 0., 0., 0. ], [ 1., 0., 0. ], [ 0., 1., 0. ], [ 0., 0., 1. ] ] ),
+        np.array( [ [ 1., 0., 0. ], [ 1., 1., 0. ], [ 0., 1., 0. ], [ 0., 0., 1. ] ] ),
+    ] ),
+    ( VTK_TRIANGLE, [
+        np.array( [ [ 0., 0., 0. ], [ 1., 0., 0. ], [ 0., 1., 0. ] ] ),
+        np.array( [ [ 1., 0., 0. ], [ 1., 1., 0. ], [ 0., 1., 0. ] ] ),
+    ] ),
+    ( VTK_QUAD, [
+        np.array( [ [ 0., 0., 0. ], [ 1., 0., 0. ], [ 1., 1., 0. ], [ 0., 1., 0. ] ] ),
+        np.array( [ [ 1., 0., 0. ], [ 2., 0., 0. ], [ 2., 1., 0. ], [ 1., 1., 0. ] ] ),
+    ] ),
+] )
+def test_splitMeshFaceConformity( cellType: int, coords: list[ npt.NDArray[ np.float64 ] ] ) -> None:
+    """Two adjacent cells sharing a face must produce topologically conforming children.
+
+    After splitting, no two distinct point IDs should occupy the same coordinates.
+    If the edge-midpoint cache works correctly, shared-edge midpoints are inserted
+    once and reused, so unique coordinates == number of points in the output.
+    """
+    cellTypeName: str = vtkCellTypes.GetClassNameFromTypeId( cellType )
+    mesh: vtkUnstructuredGrid = createMultiCellMesh( [ cellType, cellType ], coords, sharePoints=True )
+
+    splitFilter: SplitMesh = SplitMesh( mesh )
+    splitFilter.applyFilter()
+    output: vtkUnstructuredGrid = splitFilter.getOutput()
+
+    pts: npt.NDArray[ np.float64 ] = vtk_to_numpy( output.GetPoints().GetData() )
+    uniquePts = np.unique( pts, axis=0 )
+
+    assert len( uniquePts ) == output.GetNumberOfPoints(), (
+        f"{ cellTypeName }: found { output.GetNumberOfPoints() - len( uniquePts ) } duplicate coincident points after "
+        f"splitting -- shared-edge midpoints are not being reused." )
+
+
+def test_splitMeshHybrid2D3DConformity() -> None:
+    """A 3D cell with an explicit 2D face cell must split conformally.
+
+    Build a single tetrahedron (pts 0-3) whose bottom face (pts 0,1,2) is also
+    represented as an explicit triangle cell (as fracture or boundary faces are
+    in GEOS meshes). After splitting, the 2D children must share the same point
+    IDs as the corresponding face children of the 3D split -- no duplicate
+    coincident points, and the output has the correct cell counts.
+    """
+    pts = np.array( [ [ 0., 0., 0. ], [ 1., 0., 0. ], [ 0., 1., 0. ], [ 0., 0., 1. ] ], dtype=np.float64 )
+
+    mesh: vtkUnstructuredGrid = vtkUnstructuredGrid()
+    points = vtkPoints()
+    for p in pts:
+        points.InsertNextPoint( p[ 0 ], p[ 1 ], p[ 2 ] )
+    mesh.SetPoints( points )
+
+    tetIds = vtkIdList()
+    for i in [ 0, 1, 2, 3 ]:
+        tetIds.InsertNextId( i )
+    mesh.InsertNextCell( VTK_TETRA, tetIds )
+
+    triIds = vtkIdList()
+    for i in [ 0, 1, 2 ]:  # bottom face of the tet
+        triIds.InsertNextId( i )
+    mesh.InsertNextCell( VTK_TRIANGLE, triIds )
+
+    splitFilter: SplitMesh = SplitMesh( mesh )
+    splitFilter.applyFilter()
+    output: vtkUnstructuredGrid = splitFilter.getOutput()
+
+    # 1 tet -> 8 tets, 1 tri -> 4 tris
+    assert output.GetNumberOfCells() == 12, f"Expected 12 cells, got { output.GetNumberOfCells() }."
+
+    outPts: npt.NDArray[ np.float64 ] = vtk_to_numpy( output.GetPoints().GetData() )
+    uniquePts = np.unique( outPts, axis=0 )
+    assert len( uniquePts ) == output.GetNumberOfPoints(), (
+        f"Found { output.GetNumberOfPoints() - len( uniquePts ) } duplicate coincident points -- "
+        f"3D/2D shared-face edges are not sharing midpoint IDs." )

geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py

@@ -16,7 +16,7 @@
                                             VTK_POLYHEDRON )
 from vtkmodules.util.numpy_support import vtk_to_numpy, numpy_to_vtk
 
-from geos.processing.pre_processing.CellTypeCounterEnhanced import CellTypeCounterEnhanced
+from geos.mesh.stats.CellTypeCounterEnhanced import CellTypeCounterEnhanced
 from geos.mesh.model.CellTypeCounts import CellTypeCounts
 from geos.mesh.model.QualityMetricSummary import ( QualityMetricSummary, StatTypes )
 from geos.mesh.utils.arrayHelpers import getAttributesWithNumberOfComponents

geos-processing/tests/conftest.py

@@ -58,9 +58,9 @@ def _get_dataset( datasetType: str ) -> vtkMultiBlockDataSet | vtkPolyData | vtk
             vtkFilename = "geos-mesh/tests/data/hasFault.vtu"
         # Small useful meshes
         elif datasetType == "quads2_tris4":
-            vtkFilename = "geos-processing/tests/data/quads2_tris4.vtu"
+            vtkFilename = "geos-mesh/tests/data/quads2_tris4.vtu"
         elif datasetType == "hexs3_tets36_pyrs18":
-            vtkFilename = "geos-processing/tests/data/hexs3_tets36_pyrs18.vtu"
+            vtkFilename = "geos-mesh/tests/data/hexs3_tets36_pyrs18.vtu"
         elif datasetType == "meshtet1":
             vtkFilename = "geos-processing/tests/data/mesh1.vtu"
         elif datasetType == "meshtet1b":

geos-pv/src/geos/pv/plugins/generic_processing/PVSplitMesh.py

@@ -21,7 +21,7 @@
 
 update_paths()
 
-from geos.processing.generic_processing_tools.SplitMesh import SplitMesh
+from geos.mesh.utils.SplitMesh import SplitMesh
 from geos.utils.Logger import isHandlerInLogger
 from geos.pv.utils.details import ( SISOFilter, FilterCategory )