From 25455ef58c3d88b9d9f3095e1a40f229f8d636db Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Wed, 9 Sep 2026 12:56:33 -0500 Subject: [PATCH] Draft --- docs/mesh-doctor.rst | 17 + .../geos/mesh_doctor/actions/cureOneSided.py | 515 ++++++++++++++++++ .../src/geos/mesh_doctor/parsing/__init__.py | 1 + .../parsing/cureOneSidedParsing.py | 171 ++++++ mesh-doctor/src/geos/mesh_doctor/register.py | 3 +- mesh-doctor/tests/test_cureOneSided.py | 194 +++++++ 6 files changed, 900 insertions(+), 1 deletion(-) create mode 100644 mesh-doctor/src/geos/mesh_doctor/actions/cureOneSided.py create mode 100644 mesh-doctor/src/geos/mesh_doctor/parsing/cureOneSidedParsing.py create mode 100644 mesh-doctor/tests/test_cureOneSided.py diff --git a/docs/mesh-doctor.rst b/docs/mesh-doctor.rst index 06a33045..7dc5a030 100644 --- a/docs/mesh-doctor.rst +++ b/docs/mesh-doctor.rst @@ -125,6 +125,23 @@ The ``generateFractures`` module will split the mesh and generate the multi-bloc .. command-output:: mesh-doctor generateFractures --help :shell: +``cureOneSided`` +"""""""""""""""" + +Node splitting (``generateFractures``) duplicates the fault nodes so that both sides can move independently, +but it writes a single surface polygon per fault location, remapped onto whichever 3D neighbour it enumerated first. +The fault surface is then a patchwork: adjacent faces sit on different collocated node copies. +Since ``geos`` derives the seal node set from that surface, the open matrix/fracture taps end up on mixed sides +and fluid crosses the fault even though most faces look sealed. + +The ``cureOneSided`` module rewrites every tagged fault face onto the coincident real 3D-cell face of one +consistent side per fault value. Two cell arrays are added for display and QC: + ``faultSide`` (``+1`` reference side, ``-1`` residual degeneracy, ``0`` non-fault) +and ``onHole`` (``1`` on faces bordering a residual hole). A fully cured fault has ``onHole`` all ``0``. + +.. command-output:: mesh-doctor cureOneSided --help + :shell: + ``generateGlobalIds`` """"""""""""""""""""" diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/cureOneSided.py b/mesh-doctor/src/geos/mesh_doctor/actions/cureOneSided.py new file mode 100644 index 00000000..454f6687 --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/actions/cureOneSided.py @@ -0,0 +1,515 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +"""Cure the one-sided fault issue on a split (post-fracture-generation) mesh. + +Node splitting duplicates the fault nodes so that both sides can move independently, but it +writes a single surface polygon per fault location, remapped onto whichever 3D neighbour was +enumerated first. The resulting fault surface is a patchwork: adjacent faces sit on different +collocated node copies. GEOS derives its seal node set from that surface, so the open +matrix/fracture taps end up on mixed sides and fluid crosses the fault even though most faces +look sealed. + +This action rewrites every tagged fault face onto ONE consistent side per fault value, so that: + - the seal node set holds only that side's node copies (cross-fault flow is blocked, and the + fracture stays anchored to one matrix side so the flow solve converges); + - every fault face remains a real 3D-cell face (the orphan2d action still passes); + - fault junctions (nodes with 3 or 4 collocated copies) are handled by locating the coincident + reference-side 3D face geometrically, not through a two-way twin relationship. + +Collocated (split) nodes are detected by coordinate coincidence, so no separate fault or +faceBlock file is needed. Two cell arrays are added for display and QC: + - "faultSide": +1 on the reference side, -1 on a residual degeneracy (typically a fault-fault + intersection line, where no single side exists), 0 on non-fault cells. + - "onHole": 1 on faces bordering a residual hole (a boundary loop other than the fault + perimeter), 0 elsewhere. A fully cured fault has onHole all 0. +""" + +from __future__ import annotations + +from collections import defaultdict +from dataclasses import dataclass +from typing import Iterable + +import numpy as np +from numpy.typing import NDArray + +import vtk +from vtkmodules.vtkCommonCore import vtkIdList +from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkStaticPointLocator, VTK_TRIANGLE, VTK_POLYGON, + VTK_QUAD ) +from vtkmodules.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray, vtk_to_numpy + +from geos.mesh.io.vtkIO import VtkOutput, readUnstructuredGrid, writeMesh +from geos.mesh_doctor.parsing.cliParsing import setupLogger + +FAULT_SIDE_ARRAY: str = "faultSide" +ON_HOLE_ARRAY: str = "onHole" + +__CELL_TYPES_2D: tuple[ int, ...] = ( VTK_TRIANGLE, VTK_QUAD, VTK_POLYGON ) +__IS_FAULT_FACE_ARRAY: str = "__cureOneSidedIsFaultFace" +__ORIGINAL_CELL_ARRAY: str = "__cureOneSidedOriginalCell" + + +@dataclass( frozen=True ) +class Options: + """Options for the cureOneSided action. + + Attributes: + outputFile: VTK output file configuration for the cured mesh. + tagArray: Name of the cell-data array tagging the fault faces. + tagValues: Fault values to cure. Empty means every distinct non-zero value of tagArray + carried by a 2D cell. + tolerance: Distance below which two nodes are considered collocated. + """ + outputFile: VtkOutput + tagArray: str + tagValues: tuple[ int, ...] + tolerance: float + + +@dataclass( frozen=True ) +class FaultResult: + """Per-fault-value outcome of the cure. + + Attributes: + tagValue: The fault value this entry describes. + numFaces: Number of 2D cells carrying this value. + numMovedFaces: Faces rewritten onto the reference side. + numAlreadyCorrectFaces: Faces that already sat on the reference side. + numDegenerateFaces: Faces left as-is because no reference-side 3D face exists. + numSkippedFaces: Faces with no orientable normal or no 3D owner cell (orphan faces). + numHoleBorderFaces: Faces bordering a residual hole of this fault surface. + """ + tagValue: int + numFaces: int + numMovedFaces: int + numAlreadyCorrectFaces: int + numDegenerateFaces: int + numSkippedFaces: int + numHoleBorderFaces: int + + +@dataclass( frozen=True ) +class Result: + """Result of the cureOneSided action. + + Attributes: + faults: Per-fault-value outcomes, ordered by fault value. + numFaultFaces: Total number of tagged 2D fault faces. + numMovedFaces: Total number of faces rewritten onto the reference side. + numDegenerateFaces: Total number of faces left on their original side. + numSkippedFaces: Total number of faces that could not be processed. + numHoleBorderFaces: Total number of faces bordering a residual hole. + """ + faults: tuple[ FaultResult, ...] + numFaultFaces: int + numMovedFaces: int + numDegenerateFaces: int + numSkippedFaces: int + numHoleBorderFaces: int + + +class _UnionFind: + """Minimal union-find over point ids, used to group boundary edges into loops.""" + + def __init__( self ) -> None: + """Create an empty union-find structure.""" + self.m_parent: dict[ int, int ] = {} + + def find( self, a: int ) -> int: + """Return the representative of the set containing ``a``. + + Args: + a: The element to look up. It is inserted if unknown. + + Returns: + The representative element of ``a``'s set. + """ + self.m_parent.setdefault( a, a ) + while self.m_parent[ a ] != a: + self.m_parent[ a ] = self.m_parent[ self.m_parent[ a ] ] + a = self.m_parent[ a ] + return a + + def union( self, a: int, b: int ) -> None: + """Merge the sets containing ``a`` and ``b``. + + Args: + a: First element. + b: Second element. + """ + rootA, rootB = self.find( a ), self.find( b ) + if rootA != rootB: + self.m_parent[ rootA ] = rootB + + +def _toIdList( idList: vtkIdList, pointIds: NDArray[ np.int64 ] ) -> vtkIdList: + """Fill a vtkIdList with the given point ids and return it. + + Args: + idList: The list to reset and fill (reused across calls to avoid churn). + pointIds: The point ids to insert. + + Returns: + The filled list. + """ + idList.Reset() + for pointId in pointIds: + idList.InsertNextId( int( pointId ) ) + return idList + + +def buildCollocatedGroups( mesh: vtkUnstructuredGrid, pointIds: Iterable[ int ], + tolerance: float ) -> tuple[ NDArray[ np.int64 ], dict[ int, list[ int ] ] ]: + """Group the given points with every mesh point collocated with them. + + Two points closer than ``tolerance`` belong to the same group, i.e. they are two copies of a + single geometric node created by the fault split. Only the groups reachable from ``pointIds`` + are built, which is all the cure needs and much cheaper than grouping the whole mesh. + + Args: + mesh: The mesh whose points are grouped. + pointIds: The points to group, typically the corners of the fault faces. + tolerance: Distance below which two points are considered collocated. + + Returns: + A ``representative`` array mapping each point id to the lowest id of its group (to itself + for a point outside of any built group), and a mapping from that representative id to every + point id of the group. + """ + coordinates: NDArray[ np.float64 ] = vtk_to_numpy( mesh.GetPoints().GetData() ) + locator = vtkStaticPointLocator() + locator.SetDataSet( mesh ) + locator.BuildLocator() + + # Union-find rather than one bucket per query, so that a chain of near-coincident copies whose + # ends are further apart than the tolerance still ends up in a single group. + unionFind = _UnionFind() + neighbors = vtkIdList() + queried: set[ int ] = set() + for pointId in pointIds: + if pointId in queried: + continue + queried.add( pointId ) + locator.FindPointsWithinRadius( tolerance, coordinates[ pointId ], neighbors ) + unionFind.find( pointId ) + for k in range( neighbors.GetNumberOfIds() ): + unionFind.union( pointId, neighbors.GetId( k ) ) + + representative: NDArray[ np.int64 ] = np.arange( mesh.GetNumberOfPoints(), dtype=np.int64 ) + members: dict[ int, list[ int ] ] = defaultdict( list ) + for pointId in unionFind.m_parent: + members[ unionFind.find( pointId ) ].append( pointId ) + groups: dict[ int, list[ int ] ] = {} + for group in members.values(): + group.sort() + root = group[ 0 ] + representative[ group ] = root + groups[ root ] = group + return representative, groups + + +def computeFaceNormals( mesh: vtkUnstructuredGrid, tagArray: str, tagValue: int ) -> dict[ int, NDArray[ np.float64 ] ]: + """Compute a consistently oriented normal for every 2D face carrying a single fault value. + + Orientation is propagated per fault value, so that a junction between two faults does not force + the two surfaces to share an orientation. + + Args: + mesh: The mesh holding the fault faces. + tagArray: Name of the cell-data array tagging the fault faces. + tagValue: The fault value to orient. + + Returns: + A mapping from cell id to unit normal, restricted to the faces of that fault value. Faces + that the geometry extraction dropped are absent. + """ + tags = vtk_to_numpy( mesh.GetCellData().GetArray( tagArray ) ).astype( int ) + cellTypes = vtk_to_numpy( mesh.GetCellTypesArray() ) + isFaultFace = ( ( tags == tagValue ) & np.isin( cellTypes, __CELL_TYPES_2D ) ).astype( np.int8 ) + if isFaultFace.sum() == 0: + return {} + + faultFaceArray = numpy_to_vtk( isFaultFace, deep=True ) + faultFaceArray.SetName( __IS_FAULT_FACE_ARRAY ) + mesh.GetCellData().AddArray( faultFaceArray ) + originalCellArray = numpy_to_vtk( np.arange( mesh.GetNumberOfCells(), dtype=np.int64 ), deep=True ) + originalCellArray.SetName( __ORIGINAL_CELL_ARRAY ) + mesh.GetCellData().AddArray( originalCellArray ) + try: + threshold = vtk.vtkThreshold() + threshold.SetInputData( mesh ) + threshold.SetInputArrayToProcess( 0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS, __IS_FAULT_FACE_ARRAY ) + threshold.SetLowerThreshold( 0.5 ) + threshold.SetUpperThreshold( 1.5 ) + threshold.SetThresholdFunction( vtk.vtkThreshold.THRESHOLD_BETWEEN ) + threshold.Update() + + geometry = vtk.vtkGeometryFilter() + geometry.SetInputData( threshold.GetOutput() ) + geometry.Update() + + normalsFilter = vtk.vtkPolyDataNormals() + normalsFilter.SetInputData( geometry.GetOutput() ) + normalsFilter.SetConsistency( True ) + normalsFilter.SetAutoOrientNormals( False ) + normalsFilter.SetComputeCellNormals( True ) + normalsFilter.SetComputePointNormals( False ) + normalsFilter.SetSplitting( False ) + normalsFilter.Update() + oriented = normalsFilter.GetOutput() + finally: + mesh.GetCellData().RemoveArray( __IS_FAULT_FACE_ARRAY ) + mesh.GetCellData().RemoveArray( __ORIGINAL_CELL_ARRAY ) + + normals = vtk_to_numpy( oriented.GetCellData().GetNormals() ) + originalCells = vtk_to_numpy( oriented.GetCellData().GetArray( __ORIGINAL_CELL_ARRAY ) ).astype( np.int64 ) + return { int( originalCells[ i ] ): normals[ i ] for i in range( len( originalCells ) ) } + + +def markHoleBorderFaces( numCells: int, facesByValue: dict[ int, list[ int ] ], offsets: NDArray[ np.int64 ], + connectivity: NDArray[ np.int64 ] ) -> NDArray[ np.int8 ]: + """Flag the faces that border a residual hole of their fault surface. + + A hole is any boundary-edge loop of a fault other than its outer perimeter. Boundary edges are + computed from the cured raw connectivity, so coincident collocated points stay distinct and the + seams are not welded shut. + + Args: + numCells: Total number of cells in the mesh. + facesByValue: Fault value to the list of its 2D cell ids. + offsets: Cell connectivity offsets of the (cured) mesh. + connectivity: Cell connectivity of the (cured) mesh. + + Returns: + A per-cell flag, 1 on the faces bordering a residual hole and 0 elsewhere. + """ + onHole: NDArray[ np.int8 ] = np.zeros( numCells, dtype=np.int8 ) + for faces in facesByValue.values(): + edgeUseCount: dict[ tuple[ int, int ], int ] = defaultdict( int ) + edgeOwner: dict[ tuple[ int, int ], int ] = {} + for cellId in faces: + start, end = int( offsets[ cellId ] ), int( offsets[ cellId + 1 ] ) + facePoints = [ int( pointId ) for pointId in connectivity[ start:end ] ] + numFacePoints = len( facePoints ) + for k in range( numFacePoints ): + a, b = facePoints[ k ], facePoints[ ( k + 1 ) % numFacePoints ] + edge = ( a, b ) if a < b else ( b, a ) + edgeUseCount[ edge ] += 1 + edgeOwner[ edge ] = cellId + boundaryEdges = [ edge for edge, count in edgeUseCount.items() if count == 1 ] + if not boundaryEdges: + continue + + unionFind = _UnionFind() + for a, b in boundaryEdges: + unionFind.union( a, b ) + loops: dict[ int, list[ tuple[ int, int ] ] ] = defaultdict( list ) + for edge in boundaryEdges: + loops[ unionFind.find( edge[ 0 ] ) ].append( edge ) + # The longest boundary loop is the fault perimeter; any other loop is a hole. + perimeter = max( loops, key=lambda root: len( loops[ root ] ) ) + for root, edges in loops.items(): + if root == perimeter: + continue + for edge in edges: + onHole[ edgeOwner[ edge ] ] = 1 + return onHole + + +def _findTagValues( mesh: vtkUnstructuredGrid, tagArray: str, tagValues: tuple[ int, ...] ) -> list[ int ]: + """Return the fault values to cure, defaulting to every non-zero value carried by a 2D cell. + + Args: + mesh: The mesh to inspect. + tagArray: Name of the cell-data array tagging the fault faces. + tagValues: The requested values, possibly empty. + + Returns: + The sorted list of fault values to process. + """ + if tagValues: + return sorted( set( tagValues ) ) + tags = vtk_to_numpy( mesh.GetCellData().GetArray( tagArray ) ).astype( int ) + cellTypes = vtk_to_numpy( mesh.GetCellTypesArray() ) + found = np.unique( tags[ np.isin( cellTypes, __CELL_TYPES_2D ) ] ) + return [ int( value ) for value in found if value != 0 ] + + +def meshAction( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Rewrite every wrong-side fault face onto the coincident reference-side 3D cell face. + + The mesh is modified in place: the fault face connectivity is rewritten and the "faultSide" and + "onHole" cell arrays are added. + + Args: + mesh: The post-split domain mesh to cure. + options: The cure options. + + Returns: + The per-fault and overall statistics of the cure. + + Raises: + ValueError: If ``options.tagArray`` is not a cell-data array of the mesh. + """ + if mesh.GetCellData().GetArray( options.tagArray ) is None: + raise ValueError( f"Cell array \"{options.tagArray}\" is not in the mesh." ) + values: list[ int ] = _findTagValues( mesh, options.tagArray, options.tagValues ) + if not values: + setupLogger.warning( f"No non-zero value of \"{options.tagArray}\" is carried by a 2D cell." ) + + mesh.BuildLinks() + numCells: int = mesh.GetNumberOfCells() + points: NDArray[ np.float64 ] = vtk_to_numpy( mesh.GetPoints().GetData() ) + tags = vtk_to_numpy( mesh.GetCellData().GetArray( options.tagArray ) ).astype( int ) + cellTypes = vtk_to_numpy( mesh.GetCellTypesArray() ) + is2d = np.isin( cellTypes, __CELL_TYPES_2D ) + offsets: NDArray[ np.int64 ] = vtk_to_numpy( mesh.GetCells().GetOffsetsArray() ).astype( np.int64 ) + connectivity: NDArray[ np.int64 ] = vtk_to_numpy( mesh.GetCells().GetConnectivityArray() ).astype( np.int64 ).copy() + + facesByValue: dict[ int, list[ int ] ] = { + tagValue: [ int( cellId ) for cellId in np.where( ( tags == tagValue ) & is2d )[ 0 ] ] + for tagValue in values + } + faultPointIds: set[ int ] = set() + for faces in facesByValue.values(): + for cellId in faces: + faultPointIds.update( + int( pointId ) for pointId in connectivity[ int( offsets[ cellId ] ):int( offsets[ cellId + 1 ] ) ] ) + setupLogger.info( f"Grouping the {len( faultPointIds )} fault nodes collocated within {options.tolerance}." ) + representative, collocatedGroups = buildCollocatedGroups( mesh, sorted( faultPointIds ), options.tolerance ) + + faceIds, neighborIds, cellPointIds, pointCellIds = vtkIdList(), vtkIdList(), vtkIdList(), vtkIdList() + + def centroid( cellId: int ) -> NDArray[ np.float64 ]: + mesh.GetCellPoints( cellId, cellPointIds ) + return points[ [ cellPointIds.GetId( k ) for k in range( cellPointIds.GetNumberOfIds() ) ] ].mean( axis=0 ) + + faultSide: NDArray[ np.int8 ] = np.zeros( numCells, dtype=np.int8 ) + perValueCounts: dict[ int, list[ int ] ] = {} + for tagValue in values: + normals = computeFaceNormals( mesh, options.tagArray, tagValue ) + faces = facesByValue[ tagValue ] + numMoved, numAlreadyCorrect, numDegenerate, numSkipped = 0, 0, 0, 0 + setupLogger.info( f"Curing {len( faces )} faces of \"{options.tagArray}\" == {tagValue}." ) + for cellId in faces: + normal = normals.get( cellId ) + if normal is None: + numSkipped += 1 + continue + start, end = int( offsets[ cellId ] ), int( offsets[ cellId + 1 ] ) + facePointIds = connectivity[ start:end ] + faceCentroid = points[ facePointIds ].mean( axis=0 ) + + mesh.GetCellNeighbors( cellId, _toIdList( faceIds, facePointIds ), neighborIds ) + owner = -1 + for j in range( neighborIds.GetNumberOfIds() ): + neighborId = neighborIds.GetId( j ) + if mesh.GetCell( neighborId ).GetCellDimension() == 3: + owner = neighborId + break + if owner < 0: + numSkipped += 1 + continue + if np.dot( centroid( owner ) - faceCentroid, normal ) > 0: # Already on the reference side. + faultSide[ cellId ] = 1 + numAlreadyCorrect += 1 + continue + + # Look for the real 3D-cell face coincident with this one on the +normal side. Scanning + # every collocated copy of every corner also resolves fault junctions, where a node has + # 3 or 4 copies and no two-way twin exists. + target = frozenset( int( representative[ pointId ] ) for pointId in facePointIds ) + candidates: set[ int ] = set() + for pointId in facePointIds: + for collocatedId in collocatedGroups[ int( representative[ pointId ] ) ]: + mesh.GetPointCells( collocatedId, pointCellIds ) + for j in range( pointCellIds.GetNumberOfIds() ): + candidates.add( pointCellIds.GetId( j ) ) + + referenceFace: list[ int ] | None = None + for candidateId in candidates: + if candidateId == owner: + continue + candidate = mesh.GetCell( candidateId ) + if candidate.GetCellDimension() != 3: + continue + if np.dot( centroid( candidateId ) - faceCentroid, normal ) <= 0: + continue + for faceIndex in range( candidate.GetNumberOfFaces() ): + candidatePointIds = candidate.GetFace( faceIndex ).GetPointIds() + numCandidatePoints = candidatePointIds.GetNumberOfIds() + if numCandidatePoints != ( end - start ): + continue + candidateIds = [ candidatePointIds.GetId( k ) for k in range( numCandidatePoints ) ] + if frozenset( int( representative[ pointId ] ) for pointId in candidateIds ) == target: + referenceFace = candidateIds + break + if referenceFace is not None: + break + + if referenceFace is None: + # No single side exists here, typically along a fault-fault intersection line. + faultSide[ cellId ] = -1 + numDegenerate += 1 + else: + connectivity[ start:end ] = np.array( referenceFace, dtype=np.int64 ) + faultSide[ cellId ] = 1 + numMoved += 1 + perValueCounts[ tagValue ] = [ numMoved, numAlreadyCorrect, numDegenerate, numSkipped ] + + onHole = markHoleBorderFaces( numCells, facesByValue, offsets, connectivity ) + + cells = vtk.vtkCellArray() + cells.SetData( numpy_to_vtkIdTypeArray( np.ascontiguousarray( offsets ), deep=True ), + numpy_to_vtkIdTypeArray( np.ascontiguousarray( connectivity ), deep=True ) ) + mesh.SetCells( numpy_to_vtk( np.ascontiguousarray( cellTypes ), deep=True, array_type=vtk.VTK_UNSIGNED_CHAR ), + cells ) + + faultSideArray = numpy_to_vtk( faultSide, deep=True ) + faultSideArray.SetName( FAULT_SIDE_ARRAY ) + mesh.GetCellData().AddArray( faultSideArray ) + onHoleArray = numpy_to_vtk( onHole, deep=True ) + onHoleArray.SetName( ON_HOLE_ARRAY ) + mesh.GetCellData().AddArray( onHoleArray ) + for data, arrayName in ( ( mesh.GetPointData(), "GLOBAL_IDS_POINTS" ), ( mesh.GetCellData(), "GLOBAL_IDS_CELLS" ) ): + if data.GetArray( arrayName ) is not None: + data.SetGlobalIds( data.GetArray( arrayName ) ) + + faults: list[ FaultResult ] = [] + for tagValue in values: + numMoved, numAlreadyCorrect, numDegenerate, numSkipped = perValueCounts[ tagValue ] + faces = facesByValue[ tagValue ] + faults.append( + FaultResult( tagValue=tagValue, + numFaces=len( faces ), + numMovedFaces=numMoved, + numAlreadyCorrectFaces=numAlreadyCorrect, + numDegenerateFaces=numDegenerate, + numSkippedFaces=numSkipped, + numHoleBorderFaces=int( onHole[ faces ].sum() ) if faces else 0 ) ) + + return Result( faults=tuple( faults ), + numFaultFaces=sum( fault.numFaces for fault in faults ), + numMovedFaces=sum( fault.numMovedFaces for fault in faults ), + numDegenerateFaces=sum( fault.numDegenerateFaces for fault in faults ), + numSkippedFaces=sum( fault.numSkippedFaces for fault in faults ), + numHoleBorderFaces=sum( fault.numHoleBorderFaces for fault in faults ) ) + + +def action( vtuInputFile: str, options: Options ) -> Result: + """Read a split VTU mesh, cure its one-sided fault issue and write the result. + + Args: + vtuInputFile: Path to the post-split domain VTU mesh. + options: The cure options. + + Returns: + The per-fault and overall statistics of the cure. + """ + setupLogger.info( f"Reading mesh from \"{vtuInputFile}\"." ) + mesh = readUnstructuredGrid( vtuInputFile ) + result = meshAction( mesh, options ) + setupLogger.info( f"Writing cured mesh to \"{options.outputFile.output}\"." ) + writeMesh( mesh, options.outputFile ) + return result diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py b/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py index 24738719..e448206d 100644 --- a/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py @@ -21,6 +21,7 @@ EULER = "euler" CONVERT_MD2SG = "convertMD2SG" REFINE_MESH = "refineMesh" +CURE_ONE_SIDED = "cureOneSided" @dataclass( frozen=True ) diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/cureOneSidedParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/cureOneSidedParsing.py new file mode 100644 index 00000000..3cd48b37 --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/cureOneSidedParsing.py @@ -0,0 +1,171 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +"""Command line parsing for the cureOneSided action.""" + +from __future__ import annotations + +import argparse +from typing import Any + +from geos.mesh_doctor.actions.cureOneSided import Options, Result +from geos.mesh_doctor.parsing import CURE_ONE_SIDED, vtkOutputParsing +from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument + +__TAG_ARRAY = "tagArray" +__TAG_VALUES = "tagValues" +__TOLERANCE = "tolerance" +__TOLERANCE_DEFAULT = 1.e-6 + + +def parseTagValues( spec: str ) -> tuple[ int, ...]: + """Parse a fault value specification into a sorted tuple of integers. + + Args: + spec: A comma-separated list of values and inclusive ranges, e.g. "1-7", "18,19,20" or + "1,4-6". An empty string yields an empty tuple. + + Returns: + The sorted, deduplicated values. + + Raises: + ValueError: If a range is malformed or a value is not an integer. + """ + values: set[ int ] = set() + for part in spec.split( "," ): + part = part.strip() + if not part: + continue + separator = part.find( "-", 1 ) # A leading minus is a negative value, not a range separator. + if separator < 0: + values.add( int( part ) ) + else: + low, high = int( part[ :separator ] ), int( part[ separator + 1: ] ) + if high < low: + raise ValueError( f"Invalid range \"{part}\": the upper bound is below the lower bound." ) + values.update( range( low, high + 1 ) ) + return tuple( sorted( values ) ) + + +def fillSubparser( subparsers: argparse._SubParsersAction[ Any ] ) -> None: + """Fill the argument parser for the cureOneSided action. + + Args: + subparsers: The subparsers action to add the parser to. + """ + p = subparsers.add_parser( + CURE_ONE_SIDED, + help="Rewrite the fault faces of a split mesh onto one consistent side per fault.", + description="""\ +Cure the one-sided fault issue on a split (post-generateFractures) mesh. + +Node splitting duplicates the fault nodes, but writes one surface polygon per fault +location remapped onto whichever 3D neighbour was enumerated first. The fault surface +is then a patchwork: adjacent faces sit on different collocated node copies, the seal +node set GEOS derives from it holds mixed sides, and fluid crosses the fault even +though most faces look sealed. + +This action rewrites every tagged fault face onto the coincident real 3D-cell face of +ONE consistent side per fault value, so that: + - all the open matrix/fracture taps land on the same side (cross-fault flow blocked, + and the fracture stays anchored to one matrix side so the flow solve converges); + - every fault face stays a real 3D-cell face (orphan2d still passes); + - fault junctions (nodes with 3-4 collocated copies) are resolved geometrically. + +Collocated nodes are detected by coordinate coincidence, so no separate fault or +faceBlock file is needed. Global id arrays are preserved. Two cell arrays are added +for display and QC: + faultSide +1 reference side, -1 residual degeneracy (e.g. a fault-fault + intersection line, where no single side exists), 0 on non-fault cells. + onHole 1 on faces bordering a residual hole (a boundary loop other than the + fault perimeter), 0 elsewhere. A fully cured fault has onHole all 0. + +Examples: + mesh-doctor cureOneSided -i domain.vtu --output cured.vtu --tagArray FaultMask --tagValues 1-7 + mesh-doctor cureOneSided -i domain.vtu --output cured.vtu --tagArray attribute --tagValues 18-24 + mesh-doctor cureOneSided -i domain.vtu --output cured.vtu --tagArray region + +Verify the output with: + mesh-doctor orphan2d -i cured.vtu + mesh-doctor euler -i cured.vtu --mode surface --tagArray FaultMask +""", + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + addVtuInputFileArgument( p ) + vtkOutputParsing.fillVtkOutputSubparser( p ) + p.add_argument( "--" + __TAG_ARRAY, + type=str, + required=True, + metavar="NAME", + help="[string]: Cell-data array tagging the fault faces, e.g. FaultMask." ) + p.add_argument( "--" + __TAG_VALUES, + type=str, + default="", + metavar="SPEC", + help="[string]: Fault values to cure, as values and inclusive ranges, e.g. " + "\"1-7\", \"18,19,20\" or \"1,4-6\". Defaults to every distinct non-zero value " + "of --tagArray carried by a 2D cell." ) + p.add_argument( "--" + __TOLERANCE, + type=float, + default=__TOLERANCE_DEFAULT, + metavar="D", + help=f"[float]: Distance below which two nodes are considered collocated, i.e. two " + f"copies of a single split node. Default: {__TOLERANCE_DEFAULT}." ) + + +def convert( parsedOptions: dict[ str, Any ] ) -> Options: + """Convert parsed command-line options to an Options object. + + Args: + parsedOptions: Dictionary of parsed command-line options. + + Returns: + Options for the cureOneSided action. + + Raises: + ValueError: If the tolerance is negative. + """ + tolerance: float = parsedOptions.get( __TOLERANCE, __TOLERANCE_DEFAULT ) + if tolerance < 0.: + raise ValueError( f"--{__TOLERANCE} must be >= 0, got {tolerance}." ) + return Options( + outputFile=vtkOutputParsing.convert( parsedOptions ), + tagArray=parsedOptions[ __TAG_ARRAY ], + tagValues=parseTagValues( parsedOptions.get( __TAG_VALUES, "" ) ), + tolerance=tolerance, + ) + + +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the cure. + + Args: + options: The options used for the cure. + result: The result of the cureOneSided action. + """ + setupLogger.results( "=" * 80 ) + setupLogger.results( "CURE ONE-SIDED FAULTS" ) + setupLogger.results( "=" * 80 ) + setupLogger.results( f"Tag array : {options.tagArray}" ) + setupLogger.results( f"Collocation tol. : {options.tolerance}" ) + setupLogger.results( f"{'value':>8} {'faces':>10} {'moved':>10} {'kept':>10} {'degenerate':>12} " + f"{'skipped':>10} {'onHole':>10}" ) + for fault in result.faults: + setupLogger.results( f"{fault.tagValue:>8} {fault.numFaces:>10,} {fault.numMovedFaces:>10,} " + f"{fault.numAlreadyCorrectFaces:>10,} {fault.numDegenerateFaces:>12,} " + f"{fault.numSkippedFaces:>10,} {fault.numHoleBorderFaces:>10,}" ) + setupLogger.results( "-" * 80 ) + setupLogger.results( f"Fault faces : {result.numFaultFaces:,}" ) + setupLogger.results( f"Moved to reference : {result.numMovedFaces:,}" ) + setupLogger.results( f"Degenerate (kept) : {result.numDegenerateFaces:,}" ) + setupLogger.results( f"Skipped (no owner) : {result.numSkippedFaces:,}" ) + setupLogger.results( f"Hole-border faces : {result.numHoleBorderFaces:,}" ) + if result.numDegenerateFaces == 0 and result.numHoleBorderFaces == 0: + setupLogger.results( "STATUS: CURED (every fault is single-sided, no residual hole)" ) + else: + setupLogger.results( "STATUS: RESIDUAL DEGENERACY - inspect the faultSide and onHole arrays. " + "A fault-fault intersection line leaves a harmless slit, not a leak." ) + if result.numSkippedFaces > 0: + setupLogger.results( f"WARNING: {result.numSkippedFaces:,} face(s) have no 3D owner cell. " + "Run orphan2d on the input mesh." ) + setupLogger.results( f"Output written to : {options.outputFile.output}" ) + setupLogger.results( "=" * 80 ) diff --git a/mesh-doctor/src/geos/mesh_doctor/register.py b/mesh-doctor/src/geos/mesh_doctor/register.py index b8a83f6e..3785f562 100644 --- a/mesh-doctor/src/geos/mesh_doctor/register.py +++ b/mesh-doctor/src/geos/mesh_doctor/register.py @@ -59,7 +59,8 @@ def registerParsingActions( parsing.FIX_ELEMENTS_ORDERINGS, parsing.GENERATE_CUBE, parsing.GENERATE_FRACTURES, parsing.GENERATE_GLOBAL_IDS, parsing.MAIN_CHECKS, parsing.NON_CONFORMAL, parsing.SELF_INTERSECTING_ELEMENTS, parsing.SUPPORTED_ELEMENTS, parsing.ORPHAN_2D, - parsing.CHECK_INTERNAL_TAGS, parsing.EULER, parsing.CONVERT_MD2SG, parsing.REFINE_MESH ): + parsing.CHECK_INTERNAL_TAGS, parsing.EULER, parsing.CONVERT_MD2SG, parsing.REFINE_MESH, + parsing.CURE_ONE_SIDED ): __HELPERS[ actionName ] = actionName __ACTIONS[ actionName ] = actionName diff --git a/mesh-doctor/tests/test_cureOneSided.py b/mesh-doctor/tests/test_cureOneSided.py new file mode 100644 index 00000000..c6592c65 --- /dev/null +++ b/mesh-doctor/tests/test_cureOneSided.py @@ -0,0 +1,194 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +"""Tests for the cureOneSided action and its CLI parsing.""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +import pytest + +from vtkmodules.util.numpy_support import numpy_to_vtk, vtk_to_numpy +from vtkmodules.vtkCommonCore import vtkPoints +from vtkmodules.vtkCommonDataModel import ( vtkCellArray, vtkUnstructuredGrid, VTK_HEXAHEDRON, VTK_QUAD ) + +from geos.mesh.io.vtkIO import VtkOutput, readUnstructuredGrid +from geos.mesh_doctor.actions.cureOneSided import ( FAULT_SIDE_ARRAY, ON_HOLE_ARRAY, Options, action, meshAction ) +from geos.mesh_doctor.parsing import cureOneSidedParsing + +# The left block holds point ids 0..11, the right block 12..23 (see __buildSplitMesh). +NUM_LEFT_POINTS: int = 12 +TAG_ARRAY: str = "FaultMask" + + +def __leftPoint( i: int, j: int, k: int ) -> int: + """Return the id of the left-block point at grid position (i, j, k).""" + return i * 6 + j * 2 + k + + +def __rightPoint( i: int, j: int, k: int ) -> int: + """Return the id of the right-block point at grid position (i, j, k).""" + return NUM_LEFT_POINTS + i * 6 + j * 2 + k + + +def __buildSplitMesh( leftQuadFirst: bool = True ) -> vtkUnstructuredGrid: + """Build a 2x2x1 hexahedral mesh split along the x = 1 plane, tagged with two fault quads. + + The two blocks (x in [0, 1] and x in [1, 2]) carry their own copies of the x = 1 nodes, as + node splitting produces. The fault surface is a deliberate patchwork: one quad sits on the + left-block copies, the other on the right-block copies. + + Args: + leftQuadFirst: If True the y in [0, 1] quad is on the left copies and the y in [1, 2] quad + on the right ones; swapped otherwise. + + Returns: + The split mesh, with a FaultMask cell array equal to 1 on the two quads and 0 elsewhere. + """ + points = vtkPoints() + for xValue in ( 0., 1. ): # Left block. + for yValue in ( 0., 1., 2. ): + for zValue in ( 0., 1. ): + points.InsertNextPoint( xValue, yValue, zValue ) + for xValue in ( 1., 2. ): # Right block, with its own copies of the x = 1 nodes. + for yValue in ( 0., 1., 2. ): + for zValue in ( 0., 1. ): + points.InsertNextPoint( xValue, yValue, zValue ) + + connectivities: list[ list[ int ] ] = [] + cellTypes: list[ int ] = [] + for corner in ( __leftPoint, __rightPoint ): + for j in ( 0, 1 ): + connectivities.append( [ + corner( 0, j, 0 ), + corner( 1, j, 0 ), + corner( 1, j + 1, 0 ), + corner( 0, j + 1, 0 ), + corner( 0, j, 1 ), + corner( 1, j, 1 ), + corner( 1, j + 1, 1 ), + corner( 0, j + 1, 1 ), + ] ) + cellTypes.append( VTK_HEXAHEDRON ) + + lowerQuadOnLeft = leftQuadFirst + for j, onLeft in ( ( 0, lowerQuadOnLeft ), ( 1, not lowerQuadOnLeft ) ): + if onLeft: + connectivities.append( [ + __leftPoint( 1, j, 0 ), + __leftPoint( 1, j + 1, 0 ), + __leftPoint( 1, j + 1, 1 ), + __leftPoint( 1, j, 1 ) + ] ) + else: + connectivities.append( [ + __rightPoint( 0, j, 0 ), + __rightPoint( 0, j + 1, 0 ), + __rightPoint( 0, j + 1, 1 ), + __rightPoint( 0, j, 1 ) + ] ) + cellTypes.append( VTK_QUAD ) + + cells = vtkCellArray() + for connectivity in connectivities: + cells.InsertNextCell( len( connectivity ), connectivity ) + + mesh = vtkUnstructuredGrid() + mesh.SetPoints( points ) + mesh.SetCells( cellTypes, cells ) + + tags = np.array( [ 0, 0, 0, 0, 1, 1 ], dtype=np.int32 ) + tagArray = numpy_to_vtk( tags, deep=True ) + tagArray.SetName( TAG_ARRAY ) + mesh.GetCellData().AddArray( tagArray ) + return mesh + + +def __faultQuadPointIds( mesh: vtkUnstructuredGrid, cellId: int ) -> list[ int ]: + """Return the point ids of the fault quad ``cellId``.""" + pointIds = mesh.GetCell( cellId ).GetPointIds() + return [ pointIds.GetId( k ) for k in range( pointIds.GetNumberOfIds() ) ] + + +def __options( outputFile: str = "unused.vtu", tagValues: tuple[ int, ...] = ( 1, ) ) -> Options: + """Build the cure options for the test mesh.""" + return Options( outputFile=VtkOutput( output=outputFile, isDataModeBinary=True ), + tagArray=TAG_ARRAY, + tagValues=tagValues, + tolerance=1.e-8 ) + + +@pytest.mark.parametrize( "leftQuadFirst", ( True, False ) ) +def test_patchworkFaultBecomesSingleSided( leftQuadFirst: bool ) -> None: + """The two fault quads end up on the same side, whichever side the cure picks as reference.""" + mesh = __buildSplitMesh( leftQuadFirst ) + faultCellIds = ( 4, 5 ) + before = [ __faultQuadPointIds( mesh, cellId ) for cellId in faultCellIds ] + assert ( min( before[ 0 ] ) < NUM_LEFT_POINTS ) != ( min( before[ 1 ] ) < NUM_LEFT_POINTS ) + + result = meshAction( mesh, __options() ) + + assert result.numFaultFaces == 2 + assert result.numMovedFaces == 1 + assert result.numDegenerateFaces == 0 + assert result.numSkippedFaces == 0 + assert result.numHoleBorderFaces == 0 + assert len( result.faults ) == 1 + assert result.faults[ 0 ].tagValue == 1 + assert result.faults[ 0 ].numAlreadyCorrectFaces == 1 + + after = [ __faultQuadPointIds( mesh, cellId ) for cellId in faultCellIds ] + sides = [ all( pointId < NUM_LEFT_POINTS for pointId in ids ) for ids in after ] + assert sides[ 0 ] == sides[ 1 ], "The two fault quads still sit on different node copies." + # The faces keep their geometry: only the node copy they refer to changes. + coordinates = vtk_to_numpy( mesh.GetPoints().GetData() ) + for ids in after: + assert np.allclose( coordinates[ ids ][ :, 0 ], 1. ) + + +def test_curedFaceIsARealCellFace() -> None: + """Every cured fault quad matches a face of one of its 3D neighbour cells.""" + mesh = __buildSplitMesh() + meshAction( mesh, __options() ) + + for cellId in ( 4, 5 ): + target = frozenset( __faultQuadPointIds( mesh, cellId ) ) + found = False + for hexId in range( 4 ): + hexCell = mesh.GetCell( hexId ) + for faceIndex in range( hexCell.GetNumberOfFaces() ): + facePointIds = hexCell.GetFace( faceIndex ).GetPointIds() + ids = frozenset( facePointIds.GetId( k ) for k in range( facePointIds.GetNumberOfIds() ) ) + found = found or ids == target + assert found, f"Fault quad {cellId} is not a face of any 3D cell." + + +def test_qcArraysAreAdded() -> None: + """The faultSide and onHole arrays flag the fault cells and leave no residual hole.""" + mesh = __buildSplitMesh() + meshAction( mesh, __options() ) + + faultSide = vtk_to_numpy( mesh.GetCellData().GetArray( FAULT_SIDE_ARRAY ) ) + onHole = vtk_to_numpy( mesh.GetCellData().GetArray( ON_HOLE_ARRAY ) ) + assert list( faultSide ) == [ 0, 0, 0, 0, 1, 1 ] + assert onHole.sum() == 0 + + +def test_tagValuesAreDetectedWhenNotGiven() -> None: + """An empty tagValues cures every distinct non-zero value carried by a 2D cell.""" + mesh = __buildSplitMesh() + result = meshAction( mesh, __options( tagValues=() ) ) + assert [ fault.tagValue for fault in result.faults ] == [ 1 ] + assert result.numFaultFaces == 2 + + +def test_unknownTagArrayRaises() -> None: + """A tag array missing from the mesh is reported as a ValueError.""" + mesh = __buildSplitMesh() + options = Options( outputFile=VtkOutput( output="unused.vtu", isDataModeBinary=True ), + tagArray="notThere", + tagValues=( 1, ), + tolerance=1.e-8 ) + with pytest.raises( ValueError ): + meshAction( mesh, options )