feat: Adding a fault stability analysis workflow

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

@@ -0,0 +1,79 @@
+# 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.
+
+        Args:
+            arr ( npt.NDArray[np.float64]): Array to convert.
+
+        Returns:
+            npt.NDArray[np.float64]: 3x3 converted stress tensor.
+        """
+        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.
+
+        Args:
+            stressTensorArr (npt.NDArray[np.float64]): Stress tensor to rotate.
+            normal (npt.NDArray[np.float64]): Surface normal vectors.
+            tangent1 (npt.NDArray[np.float64]): Surface tangents vectors 1.
+            tangent2 (npt.NDArray[np.float64])): Surface tangents vectors 2.
+
+        Returns:
+            dict[str, Any]: Dictionary containing local stress, normal stress, shear stress and strike and shear dip.
+        """
+        # 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

@@ -1,16 +1,19 @@
 # SPDX-License-Identifier: Apache-2.0
 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
 # SPDX-FileContributor: Alexandre Benedicto
+import os
 from dataclasses import dataclass
 from enum import Enum
 from pathlib import Path
-from typing import Optional, Type, TypeAlias
-from vtkmodules.vtkCommonDataModel import vtkPointSet, vtkUnstructuredGrid
+from typing import Optional, Type, TypeAlias, Self
+from xml.etree import ElementTree as ET
+from vtkmodules.vtkCommonDataModel import vtkPointSet, vtkUnstructuredGrid, vtkDataSet
 from vtkmodules.vtkIOCore import vtkWriter
 from vtkmodules.vtkIOLegacy import vtkDataReader, vtkUnstructuredGridWriter, vtkUnstructuredGridReader
 from vtkmodules.vtkIOXML import ( vtkXMLGenericDataObjectReader, vtkXMLUnstructuredGridWriter, vtkXMLWriter,
                                   vtkXMLStructuredGridWriter )
