diff --git a/geos-geomechanics/src/geos/geomechanics/model/StressTensor.py b/geos-geomechanics/src/geos/geomechanics/model/StressTensor.py new file mode 100644 index 000000000..619fc5420 --- /dev/null +++ b/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 ] + } diff --git a/geos-mesh/src/geos/mesh/io/vtkIO.py b/geos-mesh/src/geos/mesh/io/vtkIO.py index 55efbc9f7..1b721217a 100644 --- a/geos-mesh/src/geos/mesh/io/vtkIO.py +++ b/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() ] ) \ No newline at end of file diff --git a/geos-mesh/src/geos/mesh/utils/arrayModifiers.py b/geos-mesh/src/geos/mesh/utils/arrayModifiers.py index 6bc763b32..22672f9b1 100644 --- a/geos-mesh/src/geos/mesh/utils/arrayModifiers.py +++ b/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, diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py index de11a4023..bbc422c45 100644 --- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py +++ b/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() \ No newline at end of file diff --git a/geos-processing/src/geos/processing/post_processing/FaultGeometry.py b/geos-processing/src/geos/processing/post_processing/FaultGeometry.py new file mode 100644 index 000000000..1f89be582 --- /dev/null +++ b/geos-processing/src/geos/processing/post_processing/FaultGeometry.py @@ -0,0 +1,865 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2026 TotalEnergies. +# SPDX-FileContributor: Nicolas Pillardou, Paloma Martinez +# ============================================================================ +# FAULT GEOMETRY +# ============================================================================ +import sys +import pyvista as pv +import numpy as np +from pathlib import Path +from typing_extensions import Self, Any +from vtkmodules.vtkCommonDataModel import vtkCellLocator, vtkMultiBlockDataSet +from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridWriter + +import numpy.typing as npt +from scipy.spatial import cKDTree + +from vtkmodules.util.numpy_support import vtk_to_numpy, numpy_to_vtk + +from geos.mesh.utils.arrayModifiers import createAttribute +from geos.mesh.utils.arrayHelpers import ( getArrayInObject, computeCellCenterCoordinates ) +from geos.mesh.utils.multiblockModifiers import ( mergeBlocks ) +from geos.mesh.utils.genericHelpers import ( extractCellSelection, extractSurface, computeNormals, getNormalVectors, computeCellVolumes ) +from geos.utils.pieceEnum import Piece + +from geos.mesh.io.vtkIO import writeMesh, VtkOutput + +__doc__=""" + + +""" + + +class FaultGeometry: + """Handles fault surface extraction and normal computation with optimizations.""" + + # ------------------------------------------------------------------- + def __init__( self: Self, mesh: pv.DataSet, faultValues: list[ int ], faultAttribute: str, + volumeMesh: pv.DataSet, outputDir: str = "." ) -> None: + """Initialize fault geometry with pre-computed topology. + + Args: + mesh (pv.DataSet): Input mesh + faultValues (list[int]): Config.FAULT_VALUES + faultAttribute (str): Config.FAULT_ATTRIBUTES + volumeMesh (pv.DataSet): PVD mesh + """ + self.mesh = mesh + self.faultValues = faultValues + self.faultAttribute = faultAttribute + self.volumeMesh = volumeMesh + + # These will be computed once + self.faultSurface = None + self.surfaces = None + self.adjacencyMapping = None + self.contributingCells = None + self.contributingCellsPlus = None + self.contributingCellsMinus = None + + # NEW: Pre-computed geometric properties + self.volumeCellVolumes = None # Volume of each cell + self.volumeCenters = None # Center coordinates + self.distanceToFault = None # Distance from each volume cell to nearest fault + self.faultTree = None # KDTree for fault surface + + # Config + # self.config = config + self.outputDir = Path( outputDir ) + self.outputDir.mkdir( parents=True, exist_ok=True ) + + # ------------------------------------------------------------------- + def initialize( self: Self, + scaleFactor: float = 50.0, + processFaultsSeparately: bool = True, + showPlot: bool = True, + zscale: float = 1.0, + showContributionViz: bool = True, + saveContributionCells:bool = True ) -> tuple[ pv.DataSet, dict[ int, pv.DataSet ] ]: + """One-time initialization: compute normals, adjacency topology, and geometric properties.""" + # Extract and compute normals + self.faultSurface, self.surfaces = self._extractAndComputeNormals( showPlot=showPlot, + scaleFactor=scaleFactor, + zScale=zscale ) + + # Pre-compute adjacency mapping + print( "\n🔍 Pre-computing volume-fault adjacency topology" ) + print( " Method: Face-sharing (adaptive epsilon)" ) + + self.adjacencyMapping = self._buildAdjacencyMappingFaceSharing( + processFaultsSeparately=processFaultsSeparately ) + + # Mark and optionally save contributing cells + self._markContributingCells( saveContributionCells ) + + # NEW: Pre-compute geometric properties + self._precomputeGeometricProperties() + + nMapped = len( self.adjacencyMapping ) + nWithBoth = sum( 1 for m in self.adjacencyMapping.values() + if len( m[ 'plus' ] ) > 0 and len( m[ 'minus' ] ) > 0 ) + + print( "\n✅ Adjacency topology computed:" ) + print( f" - {nMapped}/{self.faultSurface.GetNumberOfCells()} fault cells mapped" ) + print( f" - {nWithBoth} cells have neighbors on both sides" ) + + # Visualize contributions if requested + # if showContributionViz: + # self._visualizeContributions() + + return self.faultSurface, self.adjacencyMapping + + # ------------------------------------------------------------------- + def _markContributingCells( self: Self, saveContributionCells: bool = True ) -> None: + """Mark volume cells that contribute to fault stress projection.""" + print( "\n📦 Marking contributing volume cells..." ) + + nVolume = self.volumeMesh.GetNumberOfCells() + + # Collect contributing cells by side + allPlus = set() + allMinus = set() + + for _faultIdx, neighbors in self.adjacencyMapping.items(): + allPlus.update( neighbors[ 'plus' ] ) + allMinus.update( neighbors[ 'minus' ] ) + + # Create classification array + contributionSide = np.zeros( nVolume, dtype=int ) + + for idx in allPlus: + if 0 <= idx < nVolume: + contributionSide[ idx ] += 1 + + for idx in allMinus: + if 0 <= idx < nVolume: + contributionSide[ idx ] += 2 + + # Add classification to volume mesh + contribMask = contributionSide > 0 + contribMask = contribMask.astype( int ) + + createAttribute( self.volumeMesh, contributionSide, "contributionSide" ) + createAttribute( self.volumeMesh, contribMask, "contributionToFaults" ) + + # Extract subsets + maskAll = np.where( contribMask )[0] + maskPlus = np.where( ( contributionSide == 1 ) | ( contributionSide == 3 ) )[0] + maskMinus = np.where( ( contributionSide == 2 ) | ( contributionSide == 3 ) )[0] + + self.contributingCells = extractCellSelection( self.volumeMesh, maskAll ) + self.contributingCellsPlus = extractCellSelection( self.volumeMesh, maskPlus ) + self.contributingCellsMinus = extractCellSelection( self.volumeMesh, maskMinus ) + + # Statistics + nContrib = np.sum( maskAll ) + nPlus = np.sum( contributionSide == 1 ) + nMinus = np.sum( contributionSide == 2 ) + nBoth = np.sum( contributionSide == 3 ) + pctContrib = nContrib / nVolume * 100 + + print( f" ✅ Total contributing: {nContrib}/{nVolume} ({pctContrib:.1f}%)" ) + print( f" Plus side only: {nPlus} cells" ) + print( f" Minus side only: {nMinus} cells" ) + print( f" Both sides: {nBoth} cells" ) + + # Save to files if requested + if saveContributionCells: + self._saveContributingCells() + + # ------------------------------------------------------------------- + def _saveContributingCells( self: Self ) -> None: + """Save contributing volume cells to VTU files. + + Saves three files: all, plus side, minus side. + """ + # Create output directory if it doesn't exist + outputDir = self.outputDir + + # Save all contributing cells + filenameAll = outputDir / "contributing_cells_all.vtu" + + writeMesh( mesh=self.contributingCells, vtkOutput=VtkOutput(filenameAll), canOverwrite=True ) + # self.contributingCells.save( str( filenameAll ) ) + print( f"\n 💾 All contributing cells saved: {filenameAll}" ) + print( f" ({self.contributingCells.GetNumberOfCells()} cells, {self.contributingCells.GetNumberOfPoints} points)" ) + + # Save plus side + filenamePlus = outputDir / "contributingCellsPlus.vtu" + # self.contributingCellsPlus.save(str(filenamePlus)) + # print(f" 💾 Plus side cells saved: {filenamePlus}") + writeMesh( mesh=self.contributingCellsPlus, vtkOutput=VtkOutput(filenamePlus), canOverwrite=True ) + print( f" ({self.contributingCellsPlus.GetNumberOfCells()} cells, {self.contributingCellsPlus.GetNumberOfPoints} points)" ) + + # Save minus side + filenameMinus = outputDir / "contributingCellsMinus.vtu" + # self.contributingCellsMinus.save(str(filenameMinus)) + # print(f" 💾 Minus side cells saved: {filenameMinus}") + writeMesh( mesh=self.contributingCellsMinus, vtkOutput=VtkOutput(filenameMinus), canOverwrite=True ) + print( f" ({self.contributingCellsMinus.GetNumberOfCells()} cells, {self.contributingCellsMinus.GetNumberOfPoints} points)" ) + + # ------------------------------------------------------------------- + def getContributingCells( self: Self, side: str = 'all' ) -> pv.UnstructuredGrid: + """Get the extracted contributing cells. + + Parameters: + side: 'all', 'plus', or 'minus' + + Returns: + pyvista.UnstructuredGrid: Contributing volume cells + """ + if self.contributingCells is None: + raise ValueError( "Contributing cells not yet computed. Call initialize() first." ) + + if side == 'all': + return self.contributingCells + elif side == 'plus': + return self.contributingCellsPlus + elif side == 'minus': + return self.contributingCellsMinus + else: + raise ValueError( f"Invalid side '{side}'. Must be 'all', 'plus', or 'minus'." ) + + # ------------------------------------------------------------------- + def getGeometricProperties( self: Self ) -> dict[ str, Any ]: + """Get pre-computed geometric properties. + + Returns: + ------- + dict with keys: + - 'volumes': ndarray of cell volumes + - 'centers': ndarray of cell centers (nCells, 3) + - 'distances': ndarray of distances to nearest fault cell + - 'faultTree': KDTree for fault surface + """ + if self.volumeCellVolumes is None: + raise ValueError( "Geometric properties not computed. Call initialize() first." ) + + return { + 'volumes': self.volumeCellVolumes, + 'centers': self.volumeCenters, + 'distances': self.distanceToFault, + 'faultTree': self.faultTree + } + + # ------------------------------------------------------------------- + def _precomputeGeometricProperties( self: Self ) -> None: # TODO + """Pre-compute geometric properties of volume mesh for efficient stress projection. + + Computes: + - Cell volumes (for volume-weighted averaging) + - Cell centers (for distance calculations) + - Distance from each volume cell to nearest fault cell + - KDTree for fault surface + """ + print( "\n📐 Pre-computing geometric properties..." ) + + nVolume = self.volumeMesh.GetNumberOfCells() + + # 1. Compute volume centers + print( " Computing cell centers..." ) + self.volumeCenters = vtk_to_numpy( computeCellCenterCoordinates( self.volumeMesh ) ) + + # 2. Compute cell volumes + print( " Computing cell volumes..." ) + volumeWithSizes = computeCellVolumes( self.volumeMesh ) + self.volumeCellVolumes = getArrayInObject( volumeWithSizes, 'Volume', Piece.CELLS ) + + print( f" Volume range: [{np.min(self.volumeCellVolumes):.1e}, " + f"{np.max(self.volumeCellVolumes):.1e}] m³" ) + + # 3. Build KDTree for fault surface (for fast distance queries) + print( " Building KDTree for fault surface..." ) + + faultCenters = computeCellCenterCoordinates( self.faultSurface ) + self.faultTree = cKDTree( faultCenters ) + + # 4. Compute distance from each volume cell to nearest fault cell + print( " Computing distances to fault..." ) + self.distanceToFault = np.zeros( nVolume ) + + # Vectorized query for all points at once (much faster) + distances, _ = self.faultTree.query( self.volumeCenters ) + self.distanceToFault = distances + + print( f" Distance range: [{np.min(self.distanceToFault):.1f}, " + f"{np.max(self.distanceToFault):.1f}] m" ) + + # 5. Add these properties to volume mesh for reference + createAttribute ( self.volumeMesh, self.volumeCellVolumes, 'cellVolume', Piece.CELLS ) + createAttribute ( self.volumeMesh, self.distanceToFault, 'distanceToFault', Piece.CELLS ) + + print( " ✅ Geometric properties computed and cached" ) + + # ------------------------------------------------------------------- + def _buildAdjacencyMappingFaceSharing( self: Self, + processFaultsSeparately: bool = True ) -> dict[ int, pv.DataSet ]: + """Build adjacency for cells sharing faces with fault. + + Uses adaptive epsilon optimization. + """ + faultIds = np.unique( getArrayInObject ( self.faultSurface, self.faultAttribute, Piece.CELLS ) ) + nFaults = len( faultIds ) + print( f" 📋 Processing {nFaults} separate faults: {faultIds}" ) + + allMappings = {} + + for faultId in faultIds: + mask = getArrayInObject( self.faultSurface, self.faultAttribute, Piece.CELLS ) == faultId + indices = np.where( mask )[ 0 ] + singleFault = extractCellSelection( self.faultSurface, indices ) + + print( f" 🔧 Mapping Fault {faultId}..." ) + + # Build face-sharing mapping with adaptive epsilon + localMapping = self._findFaceSharingCells( singleFault ) + + # Remap local indices to global fault indices + for localIdx, neighbors in localMapping.items(): + globalIdx = indices[ localIdx ] + allMappings[ globalIdx ] = neighbors + + return allMappings + + # ------------------------------------------------------------------- + def _findFaceSharingCells( self: Self, faultSurface ) -> pv.DataSet: + """Find volume cells that share a FACE with fault cells. + + Uses FindCell with adaptive epsilon to maximize cells with both neighbors + """ + volMesh = self.volumeMesh + volCenters = vtk_to_numpy( computeCellCenterCoordinates( volMesh ) ) + faultNormals = vtk_to_numpy( faultSurface.GetCellData().GetNormals() ) + faultCenters = vtk_to_numpy( computeCellCenterCoordinates( faultSurface ) ) + + # Determine base epsilon based on mesh size + volBounds = volMesh.bounds + typicalSize = np.mean( [ + volBounds[ 1 ] - volBounds[ 0 ], volBounds[ 3 ] - volBounds[ 2 ], volBounds[ 5 ] - volBounds[ 4 ] + ] ) / 100.0 + + # Build VTK cell locator (once) + locator = vtkCellLocator() + locator.SetDataSet( volMesh ) + locator.BuildLocator() + + # Try multiple epsilon values and keep the best + epsilonCandidates = [ + typicalSize * 0.005, typicalSize * 0.01, typicalSize * 0.05, typicalSize * 0.1, typicalSize * 0.2, + typicalSize * 0.5, typicalSize * 1.0 + ] + + print( f" Testing {len(epsilonCandidates)} epsilon values..." ) + + bestEpsilon = None + bestMapping = None + bestScore = -1 + bestStats = None + + for epsilon in epsilonCandidates: + # Test this epsilon + mapping, stats = self._testEpsilon( faultSurface, locator, epsilon, faultCenters, faultNormals, volCenters ) + + # Score = percentage with both sides + penalty for no neighbors + score = stats[ 'pctBoth' ] - 2.0 * stats[ 'pctNone' ] + + print( f" ε={epsilon:.3f}m → Both: {stats['pctBoth']:.1f}%, " + f"One: {stats['pctOne']:.1f}%, None: {stats['pctNone']:.1f}%, " + f"Avg: {stats['avgNeighbors']:.2f} (score: {score:.1f})" ) + + if score > bestScore: + bestScore = score + bestEpsilon = epsilon + bestMapping = mapping + bestStats = stats + + print( f"\n ✅ Best epsilon: {bestEpsilon:.6f}m" ) + print( " ✅ Face-sharing mapping completed:" ) + print( f" Both sides: {bestStats['nBoth']} ({bestStats['pctBoth']:.1f}%)" ) + print( f" One side: {bestStats['nOne']} ({bestStats['pctOne']:.1f}%)" ) + print( f" No neighbors: {bestStats['nNone']} ({bestStats['pctNone']:.1f}%)" ) + print( f" Average neighbors per fault cell: {bestStats['avgNeighbors']:.2f}" ) + + return bestMapping + + # ------------------------------------------------------------------- + def _testEpsilon( self: Self, faultSurface: pv.DataSet, locator: vtkCellLocator, epsilon: list[ float ], + faultCenters, faultNormals: npt.NDArray[ np.float64 ], volCenters ): + """Test a specific epsilon value and return mapping + statistics. + + Statistics include: + - 'nBoth': nFoundBoth, + - 'nOne': nFoundOne, + - 'nNone': nFoundNone, + - 'pctBoth': nFoundBoth / nCells * 100, + - 'pctOne': nFoundOne / nCells * 100, + - 'pctNone': nFoundNone / nCells * 100, + - 'avgNeighbors': avgNeighbors + """ + mapping = {} + nFoundBoth = 0 + nFoundOne = 0 + nFoundNone = 0 + totalNeighbors = 0 + + for fid in range( faultSurface.GetNumberOfCells() ): + fcenter = faultCenters[ fid ] + fnormal = faultNormals[ fid ] + + plusCells = [] + minusCells = [] + + # Search on PLUS side + pointPlus = fcenter + epsilon * fnormal + cellIdPlus = locator.FindCell( pointPlus ) + # print( cellIdPlus ) + if cellIdPlus >= 0: + plusCells.append( cellIdPlus ) + + # Search on MINUS side + pointMinus = fcenter - epsilon * fnormal + cellIdMinus = locator.FindCell( pointMinus ) + if cellIdMinus >= 0: + minusCells.append( cellIdMinus ) + + mapping[ fid ] = { "plus": plusCells, "minus": minusCells } + + # Statistics + nNeighbors = len( plusCells ) + len( minusCells ) + totalNeighbors += nNeighbors + + if len( plusCells ) > 0 and len( minusCells ) > 0: + nFoundBoth += 1 + elif len( plusCells ) > 0 or len( minusCells ) > 0: + nFoundOne += 1 + else: + nFoundNone += 1 + + nCells = faultSurface.GetNumberOfCells() + avgNeighbors = totalNeighbors / nCells if nCells > 0 else 0 + + stats = { + 'nBoth': nFoundBoth, + 'nOne': nFoundOne, + 'nNone': nFoundNone, + 'pctBoth': nFoundBoth / nCells * 100, + 'pctOne': nFoundOne / nCells * 100, + 'pctNone': nFoundNone / nCells * 100, + 'avgNeighbors': avgNeighbors + } + + return mapping, stats + + # ------------------------------------------------------------------- + def _visualizeContributions( self: Self, zscale: float = 1.0, showPlots:bool = True ) -> None: + """Unified visualization of volume contributions to fault surfaces. + + 4-panel view combining full context, side classification, clip, and slice. + """ + print( "\n📊 Creating contribution visualization..." ) + + # Create plotter with 4 subplots + plotter = pv.Plotter( shape=( 2, 2 ), window_size=[ 1800, 1400 ] ) + + # ========== PLOT 1: Full context (top-left) ========== + plotter.subplot( 0, 0 ) + plotter.add_text( "Full Context - Volume & Fault", font_size=14, position='upper_edge' ) + + # All volume (transparent) + plotter.add_mesh( self.mesh, color='lightgray', opacity=0.05, show_edges=False, label='Volume' ) + + # Fault surface (red) + plotter.add_mesh( self.faultSurface, color='red', opacity=1, show_edges=True, label='Fault Surface' ) + + plotter.add_legend( loc="upper left" ) + plotter.add_axes() + plotter.set_scale( zscale=zscale ) + + # ========== PLOT 2: Contributing cells by side (top-right) ========== + plotter.subplot( 0, 1 ) + plotter.add_text( "Contributing Cells", font_size=14, position='upper_edge' ) + + if 'contributionSide' in self.volumeMesh.cell_data: + # Plus side (blue) + if self.contributingCellsPlus.n_cells > 0: + plotter.add_mesh( self.contributingCellsPlus, + color='dodgerblue', + opacity=1.0, + show_edges=True, + label=f'Plus side ({self.contributingCellsPlus.n_cells} cells)' ) + + # Minus side (orange) + if self.contributingCellsMinus.n_cells > 0: + plotter.add_mesh( self.contributingCellsMinus, + color='darkorange', + opacity=1.0, + show_edges=True, + label=f'Minus side ({self.contributingCellsMinus.n_cells} cells)' ) + + # Fault surface for reference + plotter.add_mesh( self.faultSurface, color='red', opacity=1.0, show_edges=True, label='Fault' ) + + plotter.add_legend( loc='upper right' ) + plotter.add_axes() + plotter.set_scale( zscale=zscale ) + + # ========== PLOT 3: Clipped view (bottom-left) ========== + plotter.subplot( 1, 0 ) + plotter.add_text( "Clipped View - Contributing Cells", font_size=14, position='upper_edge' ) + + # Determine clip position (middle of fault) + bounds = self.faultSurface.bounds + clip_normal = [ 0, 0, -1 ] # Clip along Z axis + clip_origin = [ 0, 0, ( bounds[ 4 ] + bounds[ 5 ] ) / 2 ] + + # Clip and show contributing cells + if self.contributingCells.n_cells > 0: + plotter.add_mesh_clip_plane( self.contributingCells, + normal=clip_normal, + origin=clip_origin, + color='blue', + opacity=1, + show_edges=True, + label='Contributing (clipped)' ) + + # Fault surface + plotter.add_mesh( self.faultSurface, color='red', opacity=1.0, show_edges=True, label='Fault' ) + + plotter.add_legend( loc='upper left' ) + plotter.add_axes() + plotter.set_scale( zscale=zscale ) + + # ========== PLOT 4: Slice view (bottom-right) ========== + plotter.subplot( 1, 1 ) + + # Determine slice position (middle of fault in Z) + slice_position = ( bounds[ 4 ] + bounds[ 5 ] ) / 2 + plotter.add_text( f"Slice View at Z={slice_position:.1f}m", font_size=14, position='upper_edge' ) + + # Create slice of volume + sliceVol = self.volumeMesh.slice( normal='z', origin=[ 0, 0, slice_position ] ) + sliceFault = self.faultSurface.slice( normal='z', origin=[ 0, 0, slice_position ] ) + + # Show contributing vs non-contributing in slice + if 'contributionSide' in sliceVol.cell_data: + # Non-contributing cells (gray) + nonContribMask = sliceVol.cell_data[ 'contributionSide' ] == 0 + if np.sum( nonContribMask ) > 0: + nonContrib = sliceVol.extract_cells( nonContribMask ) + plotter.add_mesh( nonContrib, + color='lightgray', + opacity=0.15, + show_edges=True, + line_width=1, + label='Non-contributing' ) + + # Plus side (blue) + plusMask = (sliceVol.cell_data['contributionSide'] == 1) | \ + (sliceVol.cell_data['contributionSide'] == 3) + if np.sum( plusMask ) > 0: + plusCells = sliceVol.extract_cells( plusMask ) + plotter.add_mesh( plusCells, + color='dodgerblue', + opacity=0.7, + show_edges=True, + line_width=2, + label='Plus side' ) + + # Minus side (orange) + minusMask = (sliceVol.cell_data['contributionSide'] == 2) | \ + (sliceVol.cell_data['contributionSide'] == 3) + if np.sum( minusMask ) > 0: + minusCells = sliceVol.extract_cells( minusMask ) + plotter.add_mesh( minusCells, + color='darkorange', + opacity=0.7, + show_edges=True, + line_width=2, + label='Minus side' ) + + # Fault slice (thick red line) + if sliceFault.n_cells > 0: + plotter.add_mesh( sliceFault, color='red', line_width=6, label='Fault', render_lines_as_tubes=True ) + + plotter.add_legend( loc='upper right' ) + plotter.add_axes() + plotter.set_scale( zscale=zscale ) + plotter.view_xy() + + # Link all views for synchronized rotation + plotter.link_views() + + # Show or save + if showPlots: + plotter.show() + else: + # Save screenshot + outputDir = self.outputDir + screenshot_path = outputDir / "contribution_visualization.png" + plotter.screenshot( str( screenshot_path ) ) + print( f" 💾 Visualization saved: {screenshot_path}" ) + plotter.close() + + # ------------------------------------------------------------------- + # NORMALS + # ------------------------------------------------------------------- + def _extractAndComputeNormals( self: Self, + showPlot: bool = False, + scaleFactor: float = 50.0, + zScale: float = 1.0 ) -> tuple[ pv.DataSet, list[ pv.DataSet ] ]: + """Extract fault surfaces and compute oriented normals/tangents.""" + surfaces = [] + mb = vtkMultiBlockDataSet() + mb.SetNumberOfBlocks( len( self.faultValues ) ) + + for i, faultId in enumerate( self.faultValues ): + # Extract fault cells + faultMask = np.where( getArrayInObject( self.mesh, self.faultAttribute, piece=Piece.CELLS ) == faultId )[0] + faultCells = extractCellSelection( self.mesh, ids=faultMask ) + + if faultCells.GetNumberOfCells() == 0: + print( f"⚠️ No cells for fault {faultId}" ) + continue + + # Extract surface + surf = extractSurface( faultCells ) + if surf.GetNumberOfCells() == 0: + continue + + # Compute normals + surf = computeNormals( surf, pointNormals=True ) + + # Orient normals consistently within the fault + surf = self._orientNormals( surf ) + + mb.SetBlock( i, surf) + surfaces.append( surf ) + + merged = mergeBlocks( mb, keepPartialAttributes=True) + # merged = pv.MultiBlock( surfaces ).combine() + print( f"✅ Normals computed for {merged.GetNumberOfCells()} fault cells" ) + + # if showPlot: + # self.plotGeometry( merged, scaleFactor, zScale ) + + return merged, surfaces + + # ------------------------------------------------------------------- + def _orientNormals( self: Self, surf: pv.PolyData, rotateNormals: bool = False ) -> pv.DataSet: + """Ensure normals point in consistent direction within the fault.""" + normals = vtk_to_numpy( surf.GetCellData().GetNormals() ) + meanNormal = np.mean( normals, axis=0 ) + meanNormal /= np.linalg.norm( meanNormal ) + + nCells = len( normals ) + tangents1 = np.zeros( ( nCells, 3 ) ) + tangents2 = np.zeros( ( nCells, 3 ) ) + + for i, normal in enumerate( normals ): + + # Flip if pointing opposite to mean + if np.dot( normal, meanNormal ) < 0: + normals[ i ] = -normal + + if rotateNormals: + normals[ i ] = -normal + + # Compute orthogonal tangents + normal = normals[ i ] + if abs( normal[ 0 ] ) > 1e-6 or abs( normal[ 1 ] ) > 1e-6: + t1 = np.array( [ -normal[ 1 ], normal[ 0 ], 0 ] ) + else: + t1 = np.array( [ 0, -normal[ 2 ], normal[ 1 ] ] ) + + t1 /= np.linalg.norm( t1 ) + t2 = np.cross( normal, t1 ) + t2 /= np.linalg.norm( t2 ) + + tangents1[ i ] = t1 + tangents2[ i ] = t2 + + surf.GetCellData().SetNormals( numpy_to_vtk( normals.ravel() ) ) + + createAttribute( surf, tangents1, "Tangents1" ) + createAttribute( surf, tangents2, "Tangents2" ) + surf.GetCellData().SetActiveTangents( "Tangents1" ) + + dipAngles, strikeAngles = self.computeDipStrikeFromCellBase( normals, tangents1, tangents2 ) + + createAttribute( surf, dipAngles, "dipAngle" ) + createAttribute( surf, strikeAngles, "strikeAngle" ) + + return surf + + # ------------------------------------------------------------------- + def computeDipStrikeFromCellBase( + self: Self, normals: npt.NDArray[ np.float64 ], tangent1: npt.NDArray[ np.float64 ], + tangent2: npt.NDArray[ np.float64 ] + ) -> tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ]: # TODO translate docstring + """Calcule les angles dip et strike à partir des vecteurs normaux et tangents des cellules. + + Hypothèses : + - Système de coordonnées : X=Est, Y=Nord, Z=Haut. + - Vecteurs donnés par cellule (shape: (nCells, 3)). + - Les vecteurs d'entrée sont supposés orthonormés (n = t1 x t2). + + Retourne : + dipDeg, strikeDeg (two arrays of shape (nCells,)) + """ + # 1. Identifier le vecteur strike (le plus horizontal) + t1Horizontal = tangent1 - ( tangent1[ :, 2 ][ :, np.newaxis ] * np.array( [ 0, 0, 1 ] ) ) + t2Horizontal = tangent2 - ( tangent2[ :, 2 ][ :, np.newaxis ] * np.array( [ 0, 0, 1 ] ) ) + normT1Horizontal = np.linalg.norm( t1Horizontal, axis=1 ) + normT2Horizontal = np.linalg.norm( t2Horizontal, axis=1 ) + + useT1 = normT1Horizontal > normT2Horizontal + strikeVector = np.zeros_like( tangent1 ) + strikeVector[ useT1 ] = t1Horizontal[ useT1 ] + strikeVector[ ~useT1 ] = t2Horizontal[ ~useT1 ] + + # Normaliser + strikeNorm = np.linalg.norm( strikeVector, axis=1 ) + # Éviter la division par zéro (si la faille est parfaitement verticale, le strike est bien défini par l'autre vecteur) + strikeNorm[ strikeNorm == 0 ] = 1.0 + strikeVector = strikeVector / strikeNorm[ :, np.newaxis ] + + # 2. Calculer le strike (azimut depuis le Nord, sens horaire) + strikeRad = np.arctan2( strikeVector[ :, 0 ], strikeVector[ :, 1 ] ) # atan2(E, N) + strikeDeg = np.degrees( strikeRad ) + strikeDeg = np.where( strikeDeg < 0, strikeDeg + 360, strikeDeg ) + + # 3. Calculer le dip + normHorizontal = np.linalg.norm( normals[ :, :2 ], axis=1 ) + dipRad = np.arcsin( np.clip( normHorizontal, 0, 1 ) ) # clip pour éviter les erreurs d'arrondi + dipDeg = np.degrees( dipRad ) + + return dipDeg, strikeDeg + + # ------------------------------------------------------------------- + def plotGeometry( self: Self, surface: pv.DataSet, scaleFactor: float, zScale: float ) -> None: + """Visualize fault geometry with normals.""" + plotter = pv.Plotter() + plotter.add_mesh( self.mesh, color='lightgray', opacity=0.1, label='Volume' ) + plotter.add_mesh( surface, color='darkgray', opacity=0.7, show_edges=True, label='Fault' ) + + centers = surface.cell_centers() + for name, color in [ ( 'Normals', 'red' ), ( 'tangent1', 'green' ), ( 'tangent2', 'blue' ) ]: + arrows = centers.glyph( orient=name, scale=zScale, factor=scaleFactor ) + plotter.add_mesh( arrows, color=color, label=name ) + + plotter.add_legend() + plotter.add_axes() + plotter.set_scale( zscale=zScale ) + plotter.show() + + # ------------------------------------------------------------------- + def diagnoseNormals( self: Self, scaleFactor: float = 50.0, zScale: float = 1.0 ) -> pv.DataSet: + """Diagnostic visualization to check normal quality. + + Shows orthogonality and orientation issues. + """ + surface = self.faultSurface + + print( "\n🔍 DIAGNOSTIC DES NORMALES" ) + print( "=" * 60 ) + + normals = surface.GetCellData().GetNormals() + tangent1 = surface.GetCellData().GetTangents() + tangent2 = surface.GetCellData().GetArray( "Tangents2" ) + + nCells = len( normals ) + + # Check orthogonality + dotNormT1 = np.array( [ np.dot( normals[ i ], tangent1[ i ] ) for i in range( nCells ) ] ) + dotNormT2 = np.array( [ np.dot( normals[ i ], tangent2[ i ] ) for i in range( nCells ) ] ) + dotT1T2 = np.array( [ np.dot( tangent1[ i ], tangent2[ i ] ) for i in range( nCells ) ] ) + + print( "Orthogonalité (doit être proche de 0):" ) + print( f" Normal · Tangent1 : max={np.max(np.abs(dotNormT1)):.2e}, mean={np.mean(np.abs(dotNormT1)):.2e}" ) + print( f" Normal · Tangent2 : max={np.max(np.abs(dotNormT2)):.2e}, mean={np.mean(np.abs(dotNormT2)):.2e}" ) + print( f" Tangent1 · Tangent2: max={np.max(np.abs(dotT1T2)):.2e}, mean={np.mean(np.abs(dotT1T2)):.2e}" ) + + # Check unit vectors + normN = np.linalg.norm( normals, axis=1 ) + normT1 = np.linalg.norm( tangent1, axis=1 ) + normT2 = np.linalg.norm( tangent2, axis=1 ) + + # Check unit vectors + # normN = np.array( [ np.linalg.norm( normals[ i ] ) for i in range( nCells ) ] ) # TODO + # normT1 = np.array( [ np.linalg.norm( tangent1[ i ] ) for i in range( nCells ) ] ) + # normT2 = np.array( [ np.linalg.norm( tangent2[ i ] ) for i in range( nCells ) ] ) + + print( "\nNormes (doit être proche de 1):" ) + print( f" Normals : min={np.min(normN):.6f}, max={np.max(normN):.6f}" ) + print( f" Tangent1 : min={np.min(normT1):.6f}, max={np.max(normT1):.6f}" ) + print( f" Tangent2 : min={np.min(normT2):.6f}, max={np.max(normT2):.6f}" ) + + # Check orientation consistency + meanNormal = np.mean( normals, axis=0 ) + meanNormal = meanNormal / np.linalg.norm( meanNormal ) + + dotsWithMean = np.array( [ np.dot( normals[ i ], meanNormal ) for i in range( nCells ) ] ) + nReversed = np.sum( dotsWithMean < 0 ) + + print( "\nCohérence d'orientation:" ) + print( f" Normale moyenne: [{meanNormal[0]:.3f}, {meanNormal[1]:.3f}, {meanNormal[2]:.3f}]" ) + print( f" Normales inversées: {nReversed}/{nCells} ({nReversed/nCells*100:.1f}%)" ) + + # Visual check + if nReversed > nCells * 0.1: + print( " ⚠️ Plus de 10% des normales pointent dans la direction opposée!" ) + else: + print( " ✅ Orientation cohérente" ) + + # Check for problematic cells + badOrtho = ( np.abs( dotNormT1 ) > 1e-3 ) | ( np.abs( dotNormT2 ) > 1e-3 ) | ( np.abs( dotT1T2 ) > 1e-3 ) + nBad = np.sum( badOrtho ) + + if nBad > 0: + print( f"\n⚠️ {nBad} cellules avec orthogonalité douteuse (|dot| > 1e-3)" ) + surface.cell_data[ 'orthogonality_error' ] = np.maximum.reduce( + [ np.abs( dotNormT1 ), np.abs( dotNormT2 ), + np.abs( dotT1T2 ) ] ) + else: + print( "\n✅ Toutes les cellules ont une bonne orthogonalité" ) + + print( "=" * 60 ) + + # Visualization + # plotter = pv.Plotter( shape=( 1, 2 ) ) + + # # Plot 1: Surface with normals + # plotter.subplot( 0, 0 ) + # plotter.add_mesh( surface, color='lightgray', show_edges=True, opacity=0.8 ) + + # centers = surface.cell_centers() + # arrowsNorm = centers.glyph( orient='Normals', scale=False, factor=scaleFactor ) + # plotter.add_mesh( arrowsNorm, color='red', label='Normals' ) + + # plotter.add_legend() + # plotter.add_axes() + # plotter.add_text( "Normales (Rouge)", position='upper_edge' ) + # plotter.set_scale( zscale=zScale ) + + # # Plot 2: All vectors + # plotter.subplot( 0, 1 ) + # plotter.add_mesh( surface, color='lightgray', show_edges=True, opacity=0.5 ) + + # arrowsNorm = centers.glyph( orient='Normals', scale=False, factor=scaleFactor ) + # arrowsT1 = centers.glyph( orient='tangent1', scale=False, factor=scaleFactor ) + # arrowsT2 = centers.glyph( orient='tangent2', scale=False, factor=scaleFactor ) + + # plotter.add_mesh( arrowsNorm, color='red', label='Normal' ) + # plotter.add_mesh( arrowsT1, color='green', label='Tangent1' ) + # plotter.add_mesh( arrowsT2, color='blue', label='Tangent2' ) + + # plotter.add_legend() + # plotter.add_axes() + # plotter.add_text( "Système complet (R,G,B)", position='upper_edge' ) + # plotter.set_scale( zscale=zScale ) + + # plotter.link_views() + # plotter.show() + + return surface diff --git a/geos-processing/src/geos/processing/post_processing/FaultStabilityAnalysis.py b/geos-processing/src/geos/processing/post_processing/FaultStabilityAnalysis.py new file mode 100755 index 000000000..65be56e68 --- /dev/null +++ b/geos-processing/src/geos/processing/post_processing/FaultStabilityAnalysis.py @@ -0,0 +1,199 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2026 TotalEnergies. +# SPDX-FileContributor: Nicolas Pillardou, Paloma Martinez +from pathlib import Path +import numpy as np +from typing_extensions import Self + +from geos.mesh.utils.multiblockModifiers import mergeBlocks +from geos.mesh.io.vtkIO import PVDReader +from geos.mesh.utils.arrayHelpers import (isAttributeInObject, getAttributeSet) + +from geos.processing.post_processing.FaultGeometry import FaultGeometry +from geos.processing.tools.FaultVisualizer import Visualizer +from geos.processing.post_processing.SensitivityAnalyzer import SensitivityAnalyzer +from geos.processing.post_processing.StressProjector import StressProjector +from geos.processing.post_processing.MohrCoulomb import MohrCoulomb + +from geos.utils.pieceEnum import Piece +# ============================================================================ +# TIME SERIES PROCESSING +# ============================================================================ +class TimeSeriesProcessor: + """Process multiple time steps from PVD file.""" + + # ------------------------------------------------------------------- + def __init__( self: Self, outputDir: str = ".", showPlots: bool = True, savePlots: bool = True ) -> None: + """Init.""" + self.outputDir = Path( outputDir ) + self.outputDir.mkdir( exist_ok=True ) + + self.showPlots: bool = showPlots + self.savePlots: bool = savePlots + + # ------------------------------------------------------------------- + def process( self: Self, + path: Path, + faultGeometry: FaultGeometry, + pvdFile: str, + timeIndexes: list[ int ] = [], + weightingScheme: str = "arithmetic", + cohesion: float = 0, + frictionAngle: float = 10, + runSensitivity: bool = True, + profileStartPoints: list[tuple[ float, ...]] = [], + computePrincipalStress: bool = False, + showDepthProfiles: bool = True, + stressName: str = "averageStress", + biotCoefficient: str = "rockPorosity_biotCoefficient", + profileSearchRadius=None, + minDepthProfiles=None, + maxDepthProfiles=None, + sensitivityFrictionAngles=None, + sensitivityCohesions=None ): + """Process all time steps using pre-computed fault geometry. + + Parameters: + path: base path for input files + faultGeometry: FaultGeometry object with initialized topology + pvdFile: PVD file name + """ + reader = PVDReader( path / pvdFile ) + timeValues = reader.getAllTimestepsValues() + + if timeIndexes: + timeValues = timeValues[ timeIndexes ] + + outputFiles = [] + dataInitial = None + + # Get pre-computed data from faultGeometry + surface = faultGeometry.faultSurface + adjacencyMapping = faultGeometry.adjacencyMapping + geometricProperties = faultGeometry.getGeometricProperties() + + # Initialize projector with pre-computed topology + projector = StressProjector( adjacencyMapping, geometricProperties, self.outputDir ) + + print( '\n' ) + print( "=" * 60 ) + print( "TIME SERIES PROCESSING" ) + print( "=" * 60 ) + + for i, time in enumerate( timeValues ): + print( f"\n→ Step {i+1}/{len(timeValues)}: {time/(365.25*24*3600):.2f} years" ) + + # Read time step + idx = timeIndexes[ i ] if timeIndexes else i + dataset = reader.getDataSetAtTimeIndex( idx ) + + # Merge blocks + volumeData = mergeBlocks( dataset, keepPartialAttributes=True ) + + if dataInitial is None: + dataInitial = volumeData + + + # ----------------------------------- + # Projection using pre-computed topology + # ----------------------------------- + # Projection + surfaceResult, volumeMarked, contributingCells = projector.projectStressToFault( + volumeData, + dataInitial, + surface, + time=timeValues[ i ], # Simulation time + timestep=i, # Timestep index + stressName=stressName, + biotName=biotCoefficient, + weightingScheme=weightingScheme ) + + # ----------------------------------- + # Mohr-Coulomb analysis + # ----------------------------------- + cohesion = cohesion + frictionAngle = frictionAngle + surfaceResult = MohrCoulomb.analyze( surfaceResult, cohesion, frictionAngle ) #, time ) + + # ----------------------------------- + # Visualize + # ----------------------------------- + self._plotResults( surfaceResult, contributingCells, time, self.outputDir, profileStartPoints, computePrincipalStress, showDepthProfiles, + profileSearchRadius, minDepthProfiles, maxDepthProfiles ) + + # ----------------------------------- + # Sensitivity analysis + # ----------------------------------- + if runSensitivity: + analyzer = SensitivityAnalyzer( self.outputDir, self.showPlots ) + if sensitivityFrictionAngles is None or sensitivityCohesions is None: + raise ValueError( "sensitivity friction angles and cohesions required if runSensitivity is set to True" ) + analyzer.runAnalysis( surfaceResult, time, sensitivityFrictionAngles, sensitivityCohesions, profileStartPoints, profileSearchRadius ) + + # Save + filename = f'fault_analysis_{i:04d}.vtu' + # surfaceResult.save( self.outputDir / filename ) + outputFiles.append( ( time, filename ) ) + print( f" 💾 Saved: {filename}" ) + + # Create master PVD + self._createPVD( outputFiles ) + + return surfaceResult + # ------------------------------------------------------------------- + def _plotResults( self, surface, contributingCells, time: list[ int ], + path: str, profileStartPoints: list[tuple[float, ...]], computePrincipalStress: bool = False, showDepthProfiles:bool = True, + profileSearchRadius: float|None=None, minDepthProfiles: float | None = None, + maxDepthProfiles: float | None = None, ) -> None: # TODO check type surface + Visualizer.plotMohrCoulombDiagram( surface, + time, + path, + show=self.showPlots, + save=self.savePlots, ) + + + # Profils verticaux automatiques + if showDepthProfiles: + Visualizer( profileSearchRadius).plotDepthProfiles( surface=surface, + time=time, + path=path, + show=self.showPlots, + save=self.savePlots, + profileStartPoints=profileStartPoints) + + + visualizer = Visualizer( profileSearchRadius, + minDepthProfiles, + maxDepthProfiles, + showPlots = self.showPlots, savePlots = self.savePlots ) + + if computePrincipalStress: + # Plot principal stress from volume cells + visualizer.plotVolumeStressProfiles( volumeMesh=contributingCells, + faultSurface=surface, + time=time, + path=path, + profileStartPoints=profileStartPoints ) + + # Visualize comparison analytical/numerical + visualizer.plotAnalyticalVsNumericalComparison( volumeMesh=contributingCells, + faultSurface=surface, + time=time, + path=path, + show=self.showPlots, + save=self.savePlots, + profileStartPoints=profileStartPoints ) + + # ------------------------------------------------------------------- + def _createPVD( self, outputFiles: list[ tuple[ int, str ] ] ) -> None: + """Create PVD collection file.""" + pvdPath = self.outputDir / 'fault_analysis.pvd' + with open( pvdPath, 'w' ) as f: + f.write( '\n' ) + f.write( ' \n' ) + for t, fname in outputFiles: + f.write( f' \n' ) + f.write( ' \n' ) + f.write( '\n' ) + print( f"\n✅ PVD created: {pvdPath}" ) + diff --git a/geos-processing/src/geos/processing/post_processing/MohrCoulomb.py b/geos-processing/src/geos/processing/post_processing/MohrCoulomb.py new file mode 100644 index 000000000..6a6dbcb18 --- /dev/null +++ b/geos-processing/src/geos/processing/post_processing/MohrCoulomb.py @@ -0,0 +1,111 @@ +import numpy as np +from vtkmodules.vtkCommonDataModel import ( + vtkDataSet ) +from vtkmodules.util.numpy_support import numpy_to_vtk +from geos.mesh.utils.arrayHelpers import ( getArrayInObject, isAttributeInObject ) +from geos.mesh.utils.arrayModifiers import ( createAttribute, updateAttribute ) +from geos.utils.pieceEnum import Piece +# ============================================================================ +# MOHR COULOMB +# ============================================================================ +class MohrCoulomb: + """Mohr-Coulomb failure criterion analysis.""" + + @staticmethod + def analyze( surface: vtkDataSet, cohesion: float, frictionAngleDeg: float, verbose: bool = True ) -> vtkDataSet: + """Perform Mohr-Coulomb stability analysis. + + Parameters: + surface: fault surface with stress data + cohesion: cohesion in bar + frictionAngleDeg: friction angle in degrees + verbose: print statistics + """ + mu = np.tan( np.radians( frictionAngleDeg ) ) + + # Extract stress components + sigmaN = getArrayInObject( surface, "sigmaNEffective", Piece.CELLS ) + tau = getArrayInObject( surface, "tauEffective", Piece.CELLS ) + deltaSigmaN = getArrayInObject( surface, 'deltaSigmaNEffective', Piece.CELLS ) + deltaTau = getArrayInObject( surface, 'deltaTauEffective', Piece.CELLS ) + + # Mohr-Coulomb failure envelope + tauCritical = cohesion - sigmaN * mu + + # Coulomb Failure Stress + CFS = tau - mu * sigmaN + # deltaCFS = deltaTau - mu * deltaSigmaN + + # Shear Capacity Utilization: SCU = τ / τ_crit + SCU = np.divide( tau, tauCritical, out=np.zeros_like( tau ), where=tauCritical != 0 ) + + # if "SCUInitial" not in surface.cell_data: + if not isAttributeInObject( surface, "SCUInitial", Piece.CELLS ): + # First timestep: store as initial reference + SCUInitial = SCU.copy() + CFSInitial = CFS.copy() + deltaSCU = np.zeros_like( SCU ) + deltaCFS = np.zeros_like( CFS ) + + createAttribute( surface, SCUInitial, "SCUInitial" ) + createAttribute( surface, CFSInitial, "CFSInitial" ) + + isInitial = True + else: + # Subsequent timesteps: calculate change from initial + SCUInitial = getArrayInObject( surface, "SCUInitial", Piece.CELLS ) + CFSInitial = getArrayInObject( surface, "CFSInitial", Piece.CELLS ) + deltaSCU = SCU - SCUInitial + deltaCFS = CFS - CFSInitial + isInitial = False + + # Stability classification + stability = np.zeros_like( tau, dtype=int ) + stability[ SCU >= 0.8 ] = 1 # Critical + stability[ SCU >= 1.0 ] = 2 # Unstable + + # Failure probability (sigmoid) + k = 10.0 + failureProba = 1.0 / ( 1.0 + np.exp( -k * ( SCU - 1.0 ) ) ) + + # Safety margin + safety = tauCritical - tau + + # Store results + attributes = { "mohrCohesion": np.full( surface.GetNumberOfCells(), cohesion ), + "mohrFrictionAngle": np.full( surface.GetNumberOfCells(), frictionAngleDeg ), + "mohrFrictionCoefficient": np.full( surface.GetNumberOfCells(), mu ), + "mohr_critical_shear_stress": tauCritical, + "SCU": SCU, + "deltaSCU": deltaSCU, + "CFS": CFS, + "deltaCFS": deltaCFS, + "safetyMargin": safety, + "stabilityState": stability, + "failureProbability": failureProba } + + cdata = surface.GetCellData() + for attributeName, value in attributes.items(): + updateAttribute( surface, value, attributeName, Piece.CELLS ) + + + if verbose: + nStable = np.sum( stability == 0 ) + nCritical = np.sum( stability == 1 ) + nUnstable = np.sum( stability == 2 ) + + # Additional info on deltaSCU + if not isInitial: + meanDelta = np.mean( np.abs( deltaSCU ) ) + maxIncrease = np.max( deltaSCU ) + maxDecrease = np.min( deltaSCU ) + print( f" ✅ Mohr-Coulomb: {nUnstable} unstable, {nCritical} critical, " + f"{nStable} stable cells" ) + print( f" ΔSCU: mean={meanDelta:.3f}, maxIncrease={maxIncrease:.3f}, " + f"maxDecrease={maxDecrease:.3f}" ) + else: + print( f" ✅ Mohr-Coulomb (initial): {nUnstable} unstable, {nCritical} critical, " + f"{nStable} stable cells" ) + + return surface + diff --git a/geos-processing/src/geos/processing/post_processing/ProfileExtractor.py b/geos-processing/src/geos/processing/post_processing/ProfileExtractor.py new file mode 100644 index 000000000..9ed491ddc --- /dev/null +++ b/geos-processing/src/geos/processing/post_processing/ProfileExtractor.py @@ -0,0 +1,730 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2026 TotalEnergies. +# SPDX-FileContributor: Nicolas Pillardou, Paloma Martinez +import numpy as np +import pyvista as pv +import numpy.typing as npt + +from vtkmodules.vtkCommonDataModel import vtkCellData, vtkDataSet +from vtkmodules.util.numpy_support import vtk_to_numpy + +# ============================================================================ +# PROFILE EXTRACTOR +# ============================================================================ +class ProfileExtractor: + """Utility class for extracting profiles along fault surfaces.""" + + # ------------------------------------------------------------------- + @staticmethod + def extractAdaptiveProfile( + centers: npt.NDArray[ np.float64 ], + values: npt.NDArray[ np.float64 ], + xStart: npt.NDArray[ np.float64 ], + yStart: npt.NDArray[ np.float64 ], + zStart: npt.NDArray[ np.float64 ] | None = None, + searchRadius: float | None = None, + stepSize: float = 20.0, + maxSteps: float = 500, + verbose: bool = True, + cellData: vtkCellData = None + ) -> tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], + npt.NDArray[ np.float64 ] ]: + """Extraction de profil vertical par COUCHES DE PROFONDEUR avec détection automatique de faille. + + Stratégie: + 1. Trouver le point de départ le plus proche + 2. Identifier automatiquement la faille via cellData (attribute, FaultMask, etc.) + 3. FILTRER pour ne garder QUE les cellules de cette faille + 4. Diviser en tranches Z + 5. Pour chaque tranche, prendre la cellule la plus proche en XY + + Parameters + ---------- + centers : ndarray + Cell centers (nCells, 3) + values : ndarray + Values at cells (nCells,) + xStart, yStart : float + Starting XY position + zStart : float, optional + Starting Z position (if None, uses highest point near XY) + searchRadius : float, optional + Not used (kept for compatibility) + cellData : dict, optional + Dictionary with cell data fields (e.g., {'attribute': array, 'FaultMask': array}) + Used to automatically detect and filter by fault ID + verbose : bool + Print detailed information + + Returns: + ------- + depths, profileValues, pathX, pathY : ndarrays + Extracted profile data + """ + # Convert to np arrays + centers = np.asarray( centers ) + values = np.asarray( values ) + + if len( centers ) == 0: + if verbose: + print( " ⚠️ No cells provided" ) + return np.array( [] ), np.array( [] ), np.array( [] ), np.array( [] ) + + # =================================================================== + # ÉTAPE 1: TROUVER LE POINT DE DÉPART + # =================================================================== + + if zStart is None: + # Chercher en 2D (XY), prendre le plus haut + if verbose: + print( f" Searching near ({xStart:.1f}, {yStart:.1f})" ) + + dXY = np.sqrt( ( centers[ :, 0 ] - xStart )**2 + ( centers[ :, 1 ] - yStart )**2 ) + closestIndices = np.argsort( dXY )[ :20 ] + + if len( closestIndices ) == 0: + print( " ⚠️ No cells found near start point" ) + return np.array( [] ), np.array( [] ), np.array( [] ), np.array( [] ) + + # Prendre le plus haut (plus grand Z) + closestDepths = centers[ closestIndices, 2 ] + startIdx = closestIndices[ np.argmax( closestDepths ) ] + else: + # Chercher en 3D + if verbose: + print( f" Searching near ({xStart:.1f}, {yStart:.1f}, {zStart:.1f})" ) + + d3D = np.sqrt( ( centers[ :, 0 ] - xStart )**2 + ( centers[ :, 1 ] - yStart )**2 + + ( centers[ :, 2 ] - zStart )**2 ) + startIdx = np.argmin( d3D ) + + startPoint = centers[ startIdx ] + + if verbose: + print( f" Starting point: ({startPoint[0]:.1f}, {startPoint[1]:.1f}, {startPoint[2]:.1f})" ) + print( f" Starting cell index: {startIdx}" ) + + # =================================================================== + # ÉTAPE 2: DÉTECTER AUTOMATIQUEMENT L'ID DE LA FAILLE + # =================================================================== + + faultIds = None + targetFaultId = None + + if cellData is not None: + # Chercher dans l'ordre de priorité + faultFieldNames = [ 'attribute', 'FaultMask', 'faultId', 'region' ] + + for fieldName in faultFieldNames: + if cellData.HasArray( fieldName ): + faultIds = vtk_to_numpy( cellData[ fieldName ] ) + + + if len( faultIds ) != len( centers ): + if verbose: + print( f" ⚠️ Field '{fieldName}' length mismatch, skipping" ) + continue + + # Récupérer l'ID au point de départ + targetFaultId = faultIds[ startIdx ] + + if verbose: + uniqueIds = np.unique( faultIds ) + print( f" Found fault field: '{fieldName}'" ) + print( f" Available fault IDs: {uniqueIds}" ) + print( f" Target fault ID at start point: {targetFaultId}" ) + + break + + # =================================================================== + # ÉTAPE 3: FILTRER PAR FAILLE SI DÉTECTÉE + # =================================================================== + + if targetFaultId is not None: + # FILTRER: garder SEULEMENT cette faille + maskSameFault = ( faultIds == targetFaultId ) + nTotal = len( centers ) + nOnFault = np.sum( maskSameFault ) + + if verbose: + print( + f" Filtering to fault ID={targetFaultId}: {nOnFault}/{nTotal} cells ({nOnFault/nTotal*100:.1f}%)" + ) + + if nOnFault == 0: + print( " ⚠️ No cells found on target fault" ) + return np.array( [] ), np.array( [] ), np.array( [] ), np.array( [] ) + + # REMPLACER centers et values par le subset filtré + centers = centers[ maskSameFault ].copy() + values = values[ maskSameFault ].copy() + + # Trouver le nouvel index de départ dans le subset + dToStart = np.sqrt( np.sum( ( centers - startPoint )**2, axis=1 ) ) + startIdx = np.argmin( dToStart ) + + if verbose: + print( f" ✅ Profile will stay on fault ID={targetFaultId}" ) + else: + if verbose: + print( " ⚠️ No fault identification field found" ) + if cellData is not None: + print( f" Available fields: {list(cellData.keys())}" ) + else: + print( " cellData not provided" ) + print( " Profile may jump between faults!" ) + + # À partir d'ici, centers/values ne contiennent QUE la faille cible + + # =================================================================== + # ÉTAPE 4: POSITION DE RÉFÉRENCE + # =================================================================== + + refX = centers[ startIdx, 0 ] + refY = centers[ startIdx, 1 ] + + if verbose: + print( f" Reference XY: ({refX:.1f}, {refY:.1f})" ) + + # =================================================================== + # ÉTAPE 5: GÉOMÉTRIE DE LA FAILLE + # =================================================================== + + xRange = np.max( centers[ :, 0 ] ) - np.min( centers[ :, 0 ] ) + yRange = np.max( centers[ :, 1 ] ) - np.min( centers[ :, 1 ] ) + zRange = np.max( centers[ :, 2 ] ) - np.min( centers[ :, 2 ] ) + + if zRange <= 0: + print( f" ⚠️ Invalid zRange: {zRange}" ) + return np.array( [] ), np.array( [] ), np.array( [] ), np.array( [] ) + + lateralExtent = max( xRange, yRange ) + xyTolerance = max( lateralExtent * 0.3, 100.0 ) + + if verbose: + print( f" Fault extent: X={xRange:.1f}m, Y={yRange:.1f}m, Z={zRange:.1f}m" ) + print( f" XY tolerance: {xyTolerance:.1f}m" ) + + # =================================================================== + # ÉTAPE 6: CALCUL DES TRANCHES + # =================================================================== + + zCoordsSorted = np.sort( centers[ :, 2 ] ) + zDiffs = np.diff( zCoordsSorted ) + zDiffsPositive = zDiffs[ zDiffs > 1e-6 ] + + if len( zDiffsPositive ) == 0: + if verbose: + print( " ⚠️ All cells at same Z" ) + + dXY = np.sqrt( ( centers[ :, 0 ] - refX )**2 + ( centers[ :, 1 ] - refY )**2 ) + sortedIndices = np.argsort( dXY ) + + return ( centers[ sortedIndices, 2 ], values[ sortedIndices ], centers[ sortedIndices, + 0 ], centers[ sortedIndices, 1 ] ) + + medianZSpacing = np.median( zDiffsPositive ) + + # Vérifier que medianZSpacing est raisonnable + if medianZSpacing <= 0 or medianZSpacing > zRange: + medianZSpacing = zRange / 100 # Fallback + + # Taille de tranche = espacement médian + sliceThickness = medianZSpacing + + zMin = np.min( centers[ :, 2 ] ) + zMax = np.max( centers[ :, 2 ] ) + + nSlices = int( np.ceil( zRange / sliceThickness ) ) + nSlices = min( nSlices, 10000 ) # Limiter à 10k tranches max + + if nSlices <= 0: + print( f" ⚠️ Invalid nSlices: {nSlices}" ) + return np.array( [] ), np.array( [] ), np.array( [] ), np.array( [] ) + + if verbose: + print( f" Median Z spacing: {medianZSpacing:.1f}m" ) + print( f" Creating {nSlices} slices" ) + + try: + zSlices = np.linspace( zMax, zMin, nSlices + 1 ) + except ( MemoryError, ValueError ) as e: + print( f" ⚠️ Error creating slices: {e}" ) + return np.array( [] ), np.array( [] ), np.array( [] ), np.array( [] ) + + # =================================================================== + # ÉTAPE 7: EXTRACTION PAR TRANCHES + # =================================================================== + + profileIndices = [] + + for i in range( len( zSlices ) - 1 ): + zTop = zSlices[ i ] + zBottom = zSlices[ i + 1 ] + + # Cellules dans cette tranche + maskInSlice = ( centers[ :, 2 ] <= zTop ) & ( centers[ :, 2 ] >= zBottom ) + indicesInSlice = np.where( maskInSlice )[ 0 ] + + if len( indicesInSlice ) == 0: + continue + + # Distance XY à la référence + dXYInSlice = np.sqrt( ( centers[ indicesInSlice, 0 ] - refX )**2 + + ( centers[ indicesInSlice, 1 ] - refY )**2 ) + + # Ne garder que celles dans la tolérance XY + validMask = dXYInSlice < xyTolerance + + if not np.any( validMask ): + # Aucune dans la tolérance → prendre la plus proche + closestInSlice = indicesInSlice[ np.argmin( dXYInSlice ) ] + else: + # Prendre la plus proche parmi celles dans la tolérance + validIndices = indicesInSlice[ validMask ] + dXYValid = dXYInSlice[ validMask ] + closestInSlice = validIndices[ np.argmin( dXYValid ) ] + + profileIndices.append( closestInSlice ) + + # =================================================================== + # ÉTAPE 8: SUPPRIMER DOUBLONS ET TRIER + # =================================================================== + + # Supprimer doublons + seen = set() + uniqueIndices = [] + for idx in profileIndices: + if idx not in seen: + seen.add( idx ) + uniqueIndices.append( idx ) + + if len( uniqueIndices ) == 0: + if verbose: + print( " ⚠️ No points extracted" ) + return np.array( [] ), np.array( [] ), np.array( [] ), np.array( [] ) + + profileIndicesArr = np.array( uniqueIndices ) + + # Trier par profondeur décroissante (haut → bas) + sortOrder = np.argsort( -centers[ profileIndicesArr, 2 ] ) + profileIndicesArr = profileIndicesArr[ sortOrder ] + + # Extraire résultats + depths = centers[ profileIndicesArr, 2 ] + profileValues = values[ profileIndicesArr ] + pathX = centers[ profileIndicesArr, 0 ] + pathY = centers[ profileIndicesArr, 1 ] + + # =================================================================== + # STATISTIQUES + # =================================================================== + + if verbose: + depthCoverage = ( depths.max() - depths.min() ) / zRange * 100 if zRange > 0 else 0 + xyDisplacement = np.sqrt( ( pathX[ -1 ] - pathX[ 0 ] )**2 + ( pathY[ -1 ] - pathY[ 0 ] )**2 ) + + print( f" ✅ Extracted {len(profileIndices)} points" ) + print( f" Depth range: [{depths.max():.1f}, {depths.min():.1f}]m" ) + print( f" Coverage: {depthCoverage:.1f}% of fault depth" ) + print( f" XY displacement: {xyDisplacement:.1f}m" ) + + return ( depths, profileValues, pathX, pathY ) + + # ------------------------------------------------------------------- + @staticmethod + def extractVerticalProfileTopologyBased( + surfaceMesh: pv.DataSet, + fieldName: str, + xStart: float, + yStart: float, + zStart: float | None = None, + maxSteps: int = 500, + verbose: bool = True + ) -> tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], + npt.NDArray[ np.float64 ] ]: + """Extraction de profil vertical en utilisant la TOPOLOGIE du maillage de surface.""" + if fieldName not in surfaceMesh.cell_data: + print( f" ⚠️ Field '{fieldName}' not found in mesh" ) + return np.array( [] ), np.array( [] ), np.array( [] ), np.array( [] ) + + centers = surfaceMesh.cell_centers().points + values = surfaceMesh.cell_data[ fieldName ] + + # =================================================================== + # ÉTAPE 1: TROUVER LA CELLULE DE DÉPART + # =================================================================== + + if zStart is None: + if verbose: + print( f" Searching near ({xStart:.1f}, {yStart:.1f})" ) + + dXY = np.sqrt( ( centers[ :, 0 ] - xStart )**2 + ( centers[ :, 1 ] - yStart )**2 ) + closestIndices = np.argsort( dXY )[ :20 ] + + if len( closestIndices ) == 0: + print( " ⚠️ No cells found" ) + return np.array( [] ), np.array( [] ), np.array( [] ), np.array( [] ) + + closestDepths = centers[ closestIndices, 2 ] + startIdx = closestIndices[ np.argmax( closestDepths ) ] + else: + if verbose: + print( f" Searching near ({xStart:.1f}, {yStart:.1f}, {zStart:.1f})" ) + + d3D = np.sqrt( ( centers[ :, 0 ] - xStart )**2 + ( centers[ :, 1 ] - yStart )**2 + + ( centers[ :, 2 ] - zStart )**2 ) + startIdx = np.argmin( d3D ) + + startPoint = centers[ startIdx ] + + if verbose: + print( f" Starting cell: {startIdx}" ) + print( f" Starting point: ({startPoint[0]:.1f}, {startPoint[1]:.1f}, {startPoint[2]:.1f})" ) + + # =================================================================== + # ÉTAPE 2: IDENTIFIER LA FAILLE + # =================================================================== + + targetFaultId = None + faultIds = None + faultFieldNames = [ 'attribute', 'FaultMask', 'faultId', 'region' ] + + for fieldNameCheck in faultFieldNames: + if fieldNameCheck in surfaceMesh.cell_data: + faultIds = surfaceMesh.cell_data[ fieldNameCheck ] + targetFaultId = faultIds[ startIdx ] + + if verbose: + uniqueIds = np.unique( faultIds ) + print( f" Fault field: '{fieldNameCheck}'" ) + print( f" Target fault ID: {targetFaultId} (from {uniqueIds})" ) + + break + + if targetFaultId is None and verbose: + print( " ⚠️ No fault ID found - will use all cells" ) + + # =================================================================== + # ÉTAPE 3: CONSTRUIRE LA CONNECTIVITÉ (VOISINS TOPOLOGIQUES) + # =================================================================== + + if verbose: + print( " Building cell connectivity..." ) + + nCells = surfaceMesh.n_cells + connectivity: list[ list[ int ] ] = [ [] for _ in range( nCells ) ] + + # Construire un dictionnaire arête -> cellules + edgeToCells: dict[ tuple[ int, int ], list[ int ] ] = {} + + for cellId in range( nCells ): + cell = surfaceMesh.get_cell( cellId ) + nPoints = cell.n_points + + # Pour chaque arête de la cellule + for i in range( nPoints ): + p1 = cell.point_ids[ i ] + p2 = cell.point_ids[ ( i + 1 ) % nPoints ] + + # Arête normalisée (ordre canonique) + edge = tuple( sorted( [ p1, p2 ] ) ) + + if edge not in edgeToCells: + edgeToCells[ edge ] = [] + edgeToCells[ edge ].append( cellId ) + + # Pour chaque cellule, trouver ses voisins via arêtes partagées + for cellId in range( nCells ): + cell = surfaceMesh.get_cell( cellId ) + nPoints = cell.n_points + + neighborsSet = set() + + for i in range( nPoints ): + p1 = cell.point_ids[ i ] + p2 = cell.point_ids[ ( i + 1 ) % nPoints ] + edge = tuple( sorted( [ p1, p2 ] ) ) + + # Toutes les cellules partageant cette arête sont voisines + for neighborId in edgeToCells[ edge ]: + if neighborId != cellId: + neighborsSet.add( neighborId ) + + connectivity[ cellId ] = list( neighborsSet ) + + if verbose: + avgNeighbors = np.mean( [ len( c ) for c in connectivity ] ) + maxNeighbors = np.max( [ len( c ) for c in connectivity ] ) + print( f" Connectivity built: avg={avgNeighbors:.1f} neighbors/cell, max={maxNeighbors}" ) + + # =================================================================== + # ÉTAPE 4: ALGORITHME DE DESCENTE PAR VOISINAGE TOPOLOGIQUE + # =================================================================== + + profileIndices = [ startIdx ] + visited = { startIdx } + currentIdx = startIdx + + refXY = startPoint[ :2 ] # Position XY de référence + + if verbose: + print( f" Starting descent from Z={startPoint[2]:.1f}m..." ) + + stuckCount = 0 + maxStuck = 3 + + for step in range( maxSteps ): + currentZ = centers[ currentIdx, 2 ] + + # Obtenir les voisins topologiques + neighborIndices = connectivity[ currentIdx ] + + # Filtrer les voisins: + # 1. Non visités + # 2. Même faille (si détectée) + # 3. Plus bas en Z + candidates = [] + + for idx in neighborIndices: + if idx in visited: + continue + + # Vérifier la faille + if targetFaultId is not None and faultIds is not None and faultIds[ idx ] != targetFaultId: + continue + + # Vérifier qu'on descend + if centers[ idx, 2 ] >= currentZ: + continue + + candidates.append( idx ) + + if len( candidates ) == 0: + # Si bloqué, essayer de regarder les voisins des voisins + stuckCount += 1 + + if stuckCount >= maxStuck: + if verbose: + print( f" Reached bottom at Z={currentZ:.1f}m after {step+1} steps (no more neighbors)" ) + break + + # Essayer niveau 2 (voisins des voisins) + extendedCandidates = [] + for neighborIdx in neighborIndices: + if neighborIdx in visited: + continue + + for secondNeighborIdx in connectivity[ neighborIdx ]: + if secondNeighborIdx in visited: + continue + + if targetFaultId is not None and faultIds is not None and faultIds[ + secondNeighborIdx ] != targetFaultId: + continue + + if centers[ secondNeighborIdx, 2 ] < currentZ: + extendedCandidates.append( secondNeighborIdx ) + + if len( extendedCandidates ) == 0: + if verbose: + print( f" Reached bottom at Z={currentZ:.1f}m (extended search failed)" ) + break + + candidates = extendedCandidates + if verbose: + print( f" Used extended search at step {step+1}" ) + else: + stuckCount = 0 + + # Parmi les candidats, choisir celui le plus proche en XY de la référence + bestIdx = None + bestDistanceXY = float( 'inf' ) + + for idx in candidates: + pos = centers[ idx ] + dXY = np.sqrt( ( pos[ 0 ] - refXY[ 0 ] )**2 + ( pos[ 1 ] - refXY[ 1 ] )**2 ) + + if dXY < bestDistanceXY: + bestDistanceXY = dXY + bestIdx = idx + + if bestIdx is None: + if verbose: + print( f" No valid neighbor at Z={currentZ:.1f}m" ) + break + + # Ajouter au profil + profileIndices.append( bestIdx ) + visited.add( bestIdx ) + currentIdx = bestIdx + + # Debug + if verbose and ( step + 1 ) % 100 == 0: + print( + f" Step {step+1}: Z={centers[currentIdx, 2]:.1f}m, XY=({centers[currentIdx, 0]:.1f}, {centers[currentIdx, 1]:.1f})" + ) + + # =================================================================== + # ÉTAPE 5: EXTRAIRE LES RÉSULTATS + # =================================================================== + + if len( profileIndices ) == 0: + if verbose: + print( " ⚠️ No profile extracted" ) + return np.array( [] ), np.array( [] ), np.array( [] ), np.array( [] ) + + depths = centers[ np.array( profileIndices ), 2 ] + profileValues = values[ np.array( profileIndices ) ] + pathX = centers[ np.array( profileIndices ), 0 ] + pathY = centers[ np.array( profileIndices ), 1 ] + + # =================================================================== + # STATISTIQUES + # =================================================================== + + if verbose: + zRange = np.max( centers[ :, 2 ] ) - np.min( centers[ :, 2 ] ) + depthCoverage = ( depths.max() - depths.min() ) / zRange * 100 if zRange > 0 else 0 + xyDisplacement = np.sqrt( ( pathX[ -1 ] - pathX[ 0 ] )**2 + ( pathY[ -1 ] - pathY[ 0 ] )**2 ) + + print( f" ✅ {len(profileIndices)} points extracted" ) + print( f" Depth range: [{depths.max():.1f}, {depths.min():.1f}]m" ) + print( f" Coverage: {depthCoverage:.1f}% of fault depth" ) + print( f" XY displacement: {xyDisplacement:.1f}m" ) + + return ( depths, profileValues, pathX, pathY ) + + # ------------------------------------------------------------------- + @staticmethod + def plotProfilePath3D( surface: pv.DataSet, + pathX: npt.NDArray[ np.float64 ], + pathY: npt.NDArray[ np.float64 ], + pathZ: npt.NDArray[ np.float64 ], + profileValues: npt.NDArray[ np.float64 ] | None = None, + scalarName: str = 'SCU', + savePath: str | None = None, + show: bool = True ) -> None: + """Visualize the extracted profile path on the fault surface in 3D using PyVista. + + Parameters + ---------- + surface : pyvista.PolyData + Fault surface mesh + pathX, pathY, pathZ : array-like + Coordinates of the profile path + profileValues : array-like, optional + Values along the profile (for coloring the path) + scalarName : str + Name of the scalar to display on the surface + savePath : Path or str, optional + Path to save the screenshot + show : bool + Whether to display the plot interactively + """ + if len( pathX ) == 0: + print( " ⚠️ No path to plot (empty profile)" ) + return + + print( f" 📊 Creating 3D visualization of profile path ({len(pathX)} points)" ) + + # Create plotter + plotter = pv.Plotter( window_size=[ 1600, 1200 ] ) + + # Add fault surface with scalar field + if scalarName in surface.cell_data: + plotter.add_mesh( surface, + scalars=scalarName, + cmap='RdYlGn_r', + opacity=0.7, + show_edges=False, + lighting=True, + smooth_shading=True, + scalar_bar_args={ + 'title': scalarName, + 'title_font_size': 20, + 'label_font_size': 16, + 'n_labels': 5, + 'italic': False, + 'fmt': '%.2f', + 'font_family': 'arial', + } ) + else: + plotter.add_mesh( surface, color='lightgray', opacity=0.5, show_edges=True ) + + # Create path as a polyline + pathPoints = np.column_stack( [ pathX, pathY, pathZ ] ) + pathPolyline = pv.PolyData( pathPoints ) + + # Add connectivity for line + nPoints = len( pathPoints ) + lines = np.full( ( nPoints - 1, 3 ), 2, dtype=np.int_ ) + lines[ :, 1 ] = np.arange( nPoints - 1 ) + lines[ :, 2 ] = np.arange( 1, nPoints ) + pathPolyline.lines = lines.ravel() + + # Color the path by profile values or depth + if profileValues is not None: + pathPolyline[ 'profile_value' ] = profileValues + colorField = 'profile_value' + cmapPath = 'viridis' + else: + pathPolyline[ 'depth' ] = pathZ + colorField = 'depth' + cmapPath = 'turbo_r' + + # Add path as thick tube + pathTube = pathPolyline.tube( radius=10.0 ) # Adjust radius as needed + plotter.add_mesh( pathTube, + scalars=colorField, + cmap=cmapPath, + line_width=8, + render_lines_as_tubes=True, + lighting=True, + scalar_bar_args={ + 'title': 'Path ' + colorField, + 'title_font_size': 20, + 'label_font_size': 16, + 'position_x': 0.85, + 'position_y': 0.05, + } ) + + # Add start and end markers + startPoint = pv.Sphere( radius=30, center=pathPoints[ 0 ] ) + endPoint = pv.Sphere( radius=30, center=pathPoints[ -1 ] ) + + plotter.add_mesh( startPoint, color='lime', label='Start (Top)' ) + plotter.add_mesh( endPoint, color='red', label='End (Bottom)' ) + + # Add axes and labels + plotter.add_axes( xlabel='X [m]', ylabel='Y [m]', zlabel='Z [m]', line_width=3, labels_off=False ) + + # Add legend + plotter.add_legend( labels=[ ( 'Start (Top)', 'lime' ), ( 'End (Bottom)', 'red' ) ], + bcolor='white', + border=True, + size=( 0.15, 0.1 ), + loc='upper left' ) + + # Set camera and lighting + plotter.camera_position = 'iso' + plotter.add_light( pv.Light( position=( 1, 1, 1 ), intensity=0.8 ) ) + + # Add title + pathLength = np.sum( np.sqrt( np.sum( np.diff( pathPoints, axis=0 )**2, axis=1 ) ) ) + depthRange = pathZ.max() - pathZ.min() + title = 'Profile Path Extraction\n' + title += f'Points: {len(pathX)} | Length: {pathLength:.1f}m | Depth range: {depthRange:.1f}m' + plotter.add_text( title, position='upper_edge', font_size=14, color='black' ) + + # Save screenshot + # if savePath is not None: + # screenshot_path = savePath / 'profile_path_3d.png' + # plotter.screenshot(str(screenshot_path)) + # print(f" 💾 Screenshot saved: {screenshot_path}") + + # Show + if show: + plotter.show() + else: + plotter.close() diff --git a/geos-processing/src/geos/processing/post_processing/SensitivityAnalyzer.py b/geos-processing/src/geos/processing/post_processing/SensitivityAnalyzer.py new file mode 100644 index 000000000..7fa15e026 --- /dev/null +++ b/geos-processing/src/geos/processing/post_processing/SensitivityAnalyzer.py @@ -0,0 +1,299 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2026 TotalEnergies. +# SPDX-FileContributor: Nicolas Pillardou, Paloma Martinez +# ============================================================================ +# SENSITIVITY ANALYSIS +# ============================================================================ +import pandas as pd +from pathlib import Path +import numpy as np +from matplotlib.colors import Normalize +from matplotlib.cm import ScalarMappable +import matplotlib.pyplot as plt +from typing_extensions import Any, Self + + +from vtkmodules.vtkCommonDataModel import vtkCellData, vtkDataSet, vtkPolyData, vtkUnstructuredGrid + +from geos.utils.pieceEnum import Piece +from geos.mesh.utils.arrayHelpers import ( getArrayInObject ) +from geos.processing.post_processing.MohrCoulomb import ( MohrCoulomb ) + +from geos.processing.post_processing.ProfileExtractor import ( ProfileExtractor ) + + +class SensitivityAnalyzer: + """Performs sensitivity analysis on Mohr-Coulomb parameters.""" + + # ------------------------------------------------------------------- + def __init__( self: Self, outputDir: str = ".", showPlots: bool = True ) -> None: + """Init.""" + self.outputDir = Path( outputDir ) + self.outputDir.mkdir( exist_ok=True ) + self.results: list[ dict[ str, Any ] ] = [] + self.showPlots = showPlots + + # ------------------------------------------------------------------- + def runAnalysis( self: Self, surfaceWithStress, time: float, sensitivityFrictionAngles: list[float], sensitivityCohesions: list[float], profileStartPoints: list[tuple[float]], profileSearchRadius: list[tuple[float]] ) -> list[ dict[ str, Any ] ]: + """Run sensitivity analysis for multiple friction angles and cohesions.""" + frictionAngles = sensitivityFrictionAngles + cohesions = sensitivityCohesions + + print( "\n" + "=" * 60 ) + print( "SENSITIVITY ANALYSIS" ) + print( "=" * 60 ) + print( f"Friction angles: {frictionAngles}" ) + print( f"Cohesions: {cohesions}" ) + print( f"Total combinations: {len(frictionAngles) * len(cohesions)}" ) + + results = [] + for frictionAngle in frictionAngles: + for cohesion in cohesions: + print( f"\n→ Testing φ={frictionAngle}°, C={cohesion} bar" ) + + surfaceCopy = type(surfaceWithStress)() + surfaceCopy.DeepCopy( surfaceWithStress ) + + surfaceAnalyzed = MohrCoulomb.analyze( + surfaceCopy, + cohesion, + frictionAngle, + verbose=False ) + + stats = self._extractStatistics( surfaceAnalyzed, frictionAngle, cohesion ) + results.append( stats ) + + print( f" Unstable: {stats['nUnstable']}, " + f"Critical: {stats['nCritical']}, " + f"Stable: {stats['nStable']}" ) + + self.results = results + + # Generate plots + self._plotSensitivityResults( results, time ) + + # Plot SCU vs depth + self._plotSCUDepthProfiles( results, time, surfaceWithStress, profileStartPoints, profileSearchRadius ) + + return results + + # ------------------------------------------------------------------- + def _extractStatistics( self: Self, surface: vtkPolyData, frictionAngle: float, + cohesion: float ) -> dict[ str, Any ]: + """Extract statistical metrics from analyzed surface.""" + stability = getArrayInObject( surface, "stabilityState", Piece.CELLS ) + SCU = getArrayInObject( surface, "SCU", Piece.CELLS ) + failureProba = getArrayInObject( surface, "failureProbability", Piece.CELLS ) + safetyMargin = getArrayInObject( surface, "safetyMargin", Piece.CELLS ) + + nCells = surface.GetNumberOfCells() + + + stats = { + 'frictionAngle': frictionAngle, + 'cohesion': cohesion, + 'nCells': nCells, + 'nStable': np.sum( stability == 0 ), + 'nCritical': np.sum( stability == 1 ), + 'nUnstable': np.sum( stability == 2 ), + 'pctUnstable': np.sum( stability == 2 ) / nCells * 100, + 'pctCritical': np.sum( stability == 1 ) / nCells * 100, + 'pctStable': np.sum( stability == 0 ) / nCells * 100, + 'meanSCU': np.mean( SCU ), + 'maxSCU': np.max( SCU ), + 'meanFailureProb': np.mean( failureProba ), + 'meanSafetyMargin': np.mean( safetyMargin ), + 'minSafetyMargin': np.min( safetyMargin ) + } + + return stats + + # ------------------------------------------------------------------- + def _plotSensitivityResults( self: Self, results: list[ dict[ str, Any ] ], time: float ) -> None: + """Create comprehensive sensitivity analysis plots.""" + df = pd.DataFrame( results ) + + fig, axes = plt.subplots( 2, 2, figsize=( 16, 12 ) ) + + # Plot heatmaps + self._plotHeatMap( df, 'pctUnstable', 'Unstable Cells [%]', axes[ 0, 0 ] ) + self._plotHeatMap( df, 'pctCritical', 'Critical Cells [%]', axes[ 0, 1 ] ) + self._plotHeatMap( df, 'meanSCU', 'Mean SCU [-]', axes[ 1, 0 ] ) + self._plotHeatMap( df, 'meanSafetyMargin', 'Mean Safety Margin [bar]', axes[ 1, 1 ] ) + + fig.tight_layout() + + years = time / ( 365.25 * 24 * 3600 ) + filename = f'sensitivity_analysis_{years:.0f}y.png' + fig.savefig( self.outputDir / filename, dpi=300, bbox_inches='tight' ) + print( f"\n📊 Sensitivity plot saved: {filename}" ) + + if self.showPlots: + fig.show() + + # ------------------------------------------------------------------- + def _plotHeatMap( self: Self, df: pd.DataFrame, column: str, title: str, ax: plt.Axes ) -> None: + """Create a single heatmap for sensitivity analysis.""" + pivot = df.pivot( index='cohesion', columns='frictionAngle', values=column ) + + im = ax.imshow( pivot.values, cmap='RdYlGn_r', aspect='auto', origin='lower' ) + + ax.set_xticks( np.arange( len( pivot.columns ) ) ) + ax.set_yticks( np.arange( len( pivot.index ) ) ) + ax.set_xticklabels( pivot.columns ) + ax.set_yticklabels( pivot.index ) + + ax.set_xlabel( 'Friction Angle [°]' ) + ax.set_ylabel( 'Cohesion [bar]' ) + ax.set_title( title ) + + # Add values in cells + for i in range( len( pivot.index ) ): + for j in range( len( pivot.columns ) ): + value = pivot.values[ i, j ] + textColor = 'white' if value > pivot.values.max() * 0.5 else 'black' + ax.text( j, i, f'{value:.1f}', ha='center', va='center', color=textColor, fontsize=9 ) + + plt.colorbar( im, ax=ax ) + + # ------------------------------------------------------------------- + def _plotSCUDepthProfiles( self: Self, results: list[ dict[ str, Any ] ], time: float, + surfaceWithStress: vtkDataSet, profileStartPoints=None, profileSearchRadius=None, + maxDepthProfiles=None ) -> None: + """Plot SCU depth profiles for all parameter combinations. + + Each (cohesion, friction) pair gets a unique color + """ + print( "\n 📊 Creating SCU sensitivity depth profiles..." ) + + # Extract depth data + centers = getArrayInObject( surfaceWithStress, 'elementCenter', Piece.CELLS ) + centers[ :, 2 ] + + # Get profile points from config + + # Auto-generate if not provided + if profileStartPoints is None: + print( " ⚠️ No PROFILE_START_POINTS in config, auto-generating..." ) + xMin, xMax = np.min( centers[ :, 0 ] ), np.max( centers[ :, 0 ] ) + yMin, yMax = np.min( centers[ :, 1 ] ), np.max( centers[ :, 1 ] ) + + xRange = xMax - xMin + yRange = yMax - yMin + + if xRange > yRange: + # Fault oriented in X, sample at mid-Y + xPos = ( xMin + xMax ) / 2 + yPos = ( yMin + yMax ) / 2 + else: + # Fault oriented in Y, sample at mid-X + xPos = ( xMin + xMax ) / 2 + yPos = ( yMin + yMax ) / 2 + + profileStartPoints = [ ( xPos, yPos ) ] + + # Get search radius from config or auto-compute + if profileSearchRadius is None: + xMin, xMax = np.min( centers[ :, 0 ] ), np.max( centers[ :, 0 ] ) + yMin, yMax = np.min( centers[ :, 1 ] ), np.max( centers[ :, 1 ] ) + xRange = xMax - xMin + yRange = yMax - yMin + searchRadius = min( xRange, yRange ) * 0.15 + else: + searchRadius = profileSearchRadius + + print( f" 📍 Using {len(profileStartPoints)} profile point(s) from config" ) + print( f" Search radius: {searchRadius:.1f}m" ) + + # Create colormap for parameter combinations + nCombinations = len( results ) + cmap = plt.cm.viridis # type: ignore [attr-defined] + norm = Normalize( vmin=0, vmax=nCombinations - 1 ) + ScalarMappable( norm=norm, cmap=cmap ) + + # Create figure with subplots for each profile point + nProfiles = len( profileStartPoints ) + fig, axes = plt.subplots( 1, nProfiles, figsize=( 8 * nProfiles, 10 ) ) + + # Handle single subplot case + if nProfiles == 1: + axes = [ axes ] + + # Plot each profile point + for profileIdx, ( xPos, yPos, zPos ) in enumerate( profileStartPoints ): + ax = axes[ profileIdx ] + + print( f"\n → Profile {profileIdx+1} at ({xPos:.1f}, {yPos:.1f}, {zPos:.1f}):" ) + + # Plot each parameter combination + for idx, params in enumerate( results ): + frictionAngle = params[ 'frictionAngle' ] + cohesion = params[ 'cohesion' ] + + # Re-analyze surface with these parameters + surfaceCopy = type(surfaceWithStress)() + surfaceCopy.DeepCopy( surfaceWithStress ) + surfaceAnalyzed = MohrCoulomb.analyze( + surfaceCopy, + cohesion, + frictionAngle, + verbose=False ) + + # Extract SCU + SCU = np.abs( getArrayInObject( surfaceAnalyzed, "SCU", Piece.CELLS ) ) + + # Extract profile using adaptive method + # depthsSCU, profileSCU, _, _ = ProfileExtractor.extractVerticalProfileTopologyBased( + # surfaceAnalyzed, 'SCU', xPos, yPos, zPos, verbose=False) + depthsSCU, profileSCU, _, _ = ProfileExtractor.extractAdaptiveProfile( + centers, SCU, xPos, yPos, searchRadius ) + + if len( depthsSCU ) >= 3: + color = cmap( norm( idx ) ) + label = f'φ={frictionAngle}°, C={cohesion} bar' + ax.plot( profileSCU, depthsSCU, color=color, label=label, linewidth=2, alpha=0.8 ) + + if idx == 0: # Print info only once per profile + print( f" ✅ {len(depthsSCU)} points extracted" ) + else: + if idx == 0: + print( f" ⚠️ Insufficient points ({len(depthsSCU)})" ) + + # Add critical lines + ax.axvline( x=0.8, + color='forestgreen', + linestyle='--', + linewidth=2, + label='Critical (SCU=0.8)', + zorder=100 ) + ax.axvline( x=1.0, color='red', linestyle='--', linewidth=2, label='Failure (SCU=1.0)', zorder=100 ) + + # Configure plot + ax.set_xlabel( 'Shear Capacity Utilization (SCU) [-]', fontsize=14, weight='bold' ) + ax.set_ylabel( 'Depth [m]', fontsize=14, weight='bold' ) + ax.set_title( f'Profile {profileIdx+1} @ ({xPos:.0f}, {yPos:.0f})', fontsize=14, weight='bold' ) + ax.grid( True, alpha=0.3, linestyle='--' ) + ax.set_xlim( left=0 ) + + # Change verticale scale + if maxDepthProfiles is not None: + ax.set_ylim( bottom=maxDepthProfiles ) + + # Légende en dehors à droite + ax.legend( loc='center left', bbox_to_anchor=( 1, 0.5 ), fontsize=9, ncol=1 ) + + ax.tick_params( labelsize=12 ) + + # Overall title + years = time / ( 365.25 * 24 * 3600 ) + fig.suptitle( 'SCU Depth Profiles - Sensitivity Analysis', fontsize=16, weight='bold', y=0.98 ) + + fig.tight_layout( rect=( 0, 0, 1, 0.96 ) ) + + # Save + filename = f'sensitivity_scu_profiles_{years:.0f}y.png' + fig.savefig( self.outputDir / filename, dpi=300, bbox_inches='tight' ) + print( f"\n 💾 SCU sensitivity profiles saved: {filename}" ) + + if self.showPlots: + fig.show() diff --git a/geos-processing/src/geos/processing/post_processing/StressProjector.py b/geos-processing/src/geos/processing/post_processing/StressProjector.py new file mode 100644 index 000000000..7709ce749 --- /dev/null +++ b/geos-processing/src/geos/processing/post_processing/StressProjector.py @@ -0,0 +1,688 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2026 TotalEnergies. +# SPDX-FileContributor: Nicolas Pillardou, Paloma Martinez +from pathlib import Path +import numpy as np +import numpy.typing as npt +from typing_extensions import Self, Any +from scipy.spatial import cKDTree +from xml.etree.ElementTree import ElementTree, Element, SubElement +import pyvista as pv + +from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkPolyData ) +from vtkmodules.util.numpy_support import vtk_to_numpy, numpy_to_vtk + +from geos.geomechanics.model.StressTensor import StressTensor +from geos.mesh.utils.genericHelpers import ( extractCellSelection, getLocalBasisVectors ) +from geos.mesh.utils.arrayHelpers import (isAttributeInObject, getArrayInObject, computeCellCenterCoordinates) +from geos.mesh.utils.arrayModifiers import (createAttribute, updateAttribute) +from geos.utils.pieceEnum import Piece +from geos.mesh.io.vtkIO import writeMesh, VtkOutput + + +# ============================================================================ +# STRESS PROJECTION +# ============================================================================ +class StressProjector: + """Projects volume stress onto fault surfaces and tracks principal stresses in VTU.""" + + # ------------------------------------------------------------------- + def __init__( self: Self, adjacencyMapping: dict[ int, list[ pv.DataSet ] ], + geometricProperties: dict[ str, Any ], + outputDir: str = ".", ) -> None: + """Initialize with pre-computed adjacency mapping and geometric properties. + + Parameters + ---------- + config : Configuration object + adjacencyMapping : dict + Pre-computed dict mapping fault cells to volume cells + geometricProperties : dict + Pre-computed geometric properties from FaultGeometry: + - 'volumes': cell volumes + - 'centers': cell centers + - 'distances': distances to fault + - 'faultTree': KDTree for fault + """ + # self.config = config + self.adjacencyMapping = adjacencyMapping + + # Store pre-computed geometric properties + self.volumeCellVolumes = geometricProperties[ 'volumes' ] + self.volumeCenters = geometricProperties[ 'centers' ] + self.distanceToFault = geometricProperties[ 'distances' ] + self.faultTree = geometricProperties[ 'faultTree' ] + + # Storage for time series metadata + self.timestepInfo: list[ dict[ str, Any ] ] = [] + + # Track which cells to monitor (optional) + self.monitoredCells: set[ int ] | None = None + + # Output directory for VTU files + self.vtuOutputDir = Path( outputDir ) / "principal_stresses" + + # ------------------------------------------------------------------- + def setMonitoredCells( self: Self, cellIndices: list[ int ] | None = None ) -> None: + """Set specific cells to monitor (optional). + + Parameters: + cellIndices: list of volume cell indices to track + If None, all contributing cells are tracked + """ + self.monitoredCells = set( cellIndices ) if cellIndices is not None else None + + # ------------------------------------------------------------------- + def projectStressToFault( + self: Self, + volumeData: vtkUnstructuredGrid, + volumeInitial: vtkUnstructuredGrid, + faultSurface: vtkPolyData, + time: float | None = None, + timestep: int | None = None, + weightingScheme: str = "arithmetic", + stressName: str = "averageStress", + biotName: str = "rockPorosity_biotCoefficient", + computePrincipalStresses: bool = False, + frictionAngle: float = 10, cohesion: float = 0 ) -> tuple[ vtkPolyData, vtkUnstructuredGrid, vtkUnstructuredGrid ]: + """Project stress and save principal stresses to VTU. + + Now uses pre-computed geometric properties for efficiency + """ + if not isAttributeInObject ( volumeData, stressName, Piece.CELLS ): + raise ValueError( f"No stress data '{stressName}' in dataset" ) + + # ===================================================================== + # 1. EXTRACT STRESS DATA + # ===================================================================== + pressure = getArrayInObject( volumeData, "pressure", Piece.CELLS ) / 1e5 + pressureFault = getArrayInObject( volumeInitial, "pressure", Piece.CELLS ) / 1e5 + pressureInitial = getArrayInObject( volumeInitial, "pressure", Piece.CELLS ) / 1e5 + biot = getArrayInObject( volumeData, biotName, Piece.CELLS ) + + stressEffective = StressTensor.buildFromArray( getArrayInObject( volumeData, stressName, Piece.CELLS ) / 1e5 ) + stressEffectiveInitial = StressTensor.buildFromArray( getArrayInObject( volumeInitial, stressName, Piece.CELLS ) / 1e5 ) + + # Convert effective stress to total stress + arrI = np.eye( 3 )[ None, :, : ] + stressTotal = stressEffective - biot[ :, None, None ] * pressure[ :, None, None ] * arrI + stressTotalInitial = stressEffectiveInitial - biot[ :, None, None ] * pressureInitial[ :, None, None ] * arrI + + # ===================================================================== + # 2. USE PRE-COMPUTED ADJACENCY + # ===================================================================== + + # ===================================================================== + # 3. PREPARE FAULT GEOMETRY + # ===================================================================== + # TODO fix + normalsXX, tangent1XX, tangent2XX = getLocalBasisVectors( faultSurface ) + # normals = faultSurface.cell_data[ "Normals" ] + # tangent1 = faultSurface.cell_data[ "tangent1" ] + # tangent2 = faultSurface.cell_data[ "tangent2" ] + normals = vtk_to_numpy( faultSurface.GetCellData().GetNormals() ) + tangent1 = vtk_to_numpy( faultSurface.GetCellData().GetArray( "Tangents1") ) + tangent2 = vtk_to_numpy( faultSurface.GetCellData().GetArray( "Tangents2") ) + + print( (normalsXX - normals ).max() ) + print( (tangent1XX - tangent1 ).max() ) + print( (tangent2XX - tangent2 ).max() ) + + faultCenters = computeCellCenterCoordinates( faultSurface ) + fcenters = faultSurface.GetCellData().GetArray( 'elementCenter' ) + fcenters = faultCenters + + nFault = faultSurface.GetNumberOfCells() + + # ===================================================================== + # 4. COMPUTE PRINCIPAL STRESSES FOR CONTRIBUTING CELLS + # ===================================================================== + if computePrincipalStresses and timestep is not None: + + # Collect all unique contributing cells + allContributingCells = set() + for _faultIdx, neighbors in self.adjacencyMapping.items(): + allContributingCells.update( neighbors[ 'plus' ] ) + allContributingCells.update( neighbors[ 'minus' ] ) + + # Filter by monitored cells if specified + if self.monitoredCells is not None: + cellsToTrack = allContributingCells.intersection( self.monitoredCells ) + else: + cellsToTrack = allContributingCells + + print( f" 📊 Computing principal stresses for {len(cellsToTrack)} contributing cells..." ) + + # Create mesh with only contributing cells + contributingMesh = self._createVolumicContribMesh( volumeData, faultSurface, cellsToTrack, + self.adjacencyMapping, biotName=biotName, computePrincipalStresses=computePrincipalStresses, + frictionAngle=frictionAngle, cohesion=cohesion ) + + self._savePrincipalStressVTU( contributingMesh, time, timestep ) + + else: + contributingMesh = None + + # ===================================================================== + # 6. PROJECT STRESS FOR EACH FAULT CELL + # ===================================================================== + sigmaNArr = np.zeros( nFault ) + tauArr = np.zeros( nFault ) + tauDipArr = np.zeros( nFault ) + tauStrikeArr = np.zeros( nFault ) + deltaSigmaNArr = np.zeros( nFault ) + deltaTauArr = np.zeros( nFault ) + nContributors = np.zeros( nFault, dtype=int ) + + print( f" 🔄 Projecting stress to {nFault} fault cells..." ) + print( f" Weighting scheme: {weightingScheme}" ) + + for faultIdx in range( nFault ): + if faultIdx not in self.adjacencyMapping: + continue + + volPlus = self.adjacencyMapping[ faultIdx ][ 'plus' ] + volMinus = self.adjacencyMapping[ faultIdx ][ 'minus' ] + allVol = volPlus + volMinus + + if len( allVol ) == 0: + continue + + # =================================================================== + # CALCULATE WEIGHTS (using pre-computed properties) + # =================================================================== + + if weightingScheme == 'arithmetic' or weightingScheme == 'harmonic': + weights = np.ones( len( allVol ) ) / len( allVol ) + + elif weightingScheme == 'distance': + # Use pre-computed distances + dists = np.array( [ self.distanceToFault[ v ] for v in allVol ] ) + dists = np.maximum( dists, 1e-6 ) + invDists = 1.0 / dists + weights = invDists / np.sum( invDists ) + + elif weightingScheme == 'volume': + # Use pre-computed volumes + vols = np.array( [ self.volumeCellVolumes[ v ] for v in allVol ] ) + weights = vols / np.sum( vols ) + + elif weightingScheme == 'distanceVolume': + # Use pre-computed volumes and distances + vols = np.array( [ self.volumeCellVolumes[ v ] for v in allVol ] ) + dists = np.array( [ self.distanceToFault[ v ] for v in allVol ] ) + dists = np.maximum( dists, 1e-6 ) + + weights = vols / dists + weights = weights / np.sum( weights ) + + elif weightingScheme == 'inverseSquareDistance': + # Use pre-computed distances + dists = np.array( [ self.distanceToFault[ v ] for v in allVol ] ) + dists = np.maximum( dists, 1e-6 ) + invSquareDistance = 1.0 / ( dists**2 ) + weights = invSquareDistance / np.sum( invSquareDistance ) + + else: + raise ValueError( f"Unknown weighting scheme: {weightingScheme}" ) + + # =================================================================== + # ACCUMULATE WEIGHTED CONTRIBUTIONS + # =================================================================== + + sigmaN = 0.0 + tau = 0.0 + tauDip = 0.0 + tauStrike = 0.0 + deltaSigmaN = 0.0 + deltaTau = 0.0 + + for volIdx, w in zip( allVol, weights ): + + # Total stress (with pressure) + sigmaFinal = stressTotal[ volIdx ] + pressureFault[ volIdx ] * np.eye( 3 ) + sigmaInit = stressTotalInitial[ volIdx ] + pressureInitial[ volIdx ] * np.eye( 3 ) + + # Rotate to fault frame + resFinal = StressTensor.rotateToFaultFrame( sigmaFinal, normals[ faultIdx ], tangent1[ faultIdx ], + tangent2[ faultIdx ] ) + + resInitial = StressTensor.rotateToFaultFrame( sigmaInit, normals[ faultIdx ], tangent1[ faultIdx ], + tangent2[ faultIdx ] ) + + # Accumulate weighted contributions + sigmaN += w * resFinal[ 'normalStress' ] + tau += w * resFinal[ 'shearStress' ] + tauDip += w * resFinal[ 'shearDip' ] + tauStrike += w * resFinal[ 'shearStrike' ] + deltaSigmaN += w * ( resFinal[ 'normalStress' ] - resInitial[ 'normalStress' ] ) + deltaTau += w * ( resFinal[ 'shearStress' ] - resInitial[ 'shearStress' ] ) + + sigmaNArr[ faultIdx ] = sigmaN + tauArr[ faultIdx ] = tau + tauDipArr[ faultIdx ] = tauDip + tauStrikeArr[ faultIdx ] = tauStrike + deltaSigmaNArr[ faultIdx ] = deltaSigmaN + deltaTauArr[ faultIdx ] = deltaTau + nContributors[ faultIdx ] = len( allVol ) + + # ===================================================================== + # 7. STORE RESULTS ON FAULT SURFACE + # ===================================================================== + for attributeName, value in zip( + [ "sigmaNEffective", "tauEffective", "tauStrike", "tauDip", "deltaSigmaNEffective", "deltaTauEffective" ], + [ sigmaNArr, tauDipArr, tauStrikeArr, tauDipArr, deltaSigmaNArr, deltaTauArr ] + ): + updateAttribute( faultSurface, value, attributeName, Piece.CELLS ) + + # ===================================================================== + # 8. STATISTICS + # ===================================================================== + valid = nContributors > 0 + nValid = np.sum( valid ) + + print( f" ✅ Stress projected: {nValid}/{nFault} fault cells ({nValid/nFault*100:.1f}%)" ) + + if np.sum( valid ) > 0: + print( f" Contributors per fault cell: min={np.min(nContributors[valid])}, " + f"max={np.max(nContributors[valid])}, " + f"mean={np.mean(nContributors[valid]):.1f}" ) + + return faultSurface, volumeData, contributingMesh + + # ------------------------------------------------------------------- + @staticmethod + def computePrincipalStresses( stressTensor: StressTensor ) -> dict[ str, npt.NDArray[ np.float64 ] ]: + """Compute principal stresses and directions. + + Convention: Compression is NEGATIVE + - σ1 = most compressive (most negative) + - σ3 = least compressive (least negative, or most tensile) + + Returns: + dict with eigenvalues, eigenvectors, meanStress, deviatoricStress + """ + eigenvalues, eigenvectors = np.linalg.eigh( stressTensor ) + + # Sort from MOST NEGATIVE to LEAST NEGATIVE (most compressive to least) + # Example: -600 < -450 < -200, so -600 is σ1 (most compressive) + idx = np.argsort( eigenvalues ) # Ascending order (most negative first) + eigenvaluesSorted = eigenvalues[ idx ] + eigenVectorsSorted = eigenvectors[ :, idx ] + + return { + 'sigma1': eigenvaluesSorted[ 0 ], # Most compressive (most negative) + 'sigma2': eigenvaluesSorted[ 1 ], # Intermediate + 'sigma3': eigenvaluesSorted[ 2 ], # Least compressive (least negative) + 'meanStress': np.mean( eigenvaluesSorted ), + 'deviatoricStress': eigenvaluesSorted[ 0 ] - + eigenvaluesSorted[ 2 ], # σ1 - σ3 (negative - more negative = positive or less negative) + 'direction1': eigenVectorsSorted[ :, 0 ], # Direction of σ1 + 'direction2': eigenVectorsSorted[ :, 1 ], # Direction of σ2 + 'direction3': eigenVectorsSorted[ :, 2 ] # Direction of σ3 + } + + # ------------------------------------------------------------------- + def _createVolumicContribMesh( self: Self, volumeData: pv.UnstructuredGrid, faultSurface: pv.PolyData, + cellsToTrack: set[ int ], mapping: dict[ int, list[ pv.DataSet ] ], biotName: str = "rockPorosity_biotCoefficient", + computePrincipalStresses: bool = False, frictionAngle : float = 10, cohesion: float = 0 ) -> pv.DataSet: + """Create a mesh containing only contributing cells with principal stress data and compute analytical normal/shear stresses based on fault dip angle. + + Parameters + ---------- + volumeData : pyvista.UnstructuredGrid + Volume mesh with stress data (rock_stress or averageStress) + faultSurface : pyvista.PolyData + Fault surface with dipAngle and strikeAngle per cell + cellsToTrack : set + Set of volume cell indices to include + mapping : dict + Adjacency mapping {faultIdx: {'plus': [...], 'minus': [...]}} + """ + # =================================================================== + # EXTRACT STRESS DATA FROM VOLUME + # =================================================================== + + if not isAttributeInObject( volumeData, stressName, Piece.CELLS ): + raise ValueError( f"No stress data '{stressName}' in volume dataset" ) + + print( f" 📊 Extracting stress from field: '{stressName}'" ) + + # Extract effective stress and pressure + pressure = getArrayInObject( volumeData, "pressure", Piece.CELLS ) / 1e5 # Convert to bar + biot = getArrayInObject( volumeData, biotName, Piece.CELLS ) + + stressEffective = StressTensor.buildFromArray( getArrayInObject( volumeData, stressName, Piece.CELLS ) / 1e5 ) + + # Convert effective stress to total stress + arrI = np.eye( 3 )[ None, :, : ] + stressTotal = stressEffective - biot[ :, None, None ] * pressure[ :, None, None ] * arrI + + # =================================================================== + # EXTRACT SUBSET OF CELLS + # =================================================================== + cellIndices = sorted( cellsToTrack ) + cellMask = np.zeros( volumeData.GetNumberOfCells(), dtype=bool ) + cellMask[ cellIndices ] = True + + subsetMesh = extractCellSelection ( cellMask ) + + # =================================================================== + # REBUILD MAPPING: subsetIdx -> originalIdx + # =================================================================== + originalCenters = vtk_to_numpy( computeCellCenterCoordinates( volumeData ) )[ cellIndices ] + subsetCenters = vtk_to_numpy( computeCellCenterCoordinates( subsetMesh ) ) + + tree = cKDTree( originalCenters ) + + subsetToOriginal = np.zeros( subsetMesh.GetNumberOfCells(), dtype=int ) + for subsetIdx in range( subsetMesh.GetNumberOfCells() ): + dist, idx = tree.query( subsetCenters[ subsetIdx ] ) + if dist > 1e-6: + print( f" WARNING: Cell {subsetIdx} not matched (dist={dist})" ) + subsetToOriginal[ subsetIdx ] = cellIndices[ idx ] + + # =================================================================== + # MAP VOLUME CELLS TO FAULT DIP/STRIKE ANGLES + # =================================================================== + print( " 📐 Mapping volume cells to fault dip/strike angles..." ) + + # Check if fault surface has required data + if 'dipAngle' not in faultSurface.cell_data: + print( " ⚠️ WARNING: 'dipAngle' not found in faultSurface" ) + print( f" Available fields: {list(faultSurface.cell_data.keys())}" ) + return None + + if 'strikeAngle' not in faultSurface.cell_data: + print( " ⚠️ WARNING: 'strikeAngle' not found in faultSurface" ) + + # Create mapping: volume_cell_id -> [dipAngles, strikeAngles] + volumeToDip: dict[ int, npt.NDArray[ np.float64 ] ] = {} + volumeToStrike: dict[ int, npt.NDArray[ np.float64 ] ] = {} + + for faultIdx, neighbors in mapping.items(): + # Get dip and strike angle from fault cell + faultDip = getArrayInObject( faultSurface, 'dipAngle', Piece.CELLS )[ faultIdx ] + + # Strike is optional + if isAttributeInObject ( faultSurface, 'strikeAngle', Piece.CELLS): + faultStrike = getArrayInObject( faultSurface, 'strikeAngle', Piece.CELLS )[ faultIdx ] + else: + faultStrike = np.nan + + # Assign to all contributing volume cells (plus and minus) + for volIdx in neighbors[ 'plus' ] + neighbors[ 'minus' ]: + if volIdx not in volumeToDip: + volumeToDip[ volIdx ] = [] + volumeToStrike[ volIdx ] = [] + volumeToDip[ volIdx ].append( faultDip ) + volumeToStrike[ volIdx ].append( faultStrike ) + + # Average if a volume cell contributes to multiple fault cells + volumeToDipAvg = { volIdx: np.mean( dips ) for volIdx, dips in volumeToDip.items() } + volumeToStrikeAvg = { volIdx: np.mean( strikes ) for volIdx, strikes in volumeToStrike.items() } + + print( f" ✅ Mapped {len(volumeToDipAvg)} volume cells to fault angles" ) + + # Statistics + allDips = [ np.mean( dips ) for dips in volumeToDip.values() ] + if len( allDips ) > 0: + print( f" Dip angle range: [{np.min(allDips):.1f}, {np.max(allDips):.1f}]°" ) + + # =================================================================== + # COMPUTE PRINCIPAL STRESSES AND ANALYTICAL FAULT STRESSES + # =================================================================== + nCells = subsetMesh.GetNumberOfCells() + + sigma1Arr = np.zeros( nCells ) + sigma2Arr = np.zeros( nCells ) + sigma3Arr = np.zeros( nCells ) + meanStressArr = np.zeros( nCells ) + deviatoricStressArr = np.zeros( nCells ) + pressureArr = np.zeros( nCells ) + + direction1Arr = np.zeros( ( nCells, 3 ) ) + direction2Arr = np.zeros( ( nCells, 3 ) ) + direction3Arr = np.zeros( ( nCells, 3 ) ) + + # NEW: Analytical fault stresses + sigmaNAnalyticalArr = np.zeros( nCells ) + tauAnalyticalArr = np.zeros( nCells ) + dipAngleArr = np.zeros( nCells ) + strikeAngleArr = np.zeros( nCells ) + deltaArr = np.zeros( nCells ) + + sideArr = np.zeros( nCells, dtype=int ) + nFaultCellsArr = np.zeros( nCells, dtype=int ) + + print( " 🔢 Computing principal stresses and analytical projections..." ) + + for subsetIdx in range( nCells ): + origIdx = subsetToOriginal[ subsetIdx ] + + # =============================================================== + # COMPUTE PRINCIPAL STRESSES + # =============================================================== + # Total stress = effective stress + pore pressure + sigmaTotalCell = stressTotal[ origIdx ] + pressure[ origIdx ] * np.eye( 3 ) + principal = self.computePrincipalStresses( sigmaTotalCell ) + + sigma1Arr[ subsetIdx ] = principal[ 'sigma1' ] + sigma2Arr[ subsetIdx ] = principal[ 'sigma2' ] + sigma3Arr[ subsetIdx ] = principal[ 'sigma3' ] + meanStressArr[ subsetIdx ] = principal[ 'meanStress' ] + deviatoricStressArr[ subsetIdx ] = principal[ 'deviatoricStress' ] + pressureArr[ subsetIdx ] = pressure[ origIdx ] + + direction1Arr[ subsetIdx ] = principal[ 'direction1' ] + direction2Arr[ subsetIdx ] = principal[ 'direction2' ] + direction3Arr[ subsetIdx ] = principal[ 'direction3' ] + + # =============================================================== + # COMPUTE ANALYTICAL FAULT STRESSES (Anderson formulas) + # =============================================================== + if origIdx in volumeToDipAvg: + dipDeg = volumeToDipAvg[ origIdx ] + dipAngleArr[ subsetIdx ] = dipDeg + + strikeDeg = volumeToStrikeAvg.get( origIdx, np.nan ) + strikeAngleArr[ subsetIdx ] = strikeDeg + + # δ = 90° - dip (angle from horizontal) + deltaDeg = 90.0 - dipDeg + deltaRad = np.radians( deltaDeg ) + deltaArr[ subsetIdx ] = deltaDeg + + # Extract principal stresses (compression negative) + sigma1 = principal[ 'sigma1' ] # Most compressive (most negative) + sigma3 = principal[ 'sigma3' ] # Least compressive (least negative) + + # Anderson formulas (1951) + # σ_n = (σ1 + σ3)/2 - (σ1 - σ3)/2 * cos(2δ) + # τ = |(σ1 - σ3)/2 * sin(2δ)| + + sigmaMean = ( sigma1 + sigma3 ) / 2.0 + sigmaDiff = ( sigma1 - sigma3 ) / 2.0 + + sigmaNAnalytical = sigmaMean - sigmaDiff * np.cos( 2 * deltaRad ) + tauAnalytical = sigmaDiff * np.sin( 2 * deltaRad ) + + sigmaNAnalyticalArr[ subsetIdx ] = sigmaNAnalytical + tauAnalyticalArr[ subsetIdx ] = np.abs( tauAnalytical ) + else: + # No fault association - set to NaN + dipAngleArr[ subsetIdx ] = np.nan + strikeAngleArr[ subsetIdx ] = np.nan + deltaArr[ subsetIdx ] = np.nan + sigmaNAnalyticalArr[ subsetIdx ] = np.nan + tauAnalyticalArr[ subsetIdx ] = np.nan + + # =============================================================== + # DETERMINE SIDE (plus/minus/both) + # =============================================================== + isPlus = False + isMinus = False + faultCellCount = 0 + + for _faultIdx, neighbors in mapping.items(): + if origIdx in neighbors[ 'plus' ]: + isPlus = True + faultCellCount += 1 + if origIdx in neighbors[ 'minus' ]: + isMinus = True + faultCellCount += 1 + + if isPlus and isMinus: + sideArr[ subsetIdx ] = 3 # both + elif isPlus: + sideArr[ subsetIdx ] = 1 # plus + elif isMinus: + sideArr[ subsetIdx ] = 2 # minus + else: + sideArr[ subsetIdx ] = 0 # none (should not happen) + + nFaultCellsArr[ subsetIdx ] = faultCellCount + + # =================================================================== + # ADD DATA TO MESH + # =================================================================== + for attributeName, attributeArray in zip([ + ( 'sigma1' , sigma1Arr ), + ( 'sigma2' , sigma2Arr ), + ( 'sigma3' , sigma3Arr ), + ( 'meanStress' , meanStressArr ), + ( 'deviatoricStress' , deviatoricStressArr ), + ( 'pressure_bar' , pressureArr ), + ( 'sigma1Direction' , direction1Arr ), + ( 'sigma2Direction' , direction2Arr ), + ( 'sigma3Direction' , direction3Arr ), + ( 'sigmaNAnalytical', sigmaNAnalyticalArr ), + ( 'tauAnalytical', tauAnalyticalArr ), + ( 'dipAngle', dipAngleArr ), + ( 'strikeAngle', strikeAngleArr ), + ( 'deltaAngle', deltaArr ), + ]): + updateAttribute( subsetMesh, attributeArray, attributeName, piece=Piece.CELLS ) + + # =================================================================== + # COMPUTE SCU ANALYTICALLY (Mohr-Coulomb) + # =================================================================== + mu = np.tan( np.radians( frictionAngle ) ) + + # τ_crit = C - σ_n * μ + # Note: σ_n is negative (compression), so -σ_n * μ is positive + tauCriticalArr = cohesion - sigmaNAnalyticalArr * mu + # SCU = τ / τ_crit + SCUAnalyticalArr = np.divide( tauAnalyticalArr, + tauCriticalArr, + out=np.zeros_like( tauAnalyticalArr ), + where=tauCriticalArr != 0 ) + + # CFS (Coulomb Failure Stress) + CFSAnalyticalArr = tauAnalyticalArr - mu * ( -sigmaNAnalyticalArr ) + + for attributeName, attributeArray in zip([ + ( 'tauCriticalAnalytical', tauCriticalArr ), + ( 'SCUAnalytical', SCUAnalyticalArr ), + ( 'side', sideArr ), + ( 'nFaultCells', nFaultCellsArr ), + ( 'originalCellIds', subsetToOriginal ) + ]): + createAttribute( subsetMesh, attributeArray, attributeName, piece=Piece.CELLS ) + + # =================================================================== + # STATISTICS + # =================================================================== + validAnalytical = ~np.isnan( sigmaNAnalyticalArr ) + nValid = np.sum( validAnalytical ) + + if nValid > 0: + print( f" 📊 Analytical fault stresses computed for {nValid}/{nCells} cells" ) + print( + f" σ_n range: [{np.nanmin(sigmaNAnalyticalArr):.1f}, {np.nanmax(sigmaNAnalyticalArr):.1f}] bar" ) + print( f" τ range: [{np.nanmin(tauAnalyticalArr):.1f}, {np.nanmax(tauAnalyticalArr):.1f}] bar" ) + print( f" Dip angle range: [{np.nanmin(dipAngleArr):.1f}, {np.nanmax(dipAngleArr):.1f}]°" ) + + # if hasattr( self.config, 'FRICTION_ANGLE' ) and hasattr( self.config, 'COHESION' ): + print( + f" SCU range: [{np.nanmin(SCUAnalyticalArr[validAnalytical]):.2f}, {np.nanmax(SCUAnalyticalArr[validAnalytical]):.2f}]" + ) + nCritical = np.sum( ( SCUAnalyticalArr >= 0.8 ) & ( SCUAnalyticalArr < 1.0 ) ) + nUnstable = np.sum( SCUAnalyticalArr >= 1.0 ) + print( f" Critical cells (SCU≥0.8): {nCritical} ({nCritical/nValid*100:.1f}%)" ) + print( f" Unstable cells (SCU≥1.0): {nUnstable} ({nUnstable/nValid*100:.1f}%)" ) + else: + print( " ⚠️ No analytical stresses computed (no fault mapping)" ) + + return subsetMesh + + # ------------------------------------------------------------------- + def _savePrincipalStressVTU( self: Self, mesh: pv.DataSet, time: float, timestep: int ) -> None: + """Save principal stress mesh to VTU file. + + Parameters: + mesh: PyVista mesh with principal stress data + time: Simulation time + timestep: Timestep index + """ + # Create output directory + self.vtuOutputDir.mkdir( parents=True, exist_ok=True ) + + # Generate filename + vtuFilename = f"principal_stresses_{timestep:05d}.vtu" + vtuPath = self.vtuOutputDir / vtuFilename + + # Save mesh + writeMesh( mesh=mesh, vtkOutput=VtkOutput( vtuPath ) ) + + # Store metadata for PVD + self.timestepInfo.append( { + 'time': time if time is not None else timestep, + 'timestep': timestep, + 'file': vtuFilename + } ) + + print( f" 💾 Saved principal stresses: {vtuFilename}" ) + + # ------------------------------------------------------------------- + def savePVDCollection( self: Self, filename: str = "principal_stresses.pvd" ) -> None: + """Create PVD file for time series visualization in ParaView. + + Parameters: + filename: Name of PVD file + """ + if len( self.timestepInfo ) == 0: + print( "⚠️ No timestep data to save in PVD" ) + return + + pvdPath = self.vtuOutputDir / filename + + print( f"\n💾 Creating PVD collection: {pvdPath}" ) + print( f" Timesteps: {len(self.timestepInfo)}" ) + + # Create XML structure + root = Element( 'VTKFile' ) + root.set( 'type', 'Collection' ) + root.set( 'version', '0.1' ) + root.set( 'byte_order', 'LittleEndian' ) + + collection = SubElement( root, 'Collection' ) + + for info in self.timestepInfo: + dataset = SubElement( collection, 'DataSet' ) + dataset.set( 'timestep', str( info[ 'time' ] ) ) + dataset.set( 'group', '' ) + dataset.set( 'part', '0' ) + dataset.set( 'file', info[ 'file' ] ) + + # Write to file + tree = ElementTree( root ) + tree.write( str( pvdPath ), encoding='utf-8', xml_declaration=True ) + + print( " ✅ PVD file created successfully" ) + print( f" 📂 Output directory: {self.vtuOutputDir}" ) + print( "\n 🎨 To visualize in ParaView:" ) + print( f" 1. Open: {pvdPath}" ) + print( " 2. Apply" ) + print( " 3. Color by: sigma1, sigma2, sigma3, meanStress, etc." ) + print( " 4. Use 'side' filter to show plus/minus/both" ) diff --git a/geos-processing/src/geos/processing/tools/FaultVisualizer.py b/geos-processing/src/geos/processing/tools/FaultVisualizer.py new file mode 100644 index 000000000..9934b59ee --- /dev/null +++ b/geos-processing/src/geos/processing/tools/FaultVisualizer.py @@ -0,0 +1,1595 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2026 TotalEnergies. +# SPDX-FileContributor: Nicolas Pillardou, Paloma Martinez +import os +import pandas as pd +import numpy as np +import numpy.typing as npt +from pathlib import Path +from matplotlib.lines import Line2D +import matplotlib.pyplot as plt +import pyvista as pv +from typing_extensions import Self + +from vtkmodules.vtkCommonDataModel import vtkPolyData + +from geos.processing.post_processing.ProfileExtractor import ProfileExtractor +from geos.mesh.utils.arrayHelpers import ( isAttributeInObject, getArrayInObject ) +from geos.utils.pieceEnum import Piece + + +# ============================================================================ +# VISUALIZATION +# ============================================================================ +class Visualizer: + """Visualization utilities.""" + + # ------------------------------------------------------------------- + # def __init__( self, config: Config ) -> None: + def __init__( self, profileSearchRadius: float| None = None, + minDepthProfiles: float | None = None, + maxDepthProfiles: float | None = None, + showPlots: bool = True, savePlots:bool = True ) -> None: + """Init.""" + self.profileSearchRadius = profileSearchRadius + self.minDepthProfiles = minDepthProfiles + self.maxDepthProfiles = maxDepthProfiles + # self.config = config + self.savePlots = savePlots + self.showPlots = showPlots + + # ------------------------------------------------------------------- + @staticmethod + def plotMohrCoulombDiagram( surface: vtkPolyData, + time: float, + path: Path, + show: bool = True, + save: bool = True ) -> None: + """Create Mohr-Coulomb diagram with depth coloring.""" + sigmaN = - getArrayInObject( surface, "sigmaNEffective", Piece.CELLS ) + tau = np.abs( getArrayInObject( surface, "tauEffective", Piece.CELLS ) ) + SCU = np.abs( getArrayInObject( surface, "SCU", Piece.CELLS ) ) + depth = getArrayInObject( surface, 'elementCenter', Piece.CELLS )[ :, 2 ] + + cohesion = getArrayInObject( surface, "mohrCohesion", Piece.CELLS )[ 0 ] + mu = getArrayInObject( surface, "mohrFrictionCoefficient" , Piece.CELLS )[ 0 ] + phi = getArrayInObject( surface, 'mohrFrictionAngle' , Piece.CELLS )[ 0 ] + + fig, axes = plt.subplots( 1, 2, figsize=( 16, 8 ) ) + + # Plot 1: τ vs σ_n + ax1 = axes[ 0 ] + ax1.scatter( sigmaN, tau, c=depth, cmap='turbo_r', s=20, alpha=0.8 ) + sigmaRange = np.linspace( 0, np.max( sigmaN ), 100 ) + tauCritical = cohesion + mu * sigmaRange + ax1.plot( sigmaRange, tauCritical, 'k--', linewidth=2, label=f'M-C (C={cohesion} bar, φ={phi}°)' ) + ax1.set_xlabel( 'Normal Stress [bar]' ) + ax1.set_ylabel( 'Shear Stress [bar]' ) + ax1.legend() + ax1.grid( True, alpha=0.3 ) + ax1.set_title( 'Mohr-Coulomb Diagram' ) + + # Plot 2: SCU vs σ_n + ax2 = axes[ 1 ] + sc2 = ax2.scatter( sigmaN, SCU, c=depth, cmap='turbo_r', s=20, alpha=0.8 ) + ax2.axhline( y=1.0, color='r', linestyle='--', label='Failure (SCU=1)' ) + ax2.set_xlabel( 'Normal Stress [bar]' ) + ax2.set_ylabel( 'SCU [-]' ) + ax2.legend() + ax2.grid( True, alpha=0.3 ) + ax2.set_title( 'Shear Capacity Utilization' ) + ax2.set_ylim( bottom=0 ) + + plt.colorbar( sc2, ax=ax2, label='Depth [m]' ) + plt.tight_layout() + + if save: + years = time / ( 365.25 * 24 * 3600 ) + filename = f'mohr_coulomb_phi{phi}_c{cohesion}_{years:.0f}y.png' + plt.savefig( path / filename, dpi=300, bbox_inches='tight' ) + print( f" 📊 Plot saved: {filename}" ) + + if show: + plt.show() + else: + plt.close() + + # ------------------------------------------------------------------- + @staticmethod + def loadReferenceData( time: float, + scriptDir: str | Path | None = None, + profileId: int = 1 ) -> dict[ str, npt.NDArray | None ]: + """Load GEOS and analytical reference data for comparison. + + Parameters + ---------- + time : float + Current simulation time in seconds + scriptDir : str or Path, optional + Directory containing reference data files. If None, uses current directory. + profileId : int, optional + Profile ID to extract from Excel (default: 1) + + Returns: + ------- + dict + Dictionary with keys 'geos' and 'analytical', each containing numpy arrays or None + Format: {'geos': array or None, 'analytical': array or None} + + For GEOS data from Excel, the array has columns: + [Depth_m, Normal_Stress_bar, Shear_Stress_bar, SCU, X_coordinate_m, Y_coordinate_m] + """ + if scriptDir is None: + scriptDir = os.path.dirname( os.path.abspath( __file__ ) ) + + result: dict[ str, None | npt.NDArray ] = { 'geos': None, 'analytical': None } + + # =================================================================== + # LOAD GEOS DATA - Try Excel first, then CSV + # =================================================================== + + geosFileXLSV = 'geos_data_numerical.xlsx' + geosFileCSV = 'geos_data_numerical.csv' + + # Try Excel format with time-based sheets + geosXLSVPath = os.path.join( scriptDir, geosFileXLSV ) + + if os.path.exists( geosXLSVPath ): + try: + # Generate sheet name based on current time + # Format: t_1.00e+02s + sheetName = f"t_{time:.2e}s" + + print( f" 📂 Loading GEOS data from Excel sheet: '{sheetName}'" ) + + # Try to read the specific sheet + try: + df = pd.read_excel( geosXLSVPath, sheet_name=sheetName ) + + # Filter by ProfileID if column exists + if 'ProfileID' in df.columns: + dfProfile = df[ df[ 'ProfileID' ] == profileId ] + + if len( dfProfile ) == 0: + print( f" ⚠️ ProfileID {profileId} not found in sheet '{sheetName}'" ) + print( f" Available Profile_IDs: {sorted(df['ProfileID'].unique())}" ) + # Take first profile as fallback + availableIds = sorted( df[ 'ProfileID' ].unique() ) + if len( availableIds ) > 0: + fallbackId = availableIds[ 0 ] + print( f" → Using ProfileID {fallbackId} instead" ) + dfProfile = df[ df[ 'ProfileID' ] == fallbackId ] + else: + print( f" ✅ Loaded ProfileID {profileId}: {len(dfProfile)} points" ) + + # Extract relevant columns in the expected order + # Expected: [Depth, Normal_Stress, Shear_Stress, SCU, ...] + columnsToExtract = [ 'Depth_m', 'Normal_Stress_bar', 'Shear_Stress_bar', 'SCU' ] + + # Check which columns exist + availableColumns = [ col for col in columnsToExtract if col in dfProfile.columns ] + + if len( availableColumns ) > 0: + result[ 'geos' ] = dfProfile[ availableColumns ].values + print( f" Extracted columns: {availableColumns}" ) + else: + print( " ⚠️ No expected columns found in DataFrame" ) + print( f" Available columns: {list(dfProfile.columns)}" ) + else: + # No ProfileID column, use all data + print( " ℹ️ No ProfileID column, using all data" ) + columnsToExtract = [ 'Depth_m', 'Normal_Stress_bar', 'Shear_Stress_bar', 'SCU' ] + availableColumns = [ col for col in columnsToExtract if col in df.columns ] + + if len( availableColumns ) > 0: + result[ 'geos' ] = df[ availableColumns ].values + print( f" ✅ Loaded {len(result['geos'])} points" ) + + except ValueError: + # Sheet not found, try to find closest time + print( f" ⚠️ Sheet '{sheetName}' not found, searching for closest time..." ) + + # Read all sheet names + xlFile = pd.ExcelFile( geosXLSVPath ) + sheetNames = xlFile.sheetNames + + # Extract times from sheet names + sheetTimes = [] + for sname in sheetNames: + if sname.startswith( 't_' ) and sname.endswith( 's' ): + try: + # Extract time: t_1.00e+02s -> 100.0 + timeStr = sname[ 2:-1 ] # Remove 't_' and 's' + sheetTime = float( timeStr ) + sheetTimes.append( ( sheetTime, sname ) ) + except Exception: + continue + + if sheetTimes: + # Find closest time + sheetTimes.sort( key=lambda x: abs( x[ 0 ] - time ) ) + closestTime, closestSheet = sheetTimes[ 0 ] + timeDiff = abs( closestTime - time ) + + print( f" → Using closest sheet: '{closestSheet}' (Δt={timeDiff:.2e}s)" ) + df = pd.read_excel( geosXLSVPath, sheet_name=closestSheet ) + + # Filter by ProfileID + if 'ProfileID' in df.columns: + dfProfile = df[ df[ 'ProfileID' ] == profileId ] + + if len( dfProfile ) == 0: + # Fallback to first profile + availableIds = sorted( df[ 'ProfileID' ].unique() ) + if len( availableIds ) > 0: + dfProfile = df[ df[ 'ProfileID' ] == availableIds[ 0 ] ] + print( f" → Using ProfileID {availableIds[0]}" ) + + columnsToExtract = [ 'Depth_m', 'Normal_Stress_bar', 'Shear_Stress_bar', + 'SCU' ] # TODO check + availableColumns = [ col for col in columnsToExtract if col in dfProfile.columns ] + + if len( availableColumns ) > 0: + result[ 'geos' ] = dfProfile[ availableColumns ].values + print( f" ✅ Loaded {len(result['geos'])} points" ) + else: + # Use all data + columnsToExtract = [ 'Depth_m', 'Normal_Stress_bar', 'Shear_Stress_bar', 'SCU' ] + availableColumns = [ col for col in columnsToExtract if col in df.columns ] + + if len( availableColumns ) > 0: + result[ 'geos' ] = df[ availableColumns ].values + print( f" ✅ Loaded {len(result['geos'])} points" ) + else: + print( " ⚠️ No valid time sheets found in Excel file" ) + + except ImportError: + print( " ⚠️ pandas not available, cannot read Excel file" ) + except Exception as e: + print( f" ⚠️ Error reading Excel: {e}" ) + import traceback + traceback.print_exc() + + # Fallback to CSV if Excel not found or failed + if result[ 'geos' ] is None: + geosCSVPath = os.path.join( scriptDir, geosFileCSV ) + if os.path.exists( geosCSVPath ): + try: + result[ 'geos' ] = np.loadtxt( geosCSVPath, delimiter=',', skiprows=1 ) + print( f" ✅ GEOS data loaded from CSV: {len(result['geos'])} points" ) + except Exception as e: + print( f" ⚠️ Error reading CSV: {e}" ) + + # =================================================================== + # LOAD ANALYTICAL DATA + # =================================================================== + + analyticalFile = 'analyticalData.csv' + analyticalPath = os.path.join( scriptDir, analyticalFile ) + + if os.path.exists( analyticalPath ): + try: + result[ 'analytical' ] = np.loadtxt( analyticalPath, delimiter=',', skiprows=1 ) + print( f" ✅ Analytical data loaded: {len(result['analytical'])} points" ) + except Exception as e: + print( f" ⚠️ Error loading analytical data: {e}" ) + + return result + + # ------------------------------------------------------------------- + def plotDepthProfiles( self, + surface: pv.PolyData, + time: float, + path: Path, + show: bool = True, + save: bool = True, + profileStartPoints: list[ tuple[ float, float ] ] | None = None, + maxProfilePoints: int = 1000, + referenceProfileId: int = 1, + # profileSearchRadius: float | None = None, + showProfileExtractor: bool = True, + ) -> None: + """Plot vertical profiles along the fault showing stress and SCU vs depth.""" + print( " 📊 Creating depth profiles " ) + + # Extract data + centers = surface.cell_data[ 'elementCenter' ] + depth = centers[ :, 2 ] + sigmaN = surface.cell_data[ 'sigmaNEffective' ] + tau = surface.cell_data[ 'tauEffective' ] + SCU = surface.cell_data[ 'SCU' ] + SCU = np.sqrt( SCU**2 ) + surface.cell_data[ 'deltaSCU' ] + + # Extraire les IDs de faille + faultIds = None + if 'FaultMask' in surface.cell_data: + faultIds = surface.cell_data[ 'FaultMask' ] + print( f" 📋 Detected {len(np.unique(faultIds[faultIds > 0]))} distinct faults" ) + elif 'attribute' in surface.cell_data: + faultIds = surface.cell_data[ 'attribute' ] + print( " 📋 Using 'attribute' field for fault identification" ) + else: + print( " ⚠️ No fault IDs found - profiles may jump between faults" ) + + # =================================================================== + # LOAD REFERENCE DATA (GEOS + Analytical) + # =================================================================== + scriptDir = os.path.dirname( os.path.abspath( __file__ ) ) + referenceData = Visualizer.loadReferenceData( time, scriptDir, profileId=referenceProfileId ) + + geosData = referenceData[ 'geos' ] + analyticalData = referenceData[ 'analytical' ] + + # =================================================================== + # PROFILE EXTRACTION SETUP + # =================================================================== + + # Get fault bounds + xMin, xMax = np.min( centers[ :, 0 ] ), np.max( centers[ :, 0 ] ) + yMin, yMax = np.min( centers[ :, 1 ] ), np.max( centers[ :, 1 ] ) + zMin, zMax = np.min( depth ), np.max( depth ) + + # Auto-compute search radius if not provided + xRange = xMax - xMin + yRange = yMax - yMin + zMax - zMin + + if self.profileSearchRadius is not None: + searchRadius = self.profileSearchRadius + else: + searchRadius = min( xRange, yRange ) * 0.15 + + # Auto-generate profile points if not provided + if profileStartPoints is None: + print( " ⚠️ No profileStartPoints provided, auto-generating 5 profiles..." ) + nProfiles = 5 + + # Determine dominant fault direction + if xRange > yRange: + coordName = 'X' + fixedValue = ( yMin + yMax ) / 2 + samplePositions = np.linspace( xMin, xMax, nProfiles ) + profileStartPoints = [ ( x, fixedValue ) for x in samplePositions ] + else: + coordName = 'Y' + fixedValue = ( xMin + xMax ) / 2 + samplePositions = np.linspace( yMin, yMax, nProfiles ) + profileStartPoints = [ ( fixedValue, y ) for y in samplePositions ] + + print( f" Auto-generated {nProfiles} profiles along {coordName} direction" ) + + nProfiles = len( profileStartPoints ) + + # =================================================================== + # CREATE FIGURE + # =================================================================== + + fig, axes = plt.subplots( 1, 4, figsize=( 24, 12 ) ) + colors = plt.cm.RdYlGn( np.linspace( 0, 1, nProfiles ) ) # type: ignore [attr-defined] + + print( f" 📍 Processing {nProfiles} profiles:" ) + print( f" Depth range: [{zMin:.1f}, {zMax:.1f}]m" ) + + successfulProfiles = 0 + + # =================================================================== + # EXTRACT AND PLOT PROFILES + # =================================================================== + + for i, ( xPos, yPos, zPos ) in enumerate( profileStartPoints ): + print( f" → Profile {i+1}: starting at ({xPos:.1f}, {yPos:.1f}, {zPos:.1f})" ) + + depthsSigma, profileSigmaN, PathXSigma, PathYSigma = ProfileExtractor.extractAdaptiveProfile( + centers, sigmaN, xPos, yPos, searchRadius ) + + depthsTau, profileTau, _, _ = ProfileExtractor.extractAdaptiveProfile( centers, tau, xPos, yPos, + searchRadius ) + + depthsSCU, profileSCU, _, _ = ProfileExtractor.extractAdaptiveProfile( centers, SCU, xPos, yPos, + searchRadius ) + + depthsDeltaSCU, profileDeltaSCU, _, _ = ProfileExtractor.extractAdaptiveProfile( + centers, SCU, xPos, yPos, searchRadius ) + + # Calculate path length + if len( PathXSigma ) > 1: + pathLength = np.sum( + np.sqrt( np.diff( PathXSigma )**2 + np.diff( PathYSigma )**2 + np.diff( depthsSigma )**2 ) ) + print( + f" Path length: {pathLength:.1f}m (horizontal displacement: {np.abs(PathXSigma[-1] - PathXSigma[0]):.1f}m)" + ) + + if showProfileExtractor: + ProfileExtractor.plotProfilePath3D( surface=surface, + pathX=PathXSigma, + pathY=PathYSigma, + pathZ=depthsSigma, + profileValues=profileSigmaN, + scalarName='SCU', + savePath=path, + show=show ) + + # Check if we have enough points + minPoints = 3 + nPoints = len( depthsSigma ) + + if nPoints >= minPoints: + label = f'Profile {i+1} → ({xPos:.0f}, {yPos:.0f})' + + # Plot 1: Normal stress vs depth + axes[ 0 ].plot( profileSigmaN, + depthsSigma, + color=colors[ i ], + label=label, + linewidth=2.5, + alpha=0.8, + marker='o', + markersize=3, + markevery=2 ) + + # Plot 2: Shear stress vs depth + axes[ 1 ].plot( profileTau, + depthsTau, + color=colors[ i ], + label=label, + linewidth=2.5, + alpha=0.8, + marker='o', + markersize=3, + markevery=2 ) + + # Plot 3: SCU vs depth + axes[ 2 ].plot( profileSCU, + depthsSCU, + color=colors[ i ], + label=label, + linewidth=2.5, + alpha=0.8, + marker='o', + markersize=3, + markevery=2 ) + + # Plot 4: Detla SCU vs depth + axes[ 3 ].plot( profileDeltaSCU, + depthsDeltaSCU, + color=colors[ i ], + label=label, + linewidth=2.5, + alpha=0.8, + marker='o', + markersize=3, + markevery=2 ) + + successfulProfiles += 1 + print( f" ✅ {nPoints} points found" ) + else: + print( f" ⚠️ Insufficient points ({nPoints}), skipping" ) + + if successfulProfiles == 0: + print( " ❌ No valid profiles found!" ) + plt.close() + return + + # =================================================================== + # ADD REFERENCE DATA (GEOS + Analytical) - Only once + # =================================================================== + + if geosData is not None: + # Colonnes: [Depth_m, Normal_Stress_bar, Shear_Stress_bar, SCU] + # Index: [0, 1, 2, 3] + + axes[ 0 ].plot( geosData[ :, 1 ] * 10, + geosData[ :, 0 ], + 'o', + color='blue', + markersize=6, + label='GEOS Contact Solver', + alpha=0.7, + mec='k', + mew=1, + fillstyle='none' ) + + axes[ 1 ].plot( geosData[ :, 2 ] * 10, + geosData[ :, 0 ], + 'o', + color='blue', + markersize=6, + label='GEOS Contact Solver', + alpha=0.7, + mec='k', + mew=1, + fillstyle='none' ) + + if geosData.shape[ 1 ] > 3: # SCU column exists + axes[ 2 ].plot( geosData[ :, 3 ], + geosData[ :, 0 ], + 'o', + color='blue', + markersize=6, + label='GEOS Contact Solver', + alpha=0.7, + mec='k', + mew=1, + fillstyle='none' ) + + if analyticalData is not None: + # Format analytique (peut varier) + axes[ 0 ].plot( analyticalData[ :, 1 ] * 10, + analyticalData[ :, 0 ], + '--', + color='darkorange', + linewidth=2, + label='Analytical', + alpha=0.8 ) + if analyticalData.shape[ 1 ] > 2: + axes[ 1 ].plot( analyticalData[ :, 2 ] * 10, + analyticalData[ :, 0 ], + '--', + color='darkorange', + linewidth=2, + label='Analytical', + alpha=0.8 ) + + # =================================================================== + # CONFIGURE PLOTS + # =================================================================== + + fsize = 14 + + # Plot 1: Normal Stress + axes[ 0 ].set_xlabel( 'Normal Stress σₙ [bar]', fontsize=fsize, weight="bold" ) + axes[ 0 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 0 ].set_title( 'Normal Stress Profile', fontsize=fsize + 2, weight="bold" ) + axes[ 0 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 0 ].legend( loc='upper left', fontsize=fsize - 2 ) + axes[ 0 ].tick_params( labelsize=fsize - 2 ) + + # Plot 2: Shear Stress + axes[ 1 ].set_xlabel( 'Shear Stress τ [bar]', fontsize=fsize, weight="bold" ) + axes[ 1 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 1 ].set_title( 'Shear Stress Profile', fontsize=fsize + 2, weight="bold" ) + axes[ 1 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 1 ].legend( loc='upper left', fontsize=fsize - 2 ) + axes[ 1 ].tick_params( labelsize=fsize - 2 ) + + # Plot 3: SCU + axes[ 2 ].set_xlabel( 'SCU [-]', fontsize=fsize, weight="bold" ) + axes[ 2 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 2 ].set_title( 'Shear Capacity Utilization', fontsize=fsize + 2, weight="bold" ) + axes[ 2 ].axvline( x=0.8, color='forestgreen', linestyle='--', linewidth=2, label='Critical (0.8)' ) + axes[ 2 ].axvline( x=1.0, color='red', linestyle='--', linewidth=2, label='Failure (1.0)' ) + axes[ 2 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 2 ].legend( loc='upper right', fontsize=fsize - 2 ) + axes[ 2 ].tick_params( labelsize=fsize - 2 ) + axes[ 2 ].set_xlim( left=0 ) + + # Plot 4: Delta SCU + axes[ 3 ].set_xlabel( 'Δ SCU [-]', fontsize=fsize, weight="bold" ) + axes[ 3 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 3 ].set_title( 'Delta SCU', fontsize=fsize + 2, weight="bold" ) + axes[ 3 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 3 ].legend( loc='upper right', fontsize=fsize - 2 ) + axes[ 3 ].tick_params( labelsize=fsize - 2 ) + axes[ 3 ].set_xlim( left=0, right=2 ) + + # Change verticale scale + if self.maxDepthProfiles is not None: + for i in range( len( axes ) ): + axes[ i ].set_ylim( bottom=self.maxDepthProfiles ) + + if self.minDepthProfiles is not None: + for i in range( len( axes ) ): + axes[ i ].set_ylim( top=self.minDepthProfiles ) + + # Overall title + years = time / ( 365.25 * 24 * 3600 ) + fig.suptitle( f'Fault Depth Profiles - t={years:.1f} years', fontsize=fsize + 2, fontweight='bold', y=0.98 ) + + plt.tight_layout( rect=( 0, 0, 1, 0.96 ) ) + + # Save + if save: + filename = f'depth_profiles_{years:.0f}y.png' + plt.savefig( path / filename, dpi=300, bbox_inches='tight' ) + print( f" 💾 Depth profiles saved: {filename}" ) + + # Show + if show: + plt.show() + else: + plt.close() + + # ------------------------------------------------------------------- + def plotVolumeStressProfiles( self: Self, + volumeMesh: pv.DataSet, + faultSurface: pv.PolyData, + time: float, + path: Path, + show: bool = True, + save: bool = True, + profileStartPoints: list[ tuple[ float, float, float ] ] | None = None, + maxProfilePoints: int = 1000, + # profileSearchRadius: float | None = None, + # maxDepthProfile: float | None = None, + # minDepthProfile: float | None = None + ) -> None: + """Plot stress profiles in volume cells adjacent to the fault. + + Extracts profiles through contributing cells on BOTH sides of the fault + Shows plus side and minus side on the same plots for comparison. + + NOTE: Cette fonction utilise extractAdaptiveProfile pour les VOLUMES + car volumeMesh n'est PAS un maillage surfacique. + La méthode topologique (extractVerticalProfileTopologyBased) + est réservée aux maillages SURFACIQUES (faultSurface). + """ + print( " 📊 Creating volume stress profiles (both sides)" ) + + # =================================================================== + # CHECK IF REQUIRED DATA EXISTS + # =================================================================== + + requiredFields = [ 'sigma1', 'sigma2', 'sigma3', 'side', 'elementCenter' ] + + for field in requiredFields: + if isAttributeInObject( volumeMesh, field, Piece.CELLS ): + print( f" ⚠️ Missing required field: {field}" ) + return + + # Check for pressure + if isAttributeInObject( volumeMesh, 'pressure_bar', Piece.CELLS): + pressureField = 'pressure_bar' + pressure = getArrayInObject( volumeMesh, pressureField, Piece.CELLS ) + elif isAttributeInObject( volumeMesh, 'pressure', Piece.CELLS ): + pressureField = 'pressure' + pressure = getArrayInObject( volumeMesh, pressureField, Piece.CELLS ) / 1e5 + print( " ℹ️ Converting pressure from Pa to bar" ) + else: + print( " ⚠️ No pressure field found" ) + pressure = None + + # Extract volume data + centers = getArrayInObject( volumeMesh, 'elementCenter', Piece.CELLS ) + sigma1 = getArrayInObject( volumeMesh, 'sigma1', Piece.CELLS ) + sigma2 = getArrayInObject( volumeMesh, 'sigma2', Piece.CELLS ) + sigma3 = getArrayInObject( volumeMesh, 'sigma3', Piece.CELLS ) + sideData = getArrayInObject( volumeMesh, 'side', Piece.CELLS ) + + # =================================================================== + # FILTER CELLS BY SIDE (BOTH PLUS AND MINUS) + # =================================================================== + + # Plus side (side = 1 or 3) + maskPlus = ( sideData == 1 ) | ( sideData == 3 ) + centersPlus = centers[ maskPlus ] + sigma1Plus = sigma1[ maskPlus ] + sigma2Plus = sigma2[ maskPlus ] + sigma3Plus = sigma3[ maskPlus ] + if pressure is not None: + pressurePlus = pressure[ maskPlus ] + + # Minus side (side = 2 or 3) + maskMinus = ( sideData == 2 ) | ( sideData == 3 ) + centersMinus = centers[ maskMinus ] + sigma1Minus = sigma1[ maskMinus ] + sigma2Minus = sigma2[ maskMinus ] + sigma3Minus = sigma3[ maskMinus ] + if pressure is not None: + pressureMinus = pressure[ maskMinus ] + + # Créer subset de cellData pour le côté plus + cellDataPlus = {} + cellDataMinus = {} + for key in volumeMesh.GetCellData().items(): + cellDataPlus[ key ] = getArrayInObject( volumeMesh, key )[ maskPlus ] + cellDataMinus[ key ] = getArrayInObject( volumeMesh, key )[ maskMinus ] + + + print( f" 📍 Plus side: {len(centersPlus):,} cells" ) + print( f" 📍 Minus side: {len(centersMinus):,} cells" ) + + if len( centersPlus ) == 0 and len( centersMinus ) == 0: + print( " ⚠️ No contributing cells found!" ) + return + + # =================================================================== + # GET FAULT BOUNDS + # =================================================================== + + faultCenters = faultSurface.cell_data[ 'elementCenter' ] + + xMin, xMax = np.min( faultCenters[ :, 0 ] ), np.max( faultCenters[ :, 0 ] ) + yMin, yMax = np.min( faultCenters[ :, 1 ] ), np.max( faultCenters[ :, 1 ] ) + zMin, zMax = np.min( faultCenters[ :, 2 ] ), np.max( faultCenters[ :, 2 ] ) + + xRange = xMax - xMin + yRange = yMax - yMin + zMax - zMin + + # Search radius (pour extractAdaptiveProfile sur volumes) + if self.profileSearchRadius is not None: + searchRadius = self.profileSearchRadius + else: + searchRadius = min( xRange, yRange ) * 0.2 + + # =================================================================== + # AUTO-GENERATE PROFILE POINTS IF NOT PROVIDED + # =================================================================== + + if profileStartPoints is None: + print( " ⚠️ No profileStartPoints provided, auto-generating..." ) + nProfiles = 3 + + if xRange > yRange: + coordName = 'X' + fixedValue = ( yMin + yMax ) / 2 + samplePositions = np.linspace( xMin, xMax, nProfiles ) + profileStartPoints = [ ( x, fixedValue, zMax ) for x in samplePositions ] + else: + coordName = 'Y' + fixedValue = ( xMin + xMax ) / 2 + samplePositions = np.linspace( yMin, yMax, nProfiles ) + profileStartPoints = [ ( fixedValue, y, zMax ) for y in samplePositions ] + + print( f" Auto-generated {nProfiles} profiles along {coordName}" ) + + nProfiles = len( profileStartPoints ) + + # =================================================================== + # CREATE FIGURE WITH 5 SUBPLOTS + # =================================================================== + + fig, axes = plt.subplots( 1, 5, figsize=( 22, 10 ) ) + + # Colors: different for plus and minus sides + colorsPlus = plt.cm.Reds( np.linspace( 0.4, 0.9, nProfiles ) ) # type: ignore [attr-defined] + colorsMinus = plt.cm.Blues( np.linspace( 0.4, 0.9, nProfiles ) ) # type: ignore [attr-defined] + + print( f" 📍 Processing {nProfiles} volume profiles:" ) + print( f" Depth range: [{zMin:.1f}, {zMax:.1f}]m" ) + + successfulProfiles = 0 + + # =================================================================== + # EXTRACT AND PLOT PROFILES FOR BOTH SIDES + # =================================================================== + + for i, ( xPos, yPos, zPos ) in enumerate( profileStartPoints ): + print( f"\n → Profile {i+1}: starting at ({xPos:.1f}, {yPos:.1f}, {zPos:.1f})" ) + + # ================================================================ + # PLUS SIDE + # ================================================================ + if len( centersPlus ) > 0: + print( " Processing PLUS side..." ) + + # Pour VOLUMES, utiliser extractAdaptiveProfile avec cellData + depthsSigma1Plus, profileSigma1Plus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersPlus, sigma1Plus, xPos, yPos, zPos, searchRadius, verbose=True, cellData=cellDataPlus ) + + depthsSigma2Plus, profileSigma2Plus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersPlus, sigma2Plus, xPos, yPos, zPos, searchRadius, verbose=False, cellData=cellDataPlus ) + + depthsSigma3Plus, profileSigma3Plus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersPlus, sigma3Plus, xPos, yPos, zPos, searchRadius, verbose=False, cellData=cellDataPlus ) + + if pressure is not None: + depthsPressurePlus, profilePressurePlus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersPlus, + pressurePlus, + xPos, + yPos, + zPos, + searchRadius, + verbose=False, + cellData=cellDataPlus ) + + if len( depthsSigma1Plus ) >= 3: + labelPlus = 'Plus side' + + # Plot Pressure + if pressure is not None: + axes[ 0 ].plot( profilePressurePlus, + depthsPressurePlus, + color=colorsPlus[ i ], + label=labelPlus if i == 0 else '', + linewidth=2.5, + alpha=0.8, + marker='o', + markersize=3, + markevery=2 ) + + # Plot σ1 + axes[ 1 ].plot( profileSigma1Plus, + depthsSigma1Plus, + color=colorsPlus[ i ], + label=labelPlus if i == 0 else '', + linewidth=2.5, + alpha=0.8, + marker='o', + markersize=3, + markevery=2 ) + + # Plot σ2 + axes[ 2 ].plot( profileSigma2Plus, + depthsSigma2Plus, + color=colorsPlus[ i ], + label=labelPlus if i == 0 else '', + linewidth=2.5, + alpha=0.8, + marker='o', + markersize=3, + markevery=2 ) + + # Plot σ3 + axes[ 3 ].plot( profileSigma3Plus, + depthsSigma3Plus, + color=colorsPlus[ i ], + label=labelPlus if i == 0 else '', + linewidth=2.5, + alpha=0.8, + marker='o', + markersize=3, + markevery=2 ) + + # Plot All stresses + axes[ 4 ].plot( profileSigma1Plus, + depthsSigma1Plus, + color=colorsPlus[ i ], + linewidth=2.5, + alpha=0.8, + linestyle='-', + marker="o", + markersize=2, + markevery=2 ) + axes[ 4 ].plot( profileSigma2Plus, + depthsSigma2Plus, + color=colorsPlus[ i ], + linewidth=2.0, + alpha=0.6, + linestyle='-', + marker="s", + markersize=2, + markevery=2 ) + axes[ 4 ].plot( profileSigma3Plus, + depthsSigma3Plus, + color=colorsPlus[ i ], + linewidth=2.5, + alpha=0.8, + linestyle='-', + marker="v", + markersize=2, + markevery=2 ) + + print( f" ✅ PLUS: {len(depthsSigma1Plus)} points" ) + successfulProfiles += 1 + + # ================================================================ + # MINUS SIDE + # ================================================================ + if len( centersMinus ) > 0: + print( " Processing MINUS side..." ) + + # Pour VOLUMES, utiliser extractAdaptiveProfile avec cellData + depthsSigma1Minus, profileSigma1Minus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersMinus, sigma1Minus, xPos, yPos, zPos, searchRadius, verbose=True, cellData=cellDataMinus ) + + depthsSigma2Minus, profileSigma2Minus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersMinus, sigma2Minus, xPos, yPos, zPos, searchRadius, verbose=False, cellData=cellDataMinus ) + + depthsSigma3Minus, profileSigma3Minus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersMinus, sigma3Minus, xPos, yPos, zPos, searchRadius, verbose=False, cellData=cellDataMinus ) + + if pressure is not None: + depthsPressureMinus, profilePressureMinus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersMinus, + pressureMinus, + xPos, + yPos, + zPos, + searchRadius, + verbose=False, + cellData=cellDataMinus ) + + if len( depthsSigma1Minus ) >= 3: + labelMinus = 'Minus side' + + # Plot Pressure + if pressure is not None: + axes[ 0 ].plot( profilePressureMinus, + depthsPressureMinus, + color=colorsMinus[ i ], + label=labelMinus if i == 0 else '', + linewidth=2.5, + alpha=0.8, + marker='s', + markersize=3, + markevery=2 ) + + # Plot σ1 + axes[ 1 ].plot( profileSigma1Minus, + depthsSigma1Minus, + color=colorsMinus[ i ], + label=labelMinus if i == 0 else '', + linewidth=2.5, + alpha=0.8, + marker='s', + markersize=3, + markevery=2 ) + + # Plot σ2 + axes[ 2 ].plot( profileSigma2Minus, + depthsSigma2Minus, + color=colorsMinus[ i ], + label=labelMinus if i == 0 else '', + linewidth=2.5, + alpha=0.8, + marker='s', + markersize=3, + markevery=2 ) + + # Plot σ3 + axes[ 3 ].plot( profileSigma3Minus, + depthsSigma3Minus, + color=colorsMinus[ i ], + label=labelMinus if i == 0 else '', + linewidth=2.5, + alpha=0.8, + marker='s', + markersize=3, + markevery=2 ) + + # Plot All stresses + axes[ 4 ].plot( profileSigma1Minus, + depthsSigma1Minus, + color=colorsMinus[ i ], + linewidth=2.5, + alpha=0.8, + linestyle='-', + marker="o", + markersize=2, + markevery=2 ) + axes[ 4 ].plot( profileSigma2Minus, + depthsSigma2Minus, + color=colorsMinus[ i ], + linewidth=2.0, + alpha=0.6, + linestyle='-', + marker="s", + markersize=2, + markevery=2 ) + axes[ 4 ].plot( profileSigma3Minus, + depthsSigma3Minus, + color=colorsMinus[ i ], + linewidth=2.5, + alpha=0.8, + linestyle='-', + marker='v', + markersize=2, + markevery=2 ) + + print( f" ✅ MINUS: {len(depthsSigma1Minus)} points" ) + successfulProfiles += 1 + + if successfulProfiles == 0: + print( " ❌ No valid profiles found!" ) + plt.close() + return + + # =================================================================== + # CONFIGURE PLOTS + # =================================================================== + + fsize = 14 + + # Plot 0: Pressure + axes[ 0 ].set_xlabel( 'Pressure [bar]', fontsize=fsize, weight="bold" ) + axes[ 0 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 0 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 0 ].legend( loc='best', fontsize=fsize - 2 ) + axes[ 0 ].tick_params( labelsize=fsize - 2 ) + + if pressure is None: + axes[ 0 ].text( 0.5, + 0.5, + 'No pressure data available', + ha='center', + va='center', + transform=axes[ 0 ].transAxes, + fontsize=fsize, + style='italic', + color='gray' ) + + # Plot 1: σ1 (Maximum principal stress) + axes[ 1 ].set_xlabel( 'σ₁ (Max Principal) [bar]', fontsize=fsize, weight="bold" ) + axes[ 1 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 1 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 1 ].legend( loc='best', fontsize=fsize - 2 ) + axes[ 1 ].tick_params( labelsize=fsize - 2 ) + + # Plot 2: σ2 (Intermediate principal stress) + axes[ 2 ].set_xlabel( 'σ₂ (Inter Principal) [bar]', fontsize=fsize, weight="bold" ) + axes[ 2 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 2 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 2 ].legend( loc='best', fontsize=fsize - 2 ) + axes[ 2 ].tick_params( labelsize=fsize - 2 ) + + # Plot 3: σ3 (Min principal stress) + axes[ 3 ].set_xlabel( 'σ₃ (Min Principal) [bar]', fontsize=fsize, weight="bold" ) + axes[ 3 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 3 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 3 ].legend( loc='best', fontsize=fsize - 2 ) + axes[ 3 ].tick_params( labelsize=fsize - 2 ) + + # Plot 4: All stresses together + axes[ 4 ].set_xlabel( 'Principal Stresses [bar]', fontsize=fsize, weight="bold" ) + axes[ 4 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 4 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 4 ].tick_params( labelsize=fsize - 2 ) + + # Add legend for line styles + customLines = [ + Line2D( [ 0 ], [ 0 ], color='red', linewidth=2.5, marker=None, label='Plus side', alpha=0.5 ), + Line2D( [ 0 ], [ 0 ], color='blue', linewidth=2.5, marker=None, label='Minus side', alpha=0.5 ), + Line2D( [ 0 ], [ 0 ], color='gray', linewidth=2.5, linestyle='-', marker='o', label='σ₁ (max)' ), + Line2D( [ 0 ], [ 0 ], color='gray', linewidth=2.0, linestyle='-', marker='s', label='σ₂ (inter)' ), + Line2D( [ 0 ], [ 0 ], color='gray', linewidth=2.5, linestyle='-', marker='v', label='σ₃ (min)' ) + ] + axes[ 4 ].legend( handles=customLines, loc='best', fontsize=fsize - 3, ncol=1 ) + + # Change verticale scale + if self.maxDepthProfile is not None: + for i in range( len( axes ) ): + axes[ i ].set_ylim( bottom=self.maxDepthProfiles ) + + if self.minDepthProfiles is not None: + for i in range( len( axes ) ): + axes[ i ].set_ylim( top=self.minDepthProfiles ) + + # Overall title + years = time / ( 365.25 * 24 * 3600 ) + fig.suptitle( f'Volume Stress Profiles - Both Sides Comparison - t={years:.1f} years', + fontsize=fsize + 2, + fontweight='bold', + y=0.98 ) + + plt.tight_layout( rect=( 0, 0, 1, 0.96 ) ) + + # Save + if save: + filename = f'volume_stress_profiles_both_sides_{years:.0f}y.png' + plt.savefig( path / filename, dpi=300, bbox_inches='tight' ) + print( f" 💾 Volume profiles saved: {filename}" ) + + # Show + if show: + plt.show() + else: + plt.close() + + # ------------------------------------------------------------------- + def plotAnalyticalVsNumericalComparison( self: Self, + volumeMesh: pv.PolyData, + faultSurface: pv.PolyData, + time: float, + path: Path, + show: bool = True, + save: bool = True, + profileStartPoints: list[ tuple[ int, int, int ] ] | None = None, + referenceProfileId: int = 1, + # profileSearchRadius: float| None = None, + # minDepthProfile: float | None = None, + # maxDepthProfile: float | None = None, + ) -> None: + """Plot comparison between analytical fault stresses (Anderson formulas) and numerical tensor projection - COMBINED PLOTS ONLY. + + Parameters + ---------- + volumeMesh : pyvista.UnstructuredGrid + Volume mesh with principal stresses AND analytical stresses + faultSurface : pyvista.PolyData + Fault surface mesh with projected stresses + time : float + Simulation time + path : Path + Output directory + show : bool + Show plot interactively + save : bool + Save plot to file + profileStartPoints : list of tuples + Starting points (x, y, z) for profiles + referenceProfileId : int + Which profile ID to load from Excel reference data + """ + print( "\n 📊 Creating Analytical vs Numerical Comparison" ) + + # =================================================================== + # CHECK IF ANALYTICAL DATA EXISTS + # =================================================================== + + requiredAnalytical = [ 'sigmaNAnalytical', 'tauAnalytical', 'side', 'elementCenter' ] + + for field in requiredAnalytical: + if field not in volumeMesh.cell_data: + print( f" ⚠️ Missing analytical field: {field}" ) + print( " Analytical stresses not computed in volume mesh" ) + return + + # Check numerical data on fault surface + if 'sigmaNEffective' not in faultSurface.cell_data: + print( " ⚠️ Missing numerical stress data on fault surface" ) + return + + # =================================================================== + # LOAD REFERENCE DATA (GEOS Contact Solver) + # =================================================================== + + print( " 📂 Loading GEOS Contact Solver reference data..." ) + scriptDir = os.path.dirname( os.path.abspath( __file__ ) ) + referenceData = Visualizer.loadReferenceData( time, scriptDir, profileId=referenceProfileId ) + + geosContactData = referenceData.get( 'geos', None ) + + if geosContactData is not None: + print( f" ✅ Loaded {len(geosContactData)} reference points from GEOS Contact Solver" ) + else: + print( " ⚠️ No GEOS Contact Solver reference data found" ) + + # Extraire les IDs de faille + + if 'faultId' in volumeMesh.cell_data: + volumeMesh.cell_data[ 'faultId' ] + + if 'FaultMask' in faultSurface.cell_data: + faultSurface.cell_data[ 'FaultMask' ] + elif 'attribute' in faultSurface.cell_data: + faultSurface.cell_data[ 'attribute' ] + + # =================================================================== + # EXTRACT DATA + # =================================================================== + + # Volume analytical data + centersVolume = volumeMesh.cell_data[ 'elementCenter' ] + sideData = volumeMesh.cell_data[ 'side' ] + sigmaNAnalytical = volumeMesh.cell_data[ 'sigmaNAnalytical' ] + tauAnalytical = volumeMesh.cell_data[ 'tauAnalytical' ] + + # Optional: SCU if available + hasSCUAnalytical = 'SCUAnalytical' in volumeMesh.cell_data + if hasSCUAnalytical: + SCUAnalytical = volumeMesh.cell_data[ 'SCUAnalytical' ] + + # Fault numerical data + centersFault = faultSurface.cell_data[ 'elementCenter' ] + sigmaNNumerical = faultSurface.cell_data[ 'sigmaNEffective' ] + tauNumerical = faultSurface.cell_data[ 'tauEffective' ] + + # Optional: SCU numerical + hasSCUNumerical = 'SCU' in faultSurface.cell_data + if hasSCUNumerical: + SCUNumerical = faultSurface.cell_data[ 'SCU' ] + + # Filter volume by side + maskPlus = ( sideData == 1 ) | ( sideData == 3 ) + maskMinus = ( sideData == 2 ) | ( sideData == 3 ) + + centersPlus = centersVolume[ maskPlus ] + sigmaNAnalyticalPlus = sigmaNAnalytical[ maskPlus ] + tauAnalyticalPlus = tauAnalytical[ maskPlus ] + if hasSCUAnalytical: + SCUAnalyticalPlus = SCUAnalytical[ maskPlus ] + + centersMinus = centersVolume[ maskMinus ] + sigmaNAnalyticalMinus = sigmaNAnalytical[ maskMinus ] + tauAnalyticalMinus = tauAnalytical[ maskMinus ] + if hasSCUAnalytical: + SCUAnalyticalMinus = SCUAnalytical[ maskMinus ] + + print( f" 📍 Plus side: {len(centersPlus):,} cells with analytical data" ) + print( f" 📍 Minus side: {len(centersMinus):,} cells with analytical data" ) + print( f" 📍 Fault surface: {len(centersFault):,} cells with numerical data" ) + + # =================================================================== + # GET FAULT BOUNDS AND PROFILE SETUP + # =================================================================== + + xMin, xMax = np.min( centersFault[ :, 0 ] ), np.max( centersFault[ :, 0 ] ) + yMin, yMax = np.min( centersFault[ :, 1 ] ), np.max( centersFault[ :, 1 ] ) + _zMin, zMax = np.min( centersFault[ :, 2 ] ), np.max( centersFault[ :, 2 ] ) + + xRange = xMax - xMin + yRange = yMax - yMin + + # Search radius + if self.profileSearchRadius is not None: + searchRadius = self.profileSearchRadius + else: + searchRadius = min( xRange, yRange ) * 0.2 + + # Auto-generate profile points if not provided + if profileStartPoints is None: + print( " ⚠️ No profileStartPoints provided, auto-generating..." ) + nProfiles = 3 + + if xRange > yRange: + coordName = 'X' + fixedValue = ( yMin + yMax ) / 2 + samplePositions = np.linspace( xMin, xMax, nProfiles ) + profileStartPoints = [ ( x, fixedValue, zMax ) for x in samplePositions ] + else: + coordName = 'Y' + fixedValue = ( xMin + xMax ) / 2 + samplePositions = np.linspace( yMin, yMax, nProfiles ) + profileStartPoints = [ ( fixedValue, y, zMax ) for y in samplePositions ] + + print( f" Auto-generated {nProfiles} profiles along {coordName}" ) + + nProfiles = len( profileStartPoints ) + + # =================================================================== + # CREATE FIGURE: COMBINED PLOTS ONLY + # 3 columns (σ_n, τ, SCU) x 1 row + # =================================================================== + + fig, axes = plt.subplots( 1, 3, figsize=( 18, 10 ) ) + + print( f" 📍 Processing {nProfiles} profiles for comparison:" ) + + successfulProfiles = 0 + + # =================================================================== + # EXTRACT AND PLOT PROFILES + # =================================================================== + + for i, ( xPos, yPos, zPos ) in enumerate( profileStartPoints ): + print( f"\n → Profile {i+1}: starting at ({xPos:.1f}, {yPos:.1f}, {zPos:.1f})" ) + + # ================================================================ + # PLUS SIDE - ANALYTICAL + # ================================================================ + if len( centersPlus ) > 0: + depthsSnAnaPlus, profileSnAnaPlus, _, _ = ProfileExtractor.extractAdaptiveProfile( centersPlus, + sigmaNAnalyticalPlus, + xPos, + yPos, + zPos, + searchRadius, + verbose=False ) + + depthsTauAnaPlus, profileTauAnaPlus, _, _ = ProfileExtractor.extractAdaptiveProfile( centersPlus, + tauAnalyticalPlus, + xPos, + yPos, + zPos, + searchRadius, + verbose=False ) + + if hasSCUAnalytical: + depthsSCUAnaPlus, profileSCUAnaPlus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersPlus, + SCUAnalyticalPlus, + xPos, + yPos, + zPos, + searchRadius, + verbose=False, + ) + + if len( depthsSnAnaPlus ) >= 3: + # Plot σ_n + axes[ 0 ].plot( profileSnAnaPlus, + depthsSnAnaPlus, + color='red', + linestyle='-', + linewidth=2, + alpha=0.3, + label='Analytical Side +' if i == 0 else '', + marker=None ) + + # Plot τ + axes[ 1 ].plot( profileTauAnaPlus, + depthsTauAnaPlus, + color='red', + linestyle='-', + linewidth=2, + alpha=0.3, + label='Analytical Side +' if i == 0 else '', + marker=None ) + + # Plot SCU if available + if hasSCUAnalytical and len( depthsSCUAnaPlus ) >= 3: + axes[ 2 ].plot( profileSCUAnaPlus, + depthsSCUAnaPlus, + color='red', + linestyle='-', + linewidth=2, + alpha=0.3, + label='Analytical Side +' if i == 0 else '', + marker=None ) + + # ================================================================ + # MINUS SIDE - ANALYTICAL + # ================================================================ + if len( centersMinus ) > 0: + depthsSigmaNAnaMinus, profileSigmaNAnaMinus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersMinus, sigmaNAnalyticalMinus, xPos, yPos, zPos, searchRadius, verbose=False ) + + depthsTauAnaMinus, profileTauAnaMinus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersMinus, tauAnalyticalMinus, xPos, yPos, zPos, searchRadius, verbose=False ) + + if hasSCUAnalytical: + depthsSCUAnaMinus, profileSCUAnaMinus, _, _ = ProfileExtractor.extractAdaptiveProfile( + centersMinus, SCUAnalyticalMinus, xPos, yPos, zPos, searchRadius, verbose=False ) + + if len( depthsSigmaNAnaMinus ) >= 3: + # Plot σ_n + axes[ 0 ].plot( profileSigmaNAnaMinus, + depthsSigmaNAnaMinus, + color='blue', + linestyle='-', + linewidth=2, + alpha=0.3, + label='Analytical Side -' if i == 0 else '', + marker=None ) + + # Plot τ + axes[ 1 ].plot( profileTauAnaMinus, + depthsTauAnaMinus, + color='blue', + linestyle='-', + linewidth=2, + alpha=0.3, + label='Analytical Side -' if i == 0 else '', + marker=None ) + + # Plot SCU if available + if hasSCUAnalytical and len( depthsSCUAnaMinus ) >= 3: + axes[ 2 ].plot( profileSCUAnaMinus, + depthsSCUAnaMinus, + color='blue', + linestyle='-', + linewidth=2, + alpha=0.3, + label='Analytical Side -' if i == 0 else '', + marker=None ) + + # ================================================================ + # AVERAGES - ANALYTICAL (only for first profile to avoid clutter) + # ================================================================ + if i == 0 and len( depthsSigmaNAnaMinus ) >= 3 and len( depthsSnAnaPlus ) >= 3: + # Arithmetic average + avgSigmaNArith = ( profileSigmaNAnaMinus + profileSnAnaPlus ) / 2 + avgTauArith = ( profileTauAnaMinus + profileTauAnaPlus ) / 2 + + axes[ 0 ].plot( avgSigmaNArith, + depthsSigmaNAnaMinus, + color='darkorange', + linestyle='-', + linewidth=2, + alpha=0.6, + label='Arithmetic average' ) + + axes[ 1 ].plot( avgTauArith, + depthsSigmaNAnaMinus, + color='darkorange', + linestyle='-', + linewidth=2, + alpha=0.6, + label='Arithmetic average' ) + + # Geometric average + avgTauGeom = np.sqrt( profileTauAnaMinus * profileTauAnaPlus ) + + axes[ 1 ].plot( avgTauGeom, + depthsSigmaNAnaMinus, + color='purple', + linestyle='-', + linewidth=2, + alpha=0.6, + label='Geometric average' ) + + # Harmonic average + AvgSigmaNHarm = 2 / ( 1 / profileSigmaNAnaMinus + 1 / profileSnAnaPlus ) + AvgTauHarm = 2 / ( 1 / profileTauAnaMinus + 1 / profileTauAnaPlus ) + + axes[ 0 ].plot( AvgSigmaNHarm, + depthsSigmaNAnaMinus, + color='green', + linestyle='-', + linewidth=2, + alpha=0.6, + label='Harmonic average' ) + + axes[ 1 ].plot( AvgTauHarm, + depthsSigmaNAnaMinus, + color='green', + linestyle='-', + linewidth=2, + alpha=0.6, + label='Harmonic average' ) + + # ================================================================ + # NUMERICAL DATA FROM FAULT SURFACE (Continuum) + # ================================================================ + print( " Extracting numerical data from fault surface..." ) + + depthsSigmaNNum, profileSigmaNNum, _, _ = ProfileExtractor.extractAdaptiveProfile( centersFault, + sigmaNNumerical, + xPos, + yPos, + zPos, + searchRadius, + verbose=False ) + + depthsTauNum, profileTauNum, _, _ = ProfileExtractor.extractAdaptiveProfile( centersFault, + tauNumerical, + xPos, + yPos, + zPos, + searchRadius, + verbose=False ) + + if hasSCUNumerical: + depthsSCUNum, profileSCUNum, _, _ = ProfileExtractor.extractAdaptiveProfile( centersFault, + SCUNumerical, + xPos, + yPos, + zPos, + searchRadius, + verbose=False ) + + if len( depthsSigmaNNum ) >= 3: + # Plot numerical with distinct style + axes[ 0 ].plot( profileSigmaNNum, + depthsSigmaNNum, + color='black', + linestyle='-', + linewidth=2, + alpha=0.7, + label='GEOS Continuum' if i == 0 else '', + marker='x', + markersize=5, + markevery=3 ) + + axes[ 1 ].plot( profileTauNum, + depthsTauNum, + color='black', + linestyle='-', + linewidth=2, + alpha=0.7, + label='GEOS Continuum' if i == 0 else '', + marker='x', + markersize=5, + markevery=3 ) + + if hasSCUNumerical and len( depthsSCUNum ) >= 3: + axes[ 2 ].plot( profileSCUNum, + depthsSCUNum, + color='black', + linestyle='-', + linewidth=2, + alpha=0.7, + label='GEOS Continuum' if i == 0 else '', + marker='x', + markersize=5, + markevery=3 ) + + successfulProfiles += 1 + + if successfulProfiles == 0: + print( " ❌ No valid profiles found!" ) + plt.close() + return + + # =================================================================== + # ADD GEOS CONTACT SOLVER REFERENCE DATA (only once) + # =================================================================== + + if geosContactData is not None: + # Format: [Depth_m, Normal_Stress_bar, Shear_Stress_bar, SCU] + # Index: [0, 1, 2, 3] + + print( " 📊 Adding GEOS Contact Solver reference data..." ) + + # Normal stress + axes[ 0 ].plot( geosContactData[ :, 1 ], + geosContactData[ :, 0 ], + marker='o', + color='black', + markersize=7, + label='GEOS Contact Solver', + linestyle='none', + alpha=0.8, + mec='black', + mew=1.5, + fillstyle='none' ) + + # Shear stress + axes[ 1 ].plot( geosContactData[ :, 2 ], + geosContactData[ :, 0 ], + marker='o', + color='black', + markersize=7, + label='GEOS Contact Solver', + linestyle='none', + alpha=0.8, + mec='black', + mew=1.5, + fillstyle='none' ) + + # SCU (if available) + if geosContactData.shape[ 1 ] > 3: + axes[ 2 ].plot( geosContactData[ :, 3 ], + geosContactData[ :, 0 ], + marker='o', + color='black', + markersize=7, + label='GEOS Contact Solver', + linestyle='none', + alpha=0.8, + mec='black', + mew=1.5, + fillstyle='none' ) + + # =================================================================== + # CONFIGURE PLOTS + # =================================================================== + + fsize = 14 + + # Plot 0: Normal Stress + axes[ 0 ].set_xlabel( 'Normal Stress σₙ [bar]', fontsize=fsize, weight="bold" ) + axes[ 0 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 0 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 0 ].legend( loc='best', fontsize=fsize - 2 ) + axes[ 0 ].tick_params( labelsize=fsize - 1 ) + + # Plot 1: Shear Stress + axes[ 1 ].set_xlabel( 'Shear Stress τ [bar]', fontsize=fsize, weight="bold" ) + axes[ 1 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 1 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 1 ].legend( loc='best', fontsize=fsize - 2 ) + axes[ 1 ].tick_params( labelsize=fsize - 1 ) + + # Plot 2: SCU + axes[ 2 ].set_xlabel( 'SCU [-]', fontsize=fsize, weight="bold" ) + axes[ 2 ].set_ylabel( 'Depth [m]', fontsize=fsize, weight="bold" ) + axes[ 2 ].axvline( x=0.8, color='forestgreen', linestyle='--', linewidth=2, alpha=0.5, label='Critical (0.8)' ) + axes[ 2 ].axvline( x=1.0, color='red', linestyle='--', linewidth=2, alpha=0.5, label='Failure (1.0)' ) + axes[ 2 ].grid( True, alpha=0.3, linestyle='--' ) + axes[ 2 ].legend( loc='upper right', fontsize=fsize - 2, ncol=1 ) + axes[ 2 ].tick_params( labelsize=fsize - 1 ) + axes[ 2 ].set_xlim( left=0 ) + + # Overall title + years = time / ( 365.25 * 24 * 3600 ) + fig.suptitle( f'Analytical (Anderson) vs Numerical (GEOS Continuum & Contact) - t={years:.1f} years', + fontsize=fsize + 2, + fontweight='bold', + y=0.995 ) + + # Change verticale scale + if self.maxDepthProfiles is not None: + for i in range( len( axes ) ): + axes[ i ].set_ylim( bottom=self.maxDepthProfiles ) + + if self.minDepthProfiles is not None: + for i in range( len( axes ) ): + axes[ i ].set_ylim( top=self.minDepthProfiles ) + + plt.tight_layout( rect=( 0, 0, 1, 0.99 ) ) + + # Save + if save: + filename = f'analytical_vs_numerical_comparison_{years:.0f}y.png' + plt.savefig( path / filename, dpi=300, bbox_inches='tight' ) + print( f"\n 💾 Comparison plot saved: {filename}" ) + + # Show + if show: + plt.show() + else: + plt.close()