feat: Adding a fault stability analysis workflow (#232)

Author: paloma-martinez

docs/geos_processing_docs/geos-processing.rst

@@ -12,3 +12,5 @@ The `geos-processing` package contains different tools to prepare, post-process
     pre_processing.rst
 
     post_processing.rst
+
+    tools.rst

docs/geos_processing_docs/post_processing.rst

@@ -20,6 +20,14 @@ GEOS computes many outputs including flow and geomechanic properties if coupled
     - friction angle: 10°
 
 
+FaultStabilityAnalysis
+-----------------------------
+
+.. automodule:: geos.processing.post_processing.FaultStabilityAnalysis
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 
 GeomechanicsCalculator
 --------------------------------

docs/geos_processing_docs/tools.rst

@@ -0,0 +1,50 @@
+Utilities classes for processing filters
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The `tools` folder contains utilities classes that are used in some of the processing filters.
+
+FaultGeometry
+----------------------
+.. automodule:: geos.processing.tools.FaultGeometry
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+FaultVisualizer
+----------------------
+.. automodule:: geos.processing.tools.FaultVisualizer
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+MohrCoulomb
+----------------------
+.. automodule:: geos.processing.tools.MohrCoulomb
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+ProfileExtractor
+----------------------
+.. automodule:: geos.processing.tools.ProfileExtractor
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+SensitivityAnalyzer
+----------------------
+.. automodule:: geos.processing.tools.SensitivityAnalyzer
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+StressProjector
+----------------------
+.. automodule:: geos.processing.tools.StressProjector
+    :members:
+    :undoc-members:
+    :show-inheritance:

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
+        if np.abs( np.linalg.norm( tangent1 ) - 1.0 ) >= 1e-10 or np.abs( np.linalg.norm( tangent2 ) - 1.0 ) >= 1e-10:
+            raise ValueError( "Tangents expected to be normalized." )
+        if np.abs( np.dot( normal, tangent1 ) ) >= 1e-10 or np.abs( np.dot( normal, tangent2 ) ) >= 1e-10:
+            raise ValueError( "Tangents and Normals expected to be orthogonal." )
+
+        # 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,20 @@
 # SPDX-License-Identifier: Apache-2.0
 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
 # SPDX-FileContributor: Alexandre Benedicto
+import logging
 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_extensions import Self, Union
+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, Logger )
 
 __doc__ = """
 Input and Output methods for various VTK mesh formats.
@@ -104,19 +108,28 @@ def _readData( filepath: str, readerClass: VtkReaderClass ) -> Optional[ vtkPoin
     return output
 
 
-def _writeData( mesh: vtkPointSet, writerClass: VtkWriterClass, output: str, isBinary: bool ) -> int:
+def _writeData( mesh: vtkPointSet,
+                writerClass: VtkWriterClass,
+                output: str,
+                isBinary: bool,
+                logger: Union[ Logger, None ] = None ) -> int:
     """Generic helper to write a VTK file using a specific writer class.
 
     Args:
         mesh (vtkPointSet): The grid data to write.
         writerClass (VtkWriterClass): The VTK writer class to use.
         output (str): The output file path.
         isBinary (bool): Whether to write the file in binary mode (True) or ASCII (False).
+        logger (Union[Logger, None], optional): A logger to manage the output messages.
+                    Defaults to None, an internal logger is used.
 
     Returns:
         int: 1 if success, 0 otherwise.
     """
-    ioLogger.info( f"Writing mesh to '{output}' using {writerClass.__name__}..." )
+    if logger is None:
+        logger = ioLogger
+
+    logger.info( f"Writing mesh to '{output}' using {writerClass.__name__}..." )
     writer = writerClass()
     writer.SetFileName( output )
     writer.SetInputData( mesh )
@@ -125,10 +138,10 @@ def _writeData( mesh: vtkPointSet, writerClass: VtkWriterClass, output: str, isB
     if isinstance( writer, vtkXMLWriter ):
         if isBinary:
             writer.SetDataModeToBinary()
-            ioLogger.info( "Data mode set to Binary." )
+            logger.info( "Data mode set to Binary." )
         else:
             writer.SetDataModeToAscii()
-            ioLogger.info( "Data mode set to ASCII." )
+            logger.info( "Data mode set to ASCII." )
 
     return writer.Write()
 
@@ -220,7 +233,10 @@ def readUnstructuredGrid( filepath: str ) -> vtkUnstructuredGrid:
     return mesh
 
 
-def writeMesh( mesh: vtkPointSet, vtkOutput: VtkOutput, canOverwrite: bool = False ) -> int:
+def writeMesh( mesh: vtkPointSet,
+               vtkOutput: VtkOutput,
+               canOverwrite: bool = False,
+               logger: Union[ Logger, None ] = None ) -> int:
     """
     Writes a vtkPointSet to a file.
 
