feat: Adding a fault stability analysis workflow

geos-geomechanics/src/geos/geomechanics/model/StressTensor.py

@@ -0,0 +1,62 @@
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2026 TotalEnergies.
+# SPDX-FileContributor: Nicolas Pillardou, Paloma Martinez
+
+import numpy as np
+import numpy.typing as npt
+from typing_extensions import Any
+
+
+# ============================================================================
+# STRESS TENSOR OPERATIONS
+# ============================================================================
+class StressTensor:
+    """Utility class for stress tensor operations."""
+
+    @staticmethod
+    def buildFromArray( arr: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
+        """Convert stress array to 3x3 tensor format."""
+        n = arr.shape[ 0 ]
+        tensors: npt.NDArray[ np.float64 ] = np.zeros( ( n, 3, 3 ), dtype=np.float64 )
+
+        if arr.shape[ 1 ] == 6:  # Voigt notation
+            tensors[ :, 0, 0 ] = arr[ :, 0 ]  # Sxx
+            tensors[ :, 1, 1 ] = arr[ :, 1 ]  # Syy
+            tensors[ :, 2, 2 ] = arr[ :, 2 ]  # Szz
+            tensors[ :, 1, 2 ] = tensors[ :, 2, 1 ] = arr[ :, 3 ]  # Syz
+            tensors[ :, 0, 2 ] = tensors[ :, 2, 0 ] = arr[ :, 4 ]  # Sxz
+            tensors[ :, 0, 1 ] = tensors[ :, 1, 0 ] = arr[ :, 5 ]  # Sxy
+        elif arr.shape[ 1 ] == 9:
+            tensors = arr.reshape( ( -1, 3, 3 ) )
+        else:
+            raise ValueError( f"Unsupported stress shape: {arr.shape}" )
+
+        return tensors
+
+    @staticmethod
+    def rotateToFaultFrame( stressTensorarr: npt.NDArray[ np.float64 ], normal: npt.NDArray[ np.float64 ],
+                            tangent1: npt.NDArray[ np.float64 ],
+                            tangent2: npt.NDArray[ np.float64 ] ) -> dict[ str, Any ]:
+        """Rotate stress tensor to fault local coordinate system."""
+        # Verify orthonormality
+        assert np.abs( np.linalg.norm( tangent1 ) - 1.0 ) < 1e-10, f"T1 - {np.abs( np.linalg.norm( tangent1 ) - 1.0 )}"
+        assert np.abs( np.linalg.norm( tangent2 ) - 1.0 ) < 1e-10, f"T2 - {np.abs( np.linalg.norm( tangent2 ) - 1.0 )}"
+        assert np.abs( np.dot( normal, tangent1 ) ) < 1e-10
+        assert np.abs( np.dot( normal, tangent2 ) ) < 1e-10
+
+        # Rotation matrix: columns = local directions (n, t1, t2)
+        R = np.column_stack( ( normal, tangent1, tangent2 ) )
+
+        # Rotate tensor
+        stressLocal = R.T @ stressTensorarr @ R
+
+        # Traction on fault plane (normal = [1,0,0] in local frame)
+        tractionLocal = stressLocal @ np.array( [ 1.0, 0.0, 0.0 ] )
+
+        return {
+            'stressLocal': stressLocal,
+            'normalStress': tractionLocal[ 0 ],
+            'shearStress': np.sqrt( tractionLocal[ 1 ]**2 + tractionLocal[ 2 ]**2 ),
+            'shearStrike': tractionLocal[ 1 ],
+            'shearDip': tractionLocal[ 2 ]
+        }

geos-mesh/src/geos/mesh/io/vtkIO.py

@@ -5,6 +5,7 @@
 from enum import Enum
 from pathlib import Path
 from typing import Optional, Type, TypeAlias
+from xml.etree import ElementTree as ET
 from vtkmodules.vtkCommonDataModel import vtkPointSet, vtkUnstructuredGrid
 from vtkmodules.vtkIOCore import vtkWriter
 from vtkmodules.vtkIOLegacy import vtkDataReader, vtkUnstructuredGridWriter, vtkUnstructuredGridReader
@@ -266,3 +267,31 @@ def writeMesh( mesh: vtkPointSet, vtkOutput: VtkOutput, canOverwrite: bool = Fal
     except ( ValueError, RuntimeError ) as e:
         ioLogger.error( e )
         raise
+
+
+class PVDReader:
+    def __init__( self, filename ):
+        self.filename = filename
+        self.dir = Path( filename ).parent
+        self.datasets = {}
+        self._read()
+
+    def _read( self ):
+        tree = ET.parse( self.filename )
+        root = tree.getroot()
+        datasets = root[0].findall( 'DataSet' )
+
+        n: int = 0
+        for dataset in datasets:
+            timestep = float( dataset.attrib.get( 'timestep', 0 ) )
+            datasetFile = Path( dataset.attrib.get( 'file' ) )
+            # self.datasets.update( ( n, timestep ): datasetFile )
+            self.datasets[ n ] = ( timestep, datasetFile )
+            n += 1
+
+
+    def getDataSetAtTimeIndex( self, timeIndex: int):
+        return readMesh(  self.dir / self.datasets[ timeIndex ][ 1 ] )
+
+    def getAllTimestepsValues( self ) -> list[ float ]:
+        return list( [ value[0] for _, value in self.datasets.items() ] )

geos-mesh/src/geos/mesh/utils/arrayModifiers.py

@@ -898,6 +898,32 @@ def renameAttribute(
     return
 
 
+def updateAttribute( mesh: vtkDataSet, newValue: npt.NDArray[ Any ], attributeName: str, piece: Piece = Piece.CELLS, logger: Union[ Logger, None ] = None ) -> None:
+    """Update the value of an attribute. Creates the attribute if it is not already in the dataset.
+
+    Args:
+        mesh (vtkDataSet): Input mesh.
+        attributeName (str): Name of the attribute.
+        newValue (vtkDataArray):
+    """
+    # Check if an external logger is given.
+    if logger is None:
+        logger = getLogger( "updateAttribute", True )
+
+    if isAttributeInObject( mesh, attributeName, piece ):
+        if piece == Piece.CELLS:
+            data = mesh.GetCellData()
+        elif piece == Piece.POINTS:
+            data = mesh.GetPointData()
+        else:
+            raise ValueError( "Only point and cell data handled." )
+        data.RemoveArray( attributeName )
+
+    createAttribute( mesh, newValue, attributeName, piece=piece, logger=logger )
+
+    return
+
+
 def createCellCenterAttribute(
     mesh: Union[ vtkMultiBlockDataSet, vtkDataSet ],
     cellCenterAttributeName: str,

geos-mesh/src/geos/mesh/utils/genericHelpers.py

@@ -15,6 +15,15 @@
 from vtkmodules.vtkFiltersTexture import vtkTextureMapToPlane
 from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter
 from vtkmodules.vtkFiltersGeneral import vtkDataSetTriangleFilter
+from vtkmodules.util.numpy_support import numpy_to_vtkIdTypeArray
+from vtkmodules.vtkFiltersExtraction import vtkExtractSelection
+from vtkmodules.vtkFiltersGeometry import vtkGeometryFilter
+from vtkmodules.vtkFiltersVerdict import vtkCellSizeFilter
+from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkFieldData, vtkMultiBlockDataSet, vtkDataSet,
+                                            vtkCompositeDataSet, vtkDataObject, vtkPointData, vtkCellData, vtkPolyData,
+                                            vtkCell,  , vtkSelection )
+from typing import cast
+
 
 from geos.mesh.utils.multiblockHelpers import ( getBlockElementIndexesFlatten, getBlockFromFlatIndex )
 
@@ -515,6 +524,7 @@ def getLocalBasisVectors(
 
 def computeNormals(
     surface: vtkPolyData,
+    pointNormals: bool = False,
     logger: Union[ Logger, None ] = None,
 ) -> vtkPolyData:
     """Compute and set the normals of a given surface.
@@ -665,3 +675,41 @@ def computeSurfaceTextureCoordinates(
         vtkErrorLogger.error( captured.strip() )
 
     return textureFilter.GetOutput()
+
+
+def extractCellSelection( mesh: vtkUnstructuredGrid, ids: list[ int ]) -> vtkUnstructuredGrid:
+
+    selectionNode: vtkSelectionNode = vtkSelectionNode()
+    selectionNode.SetFieldType( vtkSelectionNode.CELL )
+    selectionNode.SetContentType( vtkSelectionNode.INDICES )
+    selectionNode.SetSelectionList( numpy_to_vtkIdTypeArray (np.asarray( ids ).astype( np.int64 ) ) )
+
+    selection: vtkSelection = vtkSelection()
+    selection.AddNode( selectionNode )
+
+    extractCells = vtkExtractSelection()
+    extractCells.SetInputData(0, mesh )
+    extractCells.SetInputData(1, selection )
+    extractCells.Update()
+
+    # TODO raiseError
+    return vtkUnstructuredGrid.SafeDownCast( extractCells.GetOutputDataObject( 0 ) )
+
+
+def extractSurface( mesh: vtkUnstructuredGrid ) -> vtkUnstructuredGrid:
+    geomFilter: vtkGeometryFilter = vtkGeometryFilter()
+    geomFilter.SetInputData( mesh )
+
+    geomFilter.Update()
+
+    return geomFilter.GetOutput()
+
+
+
+def computeCellVolumes( mesh: vtkUnstructuredGrid ) -> vtkUnstructuredGrid:
+    volFilter: vtkCellSizeFilter = vtkCellSizeFilter()
+    volFilter.SetInputData( mesh )
+    volFilter.SetComputeVolume( True )
+    volFilter.Update()
+
+    return volFilter.GetOutput()