-from geos.utils.Logger import getLogger
+
+from geos.utils.Logger import ( getLogger )
 
 __doc__ = """
 Input and Output methods for various VTK mesh formats.
@@ -266,3 +269,64 @@ def writeMesh( mesh: vtkPointSet, vtkOutput: VtkOutput, canOverwrite: bool = Fal
     except ( ValueError, RuntimeError ) as e:
         ioLogger.error( e )
         raise
+
+
+class PVDReader:
+
+    def __init__( self: Self, filename: str ) -> None:
+        """PVD Reader class.
+
+        Args:
+            filename (str): PVD filename
+        """
+        self.filename = filename
+        self.dir = Path( filename ).parent
+        self.datasets = {}
+        self._read()
+
+    def _read( self ) -> None:
+        tree = ET.parse( self.filename )
+        root = tree.getroot()
+        datasets = root[ 0 ].findall( 'DataSet' )
+
+        for n, dataset in enumerate( datasets ):
+            timestep = float( dataset.attrib.get( 'timestep', 0 ) )
+            datasetFile = Path( dataset.attrib.get( 'file' ) )
+            self.datasets[ n ] = ( timestep, datasetFile )
+
+    def getDataSetAtTimeIndex( self: Self, timeIndex: int ) -> vtkDataSet:
+        """Get the dataset corresponding to requested time index.
+
+        Args:
+            timeIndex (int): Time index
+
+        Returns:
+            vtkDataSet: Dataset
+        """
+        return readMesh( self.dir / self.datasets[ timeIndex ][ 1 ] )
+
+    def getAllTimestepsValues( self: Self ) -> list[ float ]:
+        """Get the list of all timesteps values from the PVD.
+
+        Returns:
+            list[float]: List of timesteps values.
+        """
+        return [ value[ 0 ] for _, value in self.datasets.items() ]
+
+
+def createPVD( outputDir: str, outputFiles: list[ tuple[ int, str ] ] ) -> None:
+    """Create PVD collection file.
+
+    Args:
+        outputDir (str): Output directory
+        outputFiles (list[tuple[int, str]]): List containing all the filenames of the PVD files
+    """
+    pvdPath = os.path.join( outputDir, 'fault_analysis.pvd' )
+    with open( pvdPath, 'w' ) as f:
+        f.write( '<VTKFile type="Collection" version="0.1">\n' )
+        f.write( '  <Collection>\n' )
+        for t, fname in outputFiles:
+            f.write( f'    <DataSet timestep="{t}" file="{fname}"/>\n' )
+        f.write( '  </Collection>\n' )
+        f.write( '</VTKFile>\n' )
+    print( f"\n✅ PVD created: {pvdPath}" )

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

@@ -747,6 +747,40 @@ 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.
+        newValue (npt.NDArray[Any]): New value for the attribute.
+        attributeName (str): Name of the attribute.
+        piece (Piece): The piece of the attribute.
+        logger (Union[Logger, None], optional): A logger to manage the output messages.
+            Defaults to None, an internal logger is used.
+    """
+    # Check if an external logger is given.
+    if logger is None:
+        logger = getLogger( "updateAttribute", True )
+
+    if isAttributeInObject( mesh, attributeName, piece ):
+        data: Union[ vtkPointData, vtkCellData ]
+        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,
@@ -833,7 +867,6 @@ def transferPointDataToCellData(
 
     Returns:
         vtkPointSet: Output mesh where point data were transferred to cells.
-
     """
     if logger is None:
         logger = getLogger( "transferPointDataToCellData", True )

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

@@ -10,18 +10,21 @@
 from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints, reference, vtkDataArray, vtkLogger, vtkFloatArray
 from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkMultiBlockDataSet, vtkPolyData, vtkDataSet,
                                             vtkDataObject, vtkPlane, vtkCellTypes, vtkIncrementalOctreePointLocator,
-                                            VTK_TRIANGLE )
+                                            VTK_TRIANGLE, vtkSelection, vtkSelectionNode )
 from vtkmodules.vtkFiltersCore import ( vtk3DLinearGridPlaneCutter, vtkPolyDataNormals, vtkPolyDataTangents )
 from vtkmodules.vtkFiltersTexture import vtkTextureMapToPlane
-from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter
+from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter, vtkGeometryFilter
 from vtkmodules.vtkFiltersGeneral import vtkDataSetTriangleFilter
+from vtkmodules.util.numpy_support import numpy_to_vtkIdTypeArray
+from vtkmodules.vtkFiltersExtraction import vtkExtractSelection
+from vtkmodules.vtkFiltersVerdict import vtkCellSizeFilter
 
 from geos.mesh.utils.multiblockHelpers import ( getBlockElementIndexesFlatten, getBlockFromFlatIndex )
 
 from geos.utils.algebraFunctions import ( getAttributeMatrixFromVector, getAttributeVectorFromMatrix )
 from geos.utils.geometryFunctions import ( getChangeOfBasisMatrix, CANONICAL_BASIS_3D )
 from geos.utils.Logger import ( getLogger, Logger, VTKCaptureLog, RegexExceptionFilter )
-from geos.utils.Errors import VTKError
+from geos.utils.Errors import ( VTKError )
 
 __doc__ = """
 Generic VTK utilities.
@@ -515,12 +518,14 @@ def getLocalBasisVectors(
 
 def computeNormals(
     surface: vtkPolyData,
+    pointNormals: bool = False,
     logger: Union[ Logger, None ] = None,
 ) -> vtkPolyData:
     """Compute and set the normals of a given surface.
 
     Args:
         surface (vtkPolyData): The input surface.
+        pointNormals (bool): Flag to compute point normals. Defaults is False.
         logger (Union[Logger, None], optional): A logger to manage the output messages.
             Defaults to None, an internal logger is used.
 
@@ -665,3 +670,63 @@ def computeSurfaceTextureCoordinates(
         vtkErrorLogger.error( captured.strip() )
 
     return textureFilter.GetOutput()
+
+
+def extractCellSelection( mesh: vtkUnstructuredGrid, ids: list[ int ] ) -> vtkUnstructuredGrid:
+    """Extract cell selection from list of cell Ids.
+
+    Args:
+        mesh (vtkUnstructuredGrid): Initial mesh.
+        ids (list[int]): Cell Ids of the cells to extract.
+
+    Returns:
+        vtkUnstructuredGrid: Subset of the input mesh.
+    """
+    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()
+
+    return vtkUnstructuredGrid.SafeDownCast( extractCells.GetOutputDataObject( 0 ) )
+
+
+def extractSurface( mesh: vtkUnstructuredGrid ) -> vtkUnstructuredGrid:
+    """Extract surface from an input mesh.
+
+    Args:
+        mesh: Input mesh
+
+    Returns:
+        vtkUnstructuredGrid: Surface of the input mesh.
+    """
+    geomFilter: vtkGeometryFilter = vtkGeometryFilter()
+    geomFilter.SetInputData( mesh )
+
+    geomFilter.Update()
+
+    return geomFilter.GetOutput()
+
+
+def computeCellVolumes( mesh: vtkUnstructuredGrid ) -> vtkUnstructuredGrid:
+    """Compute and return the cell volumes.
+
+    Args:
+        mesh: Input mesh.
+
+    Returns:
+        vtkUnstructuredGrid: Mesh with volume attribute.
+    """
+    volFilter: vtkCellSizeFilter = vtkCellSizeFilter()
+    volFilter.SetInputData( mesh )
+    volFilter.SetComputeVolume( True )
+    volFilter.Update()
+
+    return volFilter.GetOutput()

geos-processing/pyproject.toml

@@ -35,6 +35,7 @@ dependencies = [
     "vtk >= 9.3, < 9.6",
     "numpy >= 2.2",
     "typing_extensions >= 4.12",
+    "scipy",
 ]