@@ -231,6 +247,8 @@ def writeMesh( mesh: vtkPointSet, vtkOutput: VtkOutput, canOverwrite: bool = Fal
         vtkOutput (VtkOutput): Configuration for the output file.
         canOverwrite (bool, optional): If False, raises an error if the file
                                       already exists. Defaults to False.
+        logger (Union[Logger, None], optional): A logger to manage the output messages.
+                    Defaults to None, an internal logger is used.
 
     Raises:
         FileExistsError: If the output file exists and `canOverwrite` is False.
@@ -240,6 +258,9 @@ def writeMesh( mesh: vtkPointSet, vtkOutput: VtkOutput, canOverwrite: bool = Fal
     Returns:
         int: Returns 1 on success, consistent with the VTK writer's return code.
     """
+    if logger is None:
+        logger = ioLogger
+
     outputPath = Path( vtkOutput.output )
     if outputPath.exists() and not canOverwrite:
         raise FileExistsError( f"File '{outputPath}' already exists. Set canOverwrite=True to replace it." )
@@ -256,13 +277,96 @@ def writeMesh( mesh: vtkPointSet, vtkOutput: VtkOutput, canOverwrite: bool = Fal
         if not writerClass:
             raise ValueError( f"Writing to extension '{outputPath.suffix}' is not supported." )
 
-        successCode = _writeData( mesh, writerClass, str( outputPath ), vtkOutput.isDataModeBinary )
+        successCode = _writeData( mesh, writerClass, str( outputPath ), vtkOutput.isDataModeBinary, logger )
         if not successCode:
             raise RuntimeError( f"VTK writer failed to write file '{outputPath}'." )
 
-        ioLogger.info( f"Successfully wrote mesh to '{outputPath}'." )
+        logger.info( f"Successfully wrote mesh to '{outputPath}'." )
         return successCode
 
     except ( ValueError, RuntimeError ) as e:
-        ioLogger.error( e )
+        logger.error( e )
         raise
+
+
+class PVDReader:
+
+    def __init__( self: Self, filename: str, logger: Union[ Logger, None ] = None ) -> None:
+        """PVD Reader class.
+
+        Args:
+            filename (str): PVD filename with full path.
+            logger (Union[Logger, None], optional): A logger to manage the output messages.
+                    Defaults to None, an internal logger is used.
+        """
+        self.logger: Logger
+        if logger is None:
+            self.logger = getLogger( "PVD Reader", True )
+        else:
+            self.logger = logging.getLogger( f"{logger.name}" )
+            self.logger.setLevel( logging.INFO )
+            self.logger.propagate = False
+
+        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 )
+
+        self.logger.info( "All filenames from PVD file have been read." )
+
+    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: Path,
+               pvdFilename: str,
+               outputFiles: list[ tuple[ int, str ] ],
+               logger: Union[ Logger, None ] = None ) -> None:
+    """Create PVD collection file.
+
+    Args:
+        outputDir (Path): Output directory
+        pvdFilename (str): Output PVD filename
+        outputFiles (list[tuple[int, str]]): List containing all the filenames of the PVD files
+        logger (Union[Logger, None], optional): A logger to manage the output messages.
+                    Defaults to None, an internal logger is used.
+    """
+    if logger is None:
+        logger = getLogger( "createPVD", True )
+
+    pvdPath = outputDir / pvdFilename
+    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' )
+
+    logger.info( f"PVD created: {pvdPath}." )