refactor: Move and refactor Surface Geomechanics vtk filter and plugin (#133)

* Move SurfaceGeomechanics

* Refactor the basis change computation

* Refactor of normal and tangential vectors computation

* Add tests and better error handling

Author: paloma-martinez

docs/geos_posp_docs/PVplugins.rst

@@ -62,9 +62,3 @@ PVMohrCirclePlot plugin
 ---------------------------------
 
 .. automodule:: PVplugins.PVMohrCirclePlot
-
-PVplugins.PVSurfaceGeomechanics module
---------------------------------------
-
-.. automodule:: PVplugins.PVSurfaceGeomechanics
-

docs/geos_posp_docs/filters.rst

@@ -10,11 +10,3 @@ geos_posp.filters.GeosBlockMerge module
    :members:
    :undoc-members:
    :show-inheritance:
-
-geos_posp.filters.SurfaceGeomechanics module
-------------------------------------------------
-
-.. automodule:: geos_posp.filters.SurfaceGeomechanics
-   :members:
-   :undoc-members:
-   :show-inheritance:

docs/geos_processing_docs/post_processing.rst

@@ -6,6 +6,28 @@ This module contains GEOS post-processing tools.
 Geomechanics post-processing tools
 =====================================
 
+GEOS computes many outputs including flow and geomechanic properties if coupled flow geomechanics simulations were run. Users however need additional metrics to asses geomechanical stability as part of operational studies. Two filters have been developped to compute these additional metrics: `Geos Geomechanics Calculator`` and `Geos Surfacic Geomechanics`. In addition, a third filter allows to plot Mohr's circles on selected cells and time steps.
+
+.. Note::
+
+    Several processing filters require the definition of physical parameters. The following list is non-exhaustive.
+
+    Default values:
+    - grainBulkModulus = 38e9 Pa ( quartz value )
+    - specificDensity = 1000.0 kg/m³ ( water value )
+    - rock cohesion: 0.0 Pa $( fractured case )
+    - friction angle: 10°
+
+
+geos.processing.post_processing.SurfaceGeomechanics
+-------------------------------------------------------
+
+.. automodule:: geos.processing.post_processing.SurfaceGeomechanics
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
 geos.processing.post_processing.GeomechanicsCalculator module
 --------------------------------------------------------------
 
@@ -14,6 +36,7 @@ geos.processing.post_processing.GeomechanicsCalculator module
    :undoc-members:
    :show-inheritance:
 
+
 geos.processing.post_processing.GeosBlockExtractor module
 ----------------------------------------------------------
 

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

@@ -2,14 +2,34 @@
 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
 # SPDX-FileContributor: Martin Lemay, Paloma Martinez
 import numpy as np
+import logging
 import numpy.typing as npt
-from typing import Iterator, List, Sequence, Any, Union
-from vtkmodules.util.numpy_support import numpy_to_vtk
-from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints, reference
-from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkMultiBlockDataSet, vtkPolyData, vtkDataSet, vtkDataObject, vtkPlane, vtkCellTypes, vtkIncrementalOctreePointLocator
-from vtkmodules.vtkFiltersCore import vtk3DLinearGridPlaneCutter
+from typing import Iterator, List, Sequence, Any, Union, Tuple
+from vtkmodules.util.numpy_support import ( numpy_to_vtk, vtk_to_numpy )
+from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints, reference, vtkDataArray, vtkLogger, vtkFloatArray
+from vtkmodules.vtkCommonDataModel import (
+    vtkUnstructuredGrid,
+    vtkMultiBlockDataSet,
+    vtkPolyData,
+    vtkDataSet,
+    vtkDataObject,
+    vtkPlane,
+    vtkCellTypes,
+    vtkIncrementalOctreePointLocator,
+)
+from vtkmodules.vtkFiltersCore import ( vtk3DLinearGridPlaneCutter, vtkPolyDataNormals, vtkPolyDataTangents )
+from vtkmodules.vtkFiltersTexture import vtkTextureMapToPlane
+
 from geos.mesh.utils.multiblockHelpers import ( getBlockElementIndexesFlatten, getBlockFromFlatIndex )
 
+from geos.utils.algebraFunctions import (
+    getAttributeMatrixFromVector,
+    getAttributeVectorFromMatrix,
+)
+from geos.utils.geometryFunctions import ( getChangeOfBasisMatrix, CANONICAL_BASIS_3D )
+from geos.utils.Logger import ( getLogger, Logger, VTKCaptureLog, RegexExceptionFilter )
+from geos.utils.Errors import VTKError
+
 __doc__ = """
 Generic VTK utilities.
 
@@ -267,3 +287,262 @@ def createVertices( cellPtsCoord: list[ npt.NDArray[ np.float64 ] ],
             cellVertexMap += [ ptId.get() ]  # type: ignore
         cellVertexMapAll += [ tuple( cellVertexMap ) ]
     return points, cellVertexMapAll
+
+
+def convertAttributeFromLocalToXYZForOneCell(
+    vector: npt.NDArray[ np.float64 ], localBasisVectors: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ],
+                                                                 npt.NDArray[ np.float64 ] ]
+) -> npt.NDArray[ np.float64 ]:
+    """Convert one cell attribute from local to XYZ basis.
+
+    .. Warning::
+        Vectors components are considered to be in GEOS output order such that \
+        v = ( M00, M11, M22, M12, M02, M01, M21, M20, M10 )
+
+    Args:
+        vector (npt.NDArray[np.float64]): The attribute expressed in local basis. The size of the vector must be 3, 9 or 6 (symmetrical case)
+        localBasisVectors (tuple[ npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]]): The local basis vectors expressed in the XYZ basis.
+
+    Returns:
+        vectorXYZ (npt.NDArray[np.float64]): The attribute expressed in XYZ basis.
+    """
+    # Get components matrix from GEOS attribute vector
+    matrix3x3: npt.NDArray[ np.float64 ] = getAttributeMatrixFromVector( vector )
+
+    # Get transform matrix
+    transformMatrix: npt.NDArray[ np.float64 ] = getChangeOfBasisMatrix( localBasisVectors, CANONICAL_BASIS_3D )
+
+    # Apply transformation
+    arrayXYZ: npt.NDArray[ np.float64 ] = transformMatrix @ matrix3x3 @ transformMatrix.T
+
+    # Convert back to GEOS type attribute and return
+    return getAttributeVectorFromMatrix( arrayXYZ, vector.size )
+
+
+def getNormalVectors( surface: vtkPolyData ) -> npt.NDArray[ np.float64 ]:
+    """Return the normal vectors of a surface mesh.
+
+    Args:
+        surface (vtkPolyData): The input surface.
+
+    Raises:
+        ValueError: No normal attribute found in the mesh. Use the computeNormals function beforehand.
+
+    Returns:
+        npt.NDArray[ np.float64 ]: The normal vectors of the input surface.
+    """
+    normals: Union[ npt.NDArray[ np.float64 ],
+                    None ] = surface.GetCellData().GetNormals()  # type: ignore[no-untyped-call]
+    if normals is None:
+        raise ValueError( "No normal attribute found in the mesh. Use the computeNormals function beforehand." )
+
+    return vtk_to_numpy( normals )
+
+
+def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ]:
+    """Return the tangential vectors of a surface.
+
+    Args:
+        surface (vtkPolyData): The input surface.
+
+    Raises:
+        ValueError: Tangent attribute not found in the mesh. Use computeTangents beforehand.
+        ValueError: Problem during the calculation of the second tangential vectors.
+
+    Returns:
+        Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: The tangents vectors of the input surface.
+    """
+    # Get first tangential component
+    vtkTangents: Union[ vtkFloatArray, None ] = surface.GetCellData().GetTangents()  # type: ignore[no-untyped-call]
+    tangents1: npt.NDArray[ np.float64 ]
+
+    try:
+        tangents1 = vtk_to_numpy( vtkTangents )
+    except AttributeError as err:
+        print( "No tangential attribute found in the mesh. Use the computeTangents function beforehand." )
+        raise VTKError( err ) from err
+    else:
+        # Compute second tangential component
+        normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface )
+
+        tangents2: npt.NDArray[ np.float64 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 )
+
+    return ( tangents1, tangents2 )
+
+
+def getLocalBasisVectors( surface: vtkPolyData ) -> npt.NDArray[ np.float64 ]:
+    """Return the local basis vectors for all cells of the input surface.
+
+    Args:
+        surface(vtkPolydata): The input surface.
+
+    Returns:
+        npt.NDArray[np.float64]: Array with normal, tangential 1 and tangential 2 vectors.
+    """
+    try:
+        normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface )
+        surfaceWithNormals: vtkPolyData = surface
+    # ValueError raised if no normals found in the mesh
+    except ValueError:
+        # In that case, the normals are computed.
+        surfaceWithNormals = computeNormals( surface )
+        normals = getNormalVectors( surfaceWithNormals )
+
+    # Tangents require normals to be present in the mesh
+    try:
+        tangents: Tuple[ npt.NDArray[ np.float64 ],
+                         npt.NDArray[ np.float64 ] ] = getTangentsVectors( surfaceWithNormals )
+    # If no tangents is present in the mesh, they are computed on that surface
+    except VTKError:
+        surfaceWithTangents: vtkPolyData = computeTangents( surfaceWithNormals )
+        tangents = getTangentsVectors( surfaceWithTangents )
+
+    return np.array( ( normals, *tangents ) )
+
+
+def computeNormals(
+    surface: vtkPolyData,
+    logger: Union[ Logger, None ] = None,
+) -> vtkPolyData:
+    """Compute and set the normals of a given surface.
+
+    Args:
+        surface (vtkPolyData): The input surface.
+        logger (Union[Logger, None], optional): A logger to manage the output messages.
+            Defaults to None, an internal logger is used.
+
+    Returns:
+        vtkPolyData: The surface with normal attribute.
+    """
+    if logger is None:
+        logger = getLogger( "computeSurfaceNormals" )
+    # Creation of a child logger to deal with VTKErrors without polluting parent logger
+    vtkErrorLogger: Logger = getLogger( f"{logger.name}.vtkErrorLogger" )
+    vtkErrorLogger.propagate = False
+
+    vtkLogger.SetStderrVerbosity( vtkLogger.VERBOSITY_ERROR )
+
+    vtkErrorLogger.addFilter( RegexExceptionFilter() )  # will raise VTKError if captured VTK Error
+    vtkErrorLogger.setLevel( logging.DEBUG )
+
+    with VTKCaptureLog() as capturedLog:
+        normalFilter: vtkPolyDataNormals = vtkPolyDataNormals()
+        normalFilter.SetInputData( surface )
+        normalFilter.ComputeCellNormalsOn()
+        normalFilter.ComputePointNormalsOff()
+        normalFilter.Update()
+
+        capturedLog.seek( 0 )
+        captured = capturedLog.read().decode()
+
+        vtkErrorLogger.debug( captured.strip() )
+
+    outputSurface = normalFilter.GetOutput()
+
+    if outputSurface.GetCellData().GetNormals() is None:
+        raise VTKError( "Something went wrong with VTK calculation of the normals." )
+
+    return outputSurface
+
+
+def computeTangents(
+    triangulatedSurface: vtkPolyData,
+    logger: Union[ Logger, None ] = None,
+) -> vtkPolyData:
+    """Compute and set the tangents of a given surface.
+
+    .. Warning:: The computation of tangents requires a triangulated surface.
+
+    Args:
+        triangulatedSurface (vtkPolyData): The input surface. It should be triangulated beforehand.
+        logger (Union[Logger, None], optional): A logger to manage the output messages.
+            Defaults to None, an internal logger is used.
+
+    Returns:
+        vtkPolyData: The surface with tangent attribute
+    """
+    # need to compute texture coordinates required for tangent calculation
+    surface1: vtkPolyData = computeSurfaceTextureCoordinates( triangulatedSurface )
+
+    # TODO: triangulate the surface before computation of the tangents if needed.
+
+    # compute tangents
+    if logger is None:
+        logger = getLogger( "computeSurfaceTangents" )
+    # Creation of a child logger to deal with VTKErrors without polluting parent logger
+    vtkErrorLogger: Logger = getLogger( f"{logger.name}.vtkErrorLogger" )
+    vtkErrorLogger.propagate = False
+
+    vtkLogger.SetStderrVerbosity( vtkLogger.VERBOSITY_ERROR )
+
+    vtkErrorLogger.addFilter( RegexExceptionFilter() )  # will raise VTKError if captured VTK Error
+    vtkErrorLogger.setLevel( logging.DEBUG )
+
+    with VTKCaptureLog() as capturedLog:
+
+        tangentFilter: vtkPolyDataTangents = vtkPolyDataTangents()
+        tangentFilter.SetInputData( surface1 )
+        tangentFilter.ComputeCellTangentsOn()
+        tangentFilter.ComputePointTangentsOff()
+        tangentFilter.Update()
+        surfaceOut: vtkPolyData = tangentFilter.GetOutput()
+
+        capturedLog.seek( 0 )
+        captured = capturedLog.read().decode()
+
+        vtkErrorLogger.debug( captured.strip() )
+
+    if surfaceOut is None:
+        raise VTKError( "Something went wrong in VTK calculation." )
+
+    # copy tangents attributes into filter input surface because surface attributes
+    # are not transferred to the output (bug that should be corrected by Kitware)
+    array: vtkDataArray = surfaceOut.GetCellData().GetTangents()
+    if array is None:
+        raise VTKError( "Attribute Tangents is not in the mesh." )
+
+    surface1.GetCellData().SetTangents( array )
+    surface1.GetCellData().Modified()
+    surface1.Modified()
+    return surface1
+
+
+def computeSurfaceTextureCoordinates(
+    surface: vtkPolyData,
+    logger: Union[ Logger, None ] = None,
+) -> vtkPolyData:
+    """Generate the 2D texture coordinates required for tangent computation. The dataset points are mapped onto a plane generated automatically.
+
+    Args:
+        surface (vtkPolyData): The input surface.
+        logger (Union[Logger, None], optional): A logger to manage the output messages.
+            Defaults to None, an internal logger is used.
+
+    Returns:
+        vtkPolyData: The input surface with generated texture map.
+    """
+    # Need to compute texture coordinates required for tangent calculation
+    if logger is None:
+        logger = getLogger( "computeSurfaceTextureCoordinates" )
+    # Creation of a child logger to deal with VTKErrors without polluting parent logger
+    vtkErrorLogger: Logger = getLogger( f"{logger.name}.vtkErrorLogger" )
+    vtkErrorLogger.propagate = False
+
+    vtkLogger.SetStderrVerbosity( vtkLogger.VERBOSITY_ERROR )
+
+    vtkErrorLogger.addFilter( RegexExceptionFilter() )  # will raise VTKError if captured VTK Error
+    vtkErrorLogger.setLevel( logging.DEBUG )
+
+    with VTKCaptureLog() as capturedLog:
+
+        textureFilter: vtkTextureMapToPlane = vtkTextureMapToPlane()
+        textureFilter.SetInputData( surface )
+        textureFilter.AutomaticPlaneGenerationOn()
+        textureFilter.Update()
+
+        capturedLog.seek( 0 )
+        captured = capturedLog.read().decode()
+
+        vtkErrorLogger.debug( captured.strip() )
+
+    return textureFilter.GetOutput()

geos-mesh/tests/conftest.py

@@ -1,6 +1,6 @@
 # SPDX-License-Identifier: Apache-2.0
 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
-# SPDX-FileContributor: Paloma Martinez
+# SPDX-FileContributor: Paloma Martinez, Romain Baville
 # SPDX-License-Identifier: Apache 2.0
 # ruff: noqa: E402 # disable Module level import not at top of file
 import os