diff --git a/docs/geos_posp_docs/PVplugins.rst b/docs/geos_posp_docs/PVplugins.rst
index 8f41728b..391c2632 100644
--- a/docs/geos_posp_docs/PVplugins.rst
+++ b/docs/geos_posp_docs/PVplugins.rst
@@ -62,9 +62,3 @@ PVMohrCirclePlot plugin
---------------------------------
.. automodule:: PVplugins.PVMohrCirclePlot
-
-PVplugins.PVSurfaceGeomechanics module
---------------------------------------
-
-.. automodule:: PVplugins.PVSurfaceGeomechanics
-
diff --git a/docs/geos_posp_docs/filters.rst b/docs/geos_posp_docs/filters.rst
index 350c3391..7bd6a88d 100644
--- a/docs/geos_posp_docs/filters.rst
+++ b/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:
diff --git a/docs/geos_processing_docs/post_processing.rst b/docs/geos_processing_docs/post_processing.rst
index 23891fbb..58224d9f 100644
--- a/docs/geos_processing_docs/post_processing.rst
+++ b/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
----------------------------------------------------------
diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index fb168d44..0e2ebb39 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/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()
diff --git a/geos-mesh/tests/conftest.py b/geos-mesh/tests/conftest.py
index 8c62234c..1f99f961 100644
--- a/geos-mesh/tests/conftest.py
+++ b/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
diff --git a/geos-mesh/tests/test_genericHelpers_normalAndTangents.py b/geos-mesh/tests/test_genericHelpers_normalAndTangents.py
new file mode 100644
index 00000000..5cfc3f27
--- /dev/null
+++ b/geos-mesh/tests/test_genericHelpers_normalAndTangents.py
@@ -0,0 +1,302 @@
+# SPDX-FileContributor: Paloma Martinez
+# SPDX-License-Identifier: Apache 2.0
+# ruff: noqa: E402 # disable Module level import not at top of file
+import pytest
+from enum import IntEnum
+from typing_extensions import Self, Iterator
+import numpy.typing as npt
+import numpy as np
+from dataclasses import dataclass, field
+
+from vtkmodules.vtkCommonDataModel import ( vtkPolyData, vtkTriangle, vtkCellArray )
+from vtkmodules.vtkCommonCore import vtkPoints
+from vtkmodules.util.numpy_support import vtk_to_numpy
+
+from geos.mesh.utils.genericHelpers import ( computeSurfaceTextureCoordinates, computeTangents, computeNormals,
+ getLocalBasisVectors, getTangentsVectors, getNormalVectors,
+ convertAttributeFromLocalToXYZForOneCell )
+
+from geos.utils.Errors import VTKError
+
+# yapf: disable
+pointsCoordsAll: list[ list[ list[ float ] ] ] = [
+ [ [ 0., 0., 0. ], [ 0., 1., 0. ], [ 0., 2., 0. ], [ 0., 2., 1. ], [ 0., 1., 1. ], [ 0., 0., 1. ], ],
+ [ [ 0., 0., 0. ], [ 0., 1., 0. ], [ 0., 2., 0. ], [ 0., 2., 1. ], [ 1., 1., 1.5 ], [ 0., 0., 1. ], ],
+]
+trianglesAll: list[ list[ tuple[ int, int, int ] ] ] = [
+ [ ( 0, 1, 5 ), ( 1, 4, 5 ), ( 1, 2, 3 ), ( 1, 4, 3 ) ],
+ [ ( 0, 1, 5 ), ( 1, 4, 5 ), ( 1, 2, 3 ), ( 1, 4, 3 ) ],
+]
+expectedNormalsAll: list[ npt.NDArray[ np.float64 ] ] = [
+ np.array( [ [ 1., 0., 0. ], [ 1., 0., 0. ], [ 1., 0., 0. ], [ 1., 0., 0. ], ] ),
+ np.array( [ [ 1., 0., 0. ], [ 0.7276069, -0.48507124, -0.48507124 ], [ 1., 0., 0. ],
+ [ 0.7276069, 0.48507124, -0.48507124 ] ] ),
+]
+expectedTangentsAll: list[ npt.NDArray[ np.float64 ] ] = [
+ np.array( [ [ [ 0., 2., 0. ], [ -0., 2., 0. ], [ 0., 2., 0. ], [ 0., 2., 0. ] ],
+ [ [ 0., 0., 2. ], [ 0., -0., 2. ], [ 0., 0., 2. ], [ 0., 0., 2. ],
+ ] ] ),
+ np.array( [ [ [ 0., 2., 0. ], [ 0.8301887, 2., -0.754717 ], [ 0., 2., 0. ], [ -0.8301887, 2., 0.754717 ] ],
+ [ [ 0., 0., 2. ], [ 1.33623397, 0.14643663, 1.85791445 ], [ 0., 0., 2. ],
+ [ 1.33623397, -0.14643663, 1.85791445 ] ] ] ),
+]
+expectedTexturedCoordsAll: list[ npt.NDArray[ np.float64 ] ] = [
+ np.array( [ [ 0., 0. ], [ 0.5, 0. ], [ 1., 0. ], [ 1., 1. ], [ 0.5, 1. ], [ 0., 1. ],
+ ] ),
+ np.array( [ [ [ 0., 0. ], [ 0.5, 0. ], [ 1., 0. ], [ 1., 0.41509435 ], [ 0.5, 1. ], [ 0., 0.41509435 ] ],
+ ] ),
+]
+# Same attribute for all cells, test for vector size 3, 6, 9
+attributes: list[ npt.NDArray[ np.float64 ] ] = [ np.arange( 1., 4., dtype=np.float64 ), np.arange( 1., 7., dtype=np.float64 ), np.arange( 1., 10., dtype=np.float64 )
+]
+expectedAttributesXYZ: list[ list[ npt.NDArray[ np.float64 ] ] ] = [
+ [
+ np.array( [ [ 1., 8., 12. ], [ 1., 8., 12. ], [ 1., 8., 12. ], [ 1., 8., 12. ],
+ ] ),
+ np.array( [ [ 1., 8., 12., 16., 10., 12. ], [ 1., 8., 12., 16., 10., 12. ], [ 1., 8., 12., 16., 10., 12. ],
+ [ 1., 8., 12., 16., 10., 12. ] ] ),
+ np.array( [ [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ], [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ],
+ [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ], [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ],
+ [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ], [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ] ] )
+ ], # case 1
+ [
+ np.array( [ [ 1., 8., 12. ], [ 7.26440201, 8.29962517, 11.73002787 ], [ 1., 8., 12. ],
+ [ 7.26440201, 8.29962517, 11.73002787 ] ] ),
+ np.array( [ [ 1., 8., 12., 16., 10., 12. ],
+ [ 33.11015535, -1.70942051, -4.10667954, 3.96829788, 5.78481908, 18.33796323 ],
+ [ 1., 8., 12., 16., 10., 12. ],
+ [ 0.86370966, 16.88802688, 9.54231792, 17.62557587, 12.93534579, 16.64449814 ] ] ),
+ np.array( [ [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ],
+ [ 41.16704654, -3.95432477, -9.91866643, 0.51321919, -0.39322437, 23.20275907, 13.51039649,
+ 12.82015999, 23.3879596
+ ], [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ],
+ [ -1.35966323, 18.70673795, 9.9469796, 14.59669037, 15.22437721, 25.39830602, 32.57499969,
+ 14.01099303, 21.0552047
+ ] ] )
+ ], # case 2
+]
+# yapf: enable
+
+
+class Options( IntEnum ):
+ """Enum corresponding to data required for tests.
+
+ - RAW : bare mesh
+ - NORMALS : only normals are calculated
+ - NORMALS_TEXTURED : normals and textured coordinates are computed
+ - TANGENTS : normals and tangents are present in the mesh
+ """
+ RAW = 0
+ NORMALS = 1
+ NORMALS_TEXTURED = 2
+ TANGENTS = 3
+
+
+@dataclass
+class TriangulatedSurfaceTestCase:
+ """Test case."""
+ pointsCoords: list[ list[ float ] ]
+ triangles: list[ tuple[ int, int, int ] ]
+ attribute: list[ npt.NDArray ]
+ expectedNormals: npt.NDArray
+ expectedTangents: npt.NDArray
+ expectedTexturedCoords: npt.NDArray
+ expectedAttributeInXYZ: list[ npt.NDArray ]
+ options: Options
+ mesh: vtkPolyData = field( init=False )
+
+ def __post_init__( self: Self ) -> None:
+ """Generate the mesh."""
+ points: vtkPoints = vtkPoints()
+ for coord in self.pointsCoords:
+ points.InsertNextPoint( coord )
+
+ triangles: vtkCellArray = vtkCellArray()
+ for t in self.triangles:
+ triangle = vtkTriangle()
+ triangle.GetPointIds().SetId( 0, t[ 0 ] )
+ triangle.GetPointIds().SetId( 1, t[ 1 ] )
+ triangle.GetPointIds().SetId( 2, t[ 2 ] )
+ triangles.InsertNextCell( triangle )
+
+ polydata: vtkPolyData = vtkPolyData()
+ polydata.SetPoints( points )
+ polydata.SetPolys( triangles )
+
+ # Compute additional properties depending on the tests options
+ if self.options == Options.NORMALS:
+ self.mesh = computeNormals( polydata )
+
+ elif self.options == Options.NORMALS_TEXTURED:
+ mesh: vtkPolyData = computeSurfaceTextureCoordinates( polydata )
+ mesh2: vtkPolyData = computeNormals( mesh )
+ self.mesh = mesh2
+
+ elif self.options == Options.TANGENTS:
+ mesh = computeSurfaceTextureCoordinates( polydata )
+ mesh2 = computeNormals( mesh )
+ mesh3: vtkPolyData = computeTangents( mesh2 )
+ self.mesh = mesh3
+
+ else:
+ # Unknown cases and case 0
+ self.mesh = polydata
+
+
+def __generateSurfacicTestCase( options: tuple[ Options, ...] ) -> Iterator[ TriangulatedSurfaceTestCase ]:
+ """Generate test cases with different options.
+
+ Args:
+ options (tuple[Options]): Requested additional feature.
+
+ Yields:
+ Iterator[ TriangulatedSurfaceTestCase ]: Test case containing mesh with requested optional features and expected values.
+ """
+ for opt in options:
+ for i in range( len( pointsCoordsAll ) ):
+ yield TriangulatedSurfaceTestCase( pointsCoordsAll[ i ], trianglesAll[ i ], attributes,
+ expectedNormalsAll[ i ], expectedTangentsAll[ i ],
+ expectedTexturedCoordsAll[ i ], expectedAttributesXYZ[ i ], opt )
+
+
+@pytest.mark.parametrize( "case", __generateSurfacicTestCase( options=( Options.RAW, ) ) )
+def test_computeTextureCoords( case: TriangulatedSurfaceTestCase ) -> None:
+ """Test the computation of texture coordinates."""
+ stc: vtkPolyData = computeSurfaceTextureCoordinates( case.mesh )
+ texturedMap: npt.NDArray[ np.float64 ] = vtk_to_numpy( stc.GetPointData().GetArray( "Texture Coordinates" ) )
+ assert np.allclose( texturedMap, case.expectedTexturedCoords )
+
+
+@pytest.mark.parametrize( "case", __generateSurfacicTestCase( options=( Options.NORMALS_TEXTURED, ) ) )
+def test_computeTangents( case: TriangulatedSurfaceTestCase ) -> None:
+ """Test the computation of tangents."""
+ surfaceWithTangents: vtkPolyData = computeTangents( case.mesh )
+ tangents: npt.NDArray[ np.float64 ] = vtk_to_numpy( surfaceWithTangents.GetCellData().GetTangents() )
+
+ assert np.allclose( tangents, case.expectedTangents[ 0 ] )
+
+
+@pytest.mark.parametrize( "case", __generateSurfacicTestCase( options=( Options.TANGENTS, ) ) )
+def test_getTangents( case: TriangulatedSurfaceTestCase ) -> None:
+ """Test tangents getter."""
+ tangents1: npt.NDArray[ np.float64 ]
+ tangents2: npt.NDArray[ np.float64 ]
+
+ tangents1, tangents2 = getTangentsVectors( case.mesh )
+
+ assert np.allclose( tangents1, case.expectedTangents[ 0 ], rtol=1e-3 )
+ assert np.allclose( tangents2, case.expectedTangents[ 1 ], rtol=1e-3 )
+
+
+@pytest.mark.parametrize( "case", __generateSurfacicTestCase( options=( Options.RAW, ) ) )
+def test_computeNormals( case: TriangulatedSurfaceTestCase ) -> None:
+ """Test the computation of normals."""
+ surfaceWithNormals = computeNormals( case.mesh )
+
+ normals = vtk_to_numpy( surfaceWithNormals.GetCellData().GetNormals() )
+
+ assert np.allclose( normals, case.expectedNormals, rtol=1e-3 )
+
+
+@pytest.mark.parametrize( "case", __generateSurfacicTestCase( options=( Options.NORMALS, ) ) )
+def test_getNormals( case: TriangulatedSurfaceTestCase ) -> None:
+ """Test normals getter."""
+ normals = getNormalVectors( case.mesh )
+
+ assert np.allclose( normals, case.expectedNormals, rtol=1e-3 )
+
+
+@pytest.mark.parametrize( "case",
+ __generateSurfacicTestCase( options=(
+ Options.RAW,
+ Options.NORMALS,
+ Options.TANGENTS,
+ ) ) )
+def test_getLocalBasis( case: TriangulatedSurfaceTestCase ) -> None:
+ """Test local basis getter."""
+ normals, tangents1, tangents2 = getLocalBasisVectors( case.mesh )
+
+ assert np.allclose( normals, case.expectedNormals, rtol=1e-3 )
+ assert np.allclose( tangents1, case.expectedTangents[ 0 ], rtol=1e-3 )
+ assert np.allclose( tangents2, case.expectedTangents[ 1 ], rtol=1e-3 )
+
+
+#########################################################################################
+
+
+@pytest.fixture( scope="function" )
+def emptySurface() -> vtkPolyData:
+ """Generate and return an empty vtkPolyData for tests.
+
+ Returns:
+ vtkPolyData: empty vtkPolyData
+ """
+ return vtkPolyData()
+
+
+def test_failingComputeSurfaceTextureCoords( emptySurface: vtkPolyData ) -> None:
+ """Test VTK error raising of texture coordinate calculation."""
+ with pytest.raises( VTKError ):
+ computeSurfaceTextureCoordinates( emptySurface )
+
+
+def test_failingComputeTangents( emptySurface: vtkPolyData ) -> None:
+ """Test VTK error raising during tangent calculation."""
+ with pytest.raises( VTKError ):
+ computeTangents( emptySurface )
+
+
+def test_failingComputeNormals( emptySurface: vtkPolyData ) -> None:
+ """Test VTK error raising of normals calculation."""
+ with pytest.raises( VTKError ):
+ computeNormals( emptySurface )
+
+
+def test_failingGetTangents( emptySurface: vtkPolyData ) -> None:
+ """Test error raising when getting the surface tangents."""
+ with pytest.raises( VTKError ):
+ getTangentsVectors( emptySurface )
+
+
+def test_failingGetNormals( emptySurface: vtkPolyData ) -> None:
+ """Test error raising when getting the surface normals."""
+ with pytest.raises( ValueError ):
+ getNormalVectors( emptySurface )
+
+
+###############################################################
+
+
+@dataclass( frozen=True )
+class AttributeConversionTestCase:
+ """Test case for attribute conversion from local basis to XYZ basis for one cell."""
+ vector: npt.NDArray[ np.float64 ]
+ normal: npt.NDArray[ np.float64 ]
+ tangent1: npt.NDArray[ np.float64 ]
+ tangent2: npt.NDArray[ np.float64 ]
+ expectedVectorXYZ: npt.NDArray[ np.float64 ]
+
+
+def __generateAttributeConversionTestCase() -> Iterator[ AttributeConversionTestCase ]:
+ """Generator of test cases for the conversion of attributes from local to XYZ attributes."""
+ cases: Iterator[ TriangulatedSurfaceTestCase ] = __generateSurfacicTestCase( options=( Options.RAW, ) )
+ for nc, testcase in enumerate( cases ):
+ if nc == 0:
+ # Try vector size 3, 6 and 9
+ for nvec, localAttribute in enumerate( attributes ):
+ for ncell in range( testcase.mesh.GetNumberOfCells() ):
+
+ yield AttributeConversionTestCase( localAttribute, testcase.expectedNormals[ ncell ],
+ testcase.expectedTangents[ 0 ][ ncell ],
+ testcase.expectedTangents[ 1 ][ ncell ],
+ testcase.expectedAttributeInXYZ[ nvec ][ ncell ] )
+
+
+@pytest.mark.parametrize( "testcase", __generateAttributeConversionTestCase() )
+def test_convertAttributesToXYZ( testcase: AttributeConversionTestCase ) -> None:
+ """Test the conversion of one cell attribute from local to canonic basis."""
+ localAttr: npt.NDArray[ np.float64 ] = testcase.vector
+
+ attrXYZ: npt.NDArray[ np.float64 ] = convertAttributeFromLocalToXYZForOneCell(
+ localAttr, ( testcase.normal, testcase.tangent1, testcase.tangent2 ) )
+ assert np.allclose( attrXYZ, testcase.expectedVectorXYZ, rtol=1e-6 )
diff --git a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py
index 8b559688..03024b83 100644
--- a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py
+++ b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py
@@ -36,7 +36,7 @@
from PVplugins.PVExtractMergeBlocksVolumeSurface import (
PVExtractMergeBlocksVolumeSurface, )
from geos.processing.post_processing.GeomechanicsCalculator import GeomechanicsCalculator
-from PVplugins.PVSurfaceGeomechanics import PVSurfaceGeomechanics
+from geos.pv.plugins.PVSurfaceGeomechanics import PVSurfaceGeomechanics
__doc__ = """
PVGeomechanicsWorkflowVolumeSurface is a Paraview plugin that execute
diff --git a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py
index 715de74c..52cae47b 100644
--- a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py
+++ b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py
@@ -36,7 +36,7 @@
from PVplugins.PVExtractMergeBlocksVolumeSurfaceWell import (
PVExtractMergeBlocksVolumeSurfaceWell, )
from geos.processing.post_processing.GeomechanicsCalculator import GeomechanicsCalculator
-from PVplugins.PVSurfaceGeomechanics import PVSurfaceGeomechanics
+from geos.pv.plugins.PVSurfaceGeomechanics import PVSurfaceGeomechanics
__doc__ = """
PVGeomechanicsWorkflowVolumeSurfaceWell is a Paraview plugin that execute
diff --git a/geos-posp/src/PVplugins/PVSurfaceGeomechanics.py b/geos-posp/src/PVplugins/PVSurfaceGeomechanics.py
deleted file mode 100644
index 0d0b7ed5..00000000
--- a/geos-posp/src/PVplugins/PVSurfaceGeomechanics.py
+++ /dev/null
@@ -1,250 +0,0 @@
-# SPDX-License-Identifier: Apache-2.0
-# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
-# SPDX-FileContributor: Martin Lemay
-# ruff: noqa: E402 # disable Module level import not at top of file
-import os
-import sys
-
-import numpy as np
-from typing_extensions import Self
-
-dir_path = os.path.dirname( os.path.realpath( __file__ ) )
-parent_dir_path = os.path.dirname( dir_path )
-if parent_dir_path not in sys.path:
- sys.path.append( parent_dir_path )
-
-import PVplugins # noqa: F401
-
-from geos.utils.Logger import Logger, getLogger
-from geos.utils.PhysicalConstants import (
- DEFAULT_FRICTION_ANGLE_DEG,
- DEFAULT_ROCK_COHESION,
-)
-from geos_posp.filters.SurfaceGeomechanics import SurfaceGeomechanics
-from geos.mesh.utils.multiblockHelpers import (
- getBlockElementIndexesFlatten,
- getBlockFromFlatIndex,
-)
-from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found]
- VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy,
-)
-from vtkmodules.vtkCommonCore import (
- vtkDataArray,
- vtkInformation,
- vtkInformationVector,
-)
-from vtkmodules.vtkCommonDataModel import (
- vtkDataObject,
- vtkMultiBlockDataSet,
- vtkPolyData,
-)
-
-__doc__ = r"""
-PVSurfaceGeomechanics is a Paraview plugin that allows to compute
-additional geomechanical attributes from the input surfaces.
-
-Input and output types are vtkMultiBlockDataSet.
-
-To use it:
-
-* Load the module in Paraview: Tools>Manage Plugins...>Load new>PVSurfaceGeomechanics.
-* Select any pipeline child of the second ouput from
- GeosExtractMergeBlocksVolumeSurface* filter.
-* Search and Apply PVSurfaceGeomechanics Filter.
-
-"""
-
-
-@smproxy.filter( name="PVSurfaceGeomechanics", label="Geos Surface Geomechanics" )
-@smhint.xml( '' )
-@smproperty.input( name="Input", port_index=0 )
-@smdomain.datatype( dataTypes=[ "vtkMultiBlockDataSet" ], composite_data_supported=True )
-class PVSurfaceGeomechanics( VTKPythonAlgorithmBase ):
-
- def __init__( self: Self ) -> None:
- """Paraview plugin to compute additional geomechanical surface outputs.
-
- Input is either a vtkMultiBlockDataSet that contains surfaces with
- Normals and Tangential attributes.
- """
- super().__init__(
- nInputPorts=1,
- nOutputPorts=1,
- inputType="vtkMultiBlockDataSet",
- outputType="vtkMultiBlockDataSet",
- )
-
- # rock cohesion (Pa)
- self.m_rockCohesion: float = DEFAULT_ROCK_COHESION
- # friction angle (°)
- self.m_frictionAngle: float = DEFAULT_FRICTION_ANGLE_DEG
- # logger
- self.m_logger: Logger = getLogger( "Surface Geomechanics Filter" )
-
- def SetLogger( self: Self, logger: Logger ) -> None:
- """Set filter logger.
-
- Args:
- logger (Logger): logger
- """
- self.m_logger = logger
-
- @smproperty.doublevector(
- name="RockCohesion",
- label="Rock Cohesion (Pa)",
- default_values=DEFAULT_ROCK_COHESION,
- panel_visibility="default",
- )
- @smdomain.xml( """
-
- Reference rock cohesion to compute critical pore pressure.
- The unit is Pa. Default is fractured case (i.e., 0. Pa).
-
- """ )
- def a01SetRockCohesion( self: Self, value: float ) -> None:
- """Set rock cohesion.
-
- Args:
- value (float): rock cohesion (Pa)
- """
- self.m_rockCohesion = value
- self.Modified()
-
- def getRockCohesion( self: Self ) -> float:
- """Get rock cohesion.
-
- Returns:
- float: rock cohesion.
- """
- return self.m_rockCohesion
-
- @smproperty.doublevector(
- name="FrictionAngle",
- label="Friction Angle (°)",
- default_values=DEFAULT_FRICTION_ANGLE_DEG,
- panel_visibility="default",
- )
- @smdomain.xml( """
-
- Reference friction angle to compute critical pore pressure.
- The unit is °. Default is no friction case (i.e., 0.°).
-
- """ )
- def a02SetFrictionAngle( self: Self, value: float ) -> None:
- """Set frition angle.
-
- Args:
- value (float): friction angle (°)
- """
- self.m_frictionAngle = value
- self.Modified()
-
- def getFrictionAngle( self: Self ) -> float:
- """Get friction angle in radian.
-
- Returns:
- float: friction angle.
- """
- return self.m_frictionAngle * np.pi / 180.0
-
- def FillInputPortInformation( self: Self, port: int, info: vtkInformation ) -> int:
- """Inherited from VTKPythonAlgorithmBase::RequestInformation.
-
- Args:
- port (int): input port
- info (vtkInformationVector): info
-
- Returns:
- int: 1 if calculation successfully ended, 0 otherwise.
- """
- if port == 0:
- info.Set( self.INPUT_REQUIRED_DATA_TYPE(), "vtkMultiBlockDataSet" )
- return 1
-
- def RequestInformation(
- self: Self,
- request: vtkInformation, # noqa: F841
- inInfoVec: list[ vtkInformationVector ], # noqa: F841
- outInfoVec: vtkInformationVector,
- ) -> int:
- """Inherited from VTKPythonAlgorithmBase::RequestInformation.
-
- Args:
- request (vtkInformation): request
- inInfoVec (list[vtkInformationVector]): input objects
- outInfoVec (vtkInformationVector): output objects
-
- Returns:
- int: 1 if calculation successfully ended, 0 otherwise.
- """
- executive = self.GetExecutive() # noqa: F841
- outInfo = outInfoVec.GetInformationObject( 0 ) # noqa: F841
- return 1
-
- def RequestData(
- self: Self,
- request: vtkInformation, # noqa: F841
- inInfoVec: list[ vtkInformationVector ],
- outInfoVec: vtkInformationVector,
- ) -> int:
- """Inherited from VTKPythonAlgorithmBase::RequestData.
-
- Args:
- request (vtkInformation): request
- inInfoVec (list[vtkInformationVector]): input objects
- outInfoVec (vtkInformationVector): output objects
-
- Returns:
- int: 1 if calculation successfully ended, 0 otherwise.
- """
- self.m_logger.info( f"Apply filter {__name__}" )
- try:
- input0: vtkMultiBlockDataSet = vtkMultiBlockDataSet.GetData( inInfoVec[ 0 ] )
- output: vtkMultiBlockDataSet = self.GetOutputData( outInfoVec, 0 )
-
- assert input0 is not None, "Input Surface is null."
- assert output is not None, "Output pipeline is null."
-
- output.ShallowCopy( input0 )
- self.computeSurfaceGeomecanics( input0, output )
- output.Modified()
- mess: str = ( "Surface geomechanics attributes calculation successfully ended." )
- self.m_logger.info( mess )
- except AssertionError as e:
- mess1: str = "Surface geomechanics attributes calculation failed due to:"
- self.m_logger.error( mess1 )
- self.m_logger.error( e, exc_info=True )
- return 0
- except Exception as e:
- mess0: str = "Surface geomechanics attributes calculation failed due to:"
- self.m_logger.critical( mess0 )
- self.m_logger.critical( e, exc_info=True )
- return 0
- return 1
-
- def computeSurfaceGeomecanics( self: Self, input: vtkMultiBlockDataSet, output: vtkMultiBlockDataSet ) -> None:
- """Compute surface geomechanics new attributes.
-
- Args:
- input (vtkMultiBlockDataSet): input multiBlockDataSet
- output (vtkMultiBlockDataSet): output multiBlockDataSet
- """
- surfaceBlockIndexes: list[ int ] = getBlockElementIndexesFlatten( input )
- for blockIndex in surfaceBlockIndexes:
- surfaceBlock0: vtkDataObject = getBlockFromFlatIndex( output, blockIndex )
- assert surfaceBlock0 is not None, "Surface is undefined."
- surfaceBlock: vtkPolyData = vtkPolyData.SafeDownCast( surfaceBlock0 )
- filter: SurfaceGeomechanics = SurfaceGeomechanics()
- filter.AddInputDataObject( surfaceBlock )
- filter.SetRockCohesion( self.getRockCohesion() )
- filter.SetFrictionAngle( self.getFrictionAngle() )
- filter.Update()
- outputSurface: vtkPolyData = filter.GetOutputDataObject( 0 )
-
- # add attributes to output surface mesh
- for attributeName in filter.GetNewAttributeNames():
- attr: vtkDataArray = outputSurface.GetCellData().GetArray( attributeName )
- surfaceBlock.GetCellData().AddArray( attr )
- surfaceBlock.GetCellData().Modified()
- surfaceBlock.Modified()
- output.Modified()
diff --git a/geos-posp/src/geos_posp/filters/GeosBlockMerge.py b/geos-posp/src/geos_posp/filters/GeosBlockMerge.py
index 56021ecb..4b76d9a2 100644
--- a/geos-posp/src/geos_posp/filters/GeosBlockMerge.py
+++ b/geos-posp/src/geos_posp/filters/GeosBlockMerge.py
@@ -14,7 +14,6 @@
from typing_extensions import Self
from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase
from vtkmodules.vtkCommonCore import (
- vtkDataArray,
vtkInformation,
vtkInformationVector,
)
@@ -25,19 +24,18 @@
vtkPolyData,
vtkUnstructuredGrid,
)
-from vtkmodules.vtkFiltersCore import (
- vtkArrayRename,
- vtkPolyDataNormals,
- vtkPolyDataTangents,
-)
+from vtkmodules.vtkFiltersCore import vtkArrayRename
from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter
-from vtkmodules.vtkFiltersTexture import vtkTextureMapToPlane
from geos.mesh.utils.multiblockHelpers import getElementaryCompositeBlockIndexes
from geos.mesh.utils.arrayHelpers import getAttributeSet
from geos.mesh.utils.arrayModifiers import createConstantAttribute, fillAllPartialAttributes
from geos.mesh.utils.multiblockHelpers import extractBlock
from geos.mesh.utils.multiblockModifiers import mergeBlocks
+from geos.mesh.utils.genericHelpers import (
+ computeNormals,
+ computeTangents,
+)
__doc__ = """
GeosBlockMerge module is a vtk filter that allows to merge Geos ranks, rename
@@ -377,10 +375,10 @@ def convertFaultsToSurfaces( self: Self ) -> bool:
surface1: vtkPolyData = self.convertBlockToSurface( surface0 )
assert surface1 is not None, "Surface extraction from block failed."
# compute normals
- surface2: vtkPolyData = self.computeNormals( surface1 )
+ surface2: vtkPolyData = computeNormals( surface1 )
assert surface2 is not None, "Normal calculation failed."
# compute tangents
- surface3: vtkPolyData = self.computeTangents( surface2 )
+ surface3: vtkPolyData = computeTangents( surface2 )
assert surface3 is not None, "Tangent calculation failed."
# set surface to output multiBlockDataSet
@@ -417,56 +415,3 @@ def convertBlockToSurface( self: Self, block: vtkUnstructuredGrid ) -> vtkPolyDa
extractSurfaceFilter.Update()
output: vtkPolyData = extractSurfaceFilter.GetOutput()
return output
-
- def computeNormals( self: Self, surface: vtkPolyData ) -> vtkPolyData:
- """Compute normals of the given surface.
-
- Args:
- surface (vtkPolyData): surface to compute the normals
-
- Returns:
- vtkPolyData: surface with normal attribute
- """
- normalFilter: vtkPolyDataNormals = vtkPolyDataNormals()
- normalFilter.SetInputData( surface )
- normalFilter.ComputeCellNormalsOn()
- normalFilter.ComputePointNormalsOff()
- normalFilter.Update()
- output: vtkPolyData = normalFilter.GetOutput()
- return output
-
- def computeTangents( self: Self, surface: vtkPolyData ) -> vtkPolyData:
- """Compute tangents of the given surface.
-
- Args:
- surface (vtkPolyData): surface to compute the tangents
-
- Returns:
- vtkPolyData: surface with tangent attribute
- """
- # need to compute texture coordinates required for tangent calculation
- textureFilter: vtkTextureMapToPlane = vtkTextureMapToPlane()
- textureFilter.SetInputData( surface )
- textureFilter.AutomaticPlaneGenerationOn()
- textureFilter.Update()
- surface1: vtkPolyData = textureFilter.GetOutput()
- assert ( surface1 is not None ), "Texture calculation during Tangent calculation failed."
-
- # compute tangents
- tangentFilter: vtkPolyDataTangents = vtkPolyDataTangents()
- tangentFilter.SetInputData( surface1 )
- tangentFilter.ComputeCellTangentsOn()
- tangentFilter.ComputePointTangentsOff()
- tangentFilter.Update()
- surfaceOut: vtkPolyData = tangentFilter.GetOutput()
- assert surfaceOut is not None, "Tangent components calculation failed."
-
- # 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()
- assert array is not None, "Attribute Tangents is not in the mesh."
-
- surface1.GetCellData().SetTangents( array )
- surface1.GetCellData().Modified()
- surface1.Modified()
- return surface1
diff --git a/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py b/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py
deleted file mode 100644
index 7e50cbf7..00000000
--- a/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py
+++ /dev/null
@@ -1,457 +0,0 @@
-# SPDX-License-Identifier: Apache-2.0
-# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
-# SPDX-FileContributor: Martin Lemay
-# ruff: noqa: E402 # disable Module level import not at top of file
-import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts
-import geos.utils.geometryFunctions as geom
-import numpy as np
-import numpy.typing as npt
-import vtkmodules.util.numpy_support as vnp
-from geos.utils.algebraFunctions import (
- getAttributeMatrixFromVector,
- getAttributeVectorFromMatrix,
-)
-from geos.utils.GeosOutputsConstants import (
- ComponentNameEnum,
- GeosMeshOutputsEnum,
- PostProcessingOutputsEnum,
- getAttributeToConvertFromLocalToXYZ,
-)
-from geos.utils.Logger import Logger, getLogger
-from geos.utils.PhysicalConstants import (
- DEFAULT_FRICTION_ANGLE_RAD,
- DEFAULT_ROCK_COHESION,
-)
-from typing_extensions import Self
-from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase
-from vtkmodules.vtkCommonCore import (
- vtkDataArray,
- vtkDoubleArray,
- vtkInformation,
- vtkInformationVector,
-)
-from vtkmodules.vtkCommonDataModel import (
- vtkPolyData, )
-
-from geos.mesh.utils.arrayModifiers import createAttribute
-from geos.mesh.utils.arrayHelpers import (
- getArrayInObject,
- getAttributeSet,
- isAttributeInObject,
-)
-
-__doc__ = """
-SurfaceGeomechanics module is a vtk filter that allows to compute Geomechanical
-properties on surfaces.
-
-Filter input and output types are vtkPolyData.
-
-To use the filter:
-
-.. code-block:: python
-
- from filters.SurfaceGeomechanics import SurfaceGeomechanics
-
- # filter inputs
- logger :Logger
- input :vtkPolyData
- rockCohesion :float
- frictionAngle :float # angle in radian
-
- # instanciate the filter
- surfaceGeomechanicsFilter :SurfaceGeomechanics = SurfaceGeomechanics()
- # set the logger
- surfaceGeomechanicsFilter.SetLogger(logger)
- # set input data object
- surfaceGeomechanicsFilter.SetInputDataObject(input)
- # set rock cohesion anf friction angle
- surfaceGeomechanicsFilter.SetRockCohesion(rockCohesion)
- surfaceGeomechanicsFilter.SetFrictionAngle(frictionAngle)
- # do calculations
- surfaceGeomechanicsFilter.Update()
- # get output object
- output :vtkPolyData = surfaceGeomechanicsFilter.GetOutputDataObject(0)
- # get created attribute names
- newAttributeNames :set[str] = surfaceGeomechanicsFilter.GetNewAttributeNames()
-"""
-
-
-class SurfaceGeomechanics( VTKPythonAlgorithmBase ):
-
- def __init__( self: Self ) -> None:
- """Vtk filter to compute geomechanical surfacic attributes.
-
- Input and Output objects are a vtkMultiBlockDataSet containing surface
- objects with Normals and Tangential attributes.
- """
- super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkPolyData" ) # type: ignore[no-untyped-call]
-
- # output surface mesh
- self.m_output: vtkPolyData
- # attributes are either on points or on cells
- self.m_attributeOnPoints: bool = False
- # rock cohesion (Pa)
- self.m_rockCohesion: float = DEFAULT_ROCK_COHESION
- # friction angle (rad)
- self.m_frictionAngle: float = DEFAULT_FRICTION_ANGLE_RAD
- # new created attributes names
- self.m_newAttributeNames: set[ str ] = set()
-
- # logger
- self.m_logger: Logger = getLogger( "Surface Geomechanics Filter" )
-
- def SetRockCohesion( self: Self, rockCohesion: float ) -> None:
- """Set rock cohesion value.
-
- Args:
- rockCohesion (float): rock cohesion (Pa)
- """
- self.m_rockCohesion = rockCohesion
-
- def SetFrictionAngle( self: Self, frictionAngle: float ) -> None:
- """Set friction angle value.
-
- Args:
- frictionAngle (float): friction angle (rad)
- """
- self.m_frictionAngle = frictionAngle
-
- def SetLogger( self: Self, logger: Logger ) -> None:
- """Set the logger.
-
- Args:
- logger (Logger): logger
- """
- self.m_logger = logger
-
- def GetNewAttributeNames( self: Self ) -> set[ str ]:
- """Get the set of new attribute names that were created.
-
- Returns:
- set[str]: set of new attribute names
- """
- return self.m_newAttributeNames
-
- def FillInputPortInformation( self: Self, port: int, info: vtkInformation ) -> int:
- """Inherited from VTKPythonAlgorithmBase::RequestInformation.
-
- Args:
- port (int): input port
- info (vtkInformationVector): info
-
- Returns:
- int: 1 if calculation successfully ended, 0 otherwise.
- """
- if port == 0:
- info.Set( self.INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData" )
- return 1
-
- def RequestInformation(
- self: Self,
- request: vtkInformation, # noqa: F841
- inInfoVec: list[ vtkInformationVector ], # noqa: F841
- outInfoVec: vtkInformationVector,
- ) -> int:
- """Inherited from VTKPythonAlgorithmBase::RequestInformation.
-
- Args:
- request (vtkInformation): request
- inInfoVec (list[vtkInformationVector]): input objects
- outInfoVec (vtkInformationVector): output objects
-
- Returns:
- int: 1 if calculation successfully ended, 0 otherwise.
- """
- executive = self.GetExecutive() # noqa: F841
- outInfo = outInfoVec.GetInformationObject( 0 ) # noqa: F841
- return 1
-
- def RequestData(
- self: Self,
- request: vtkInformation, # noqa: F841
- inInfoVec: list[ vtkInformationVector ],
- outInfoVec: vtkInformationVector,
- ) -> int:
- """Inherited from VTKPythonAlgorithmBase::RequestData.
-
- Args:
- request (vtkInformation): request
- inInfoVec (list[vtkInformationVector]): input objects
- outInfoVec (vtkInformationVector): output objects
-
- Returns:
- int: 1 if calculation successfully ended, 0 otherwise.
- """
- try:
- input = vtkPolyData.GetData( inInfoVec[ 0 ] )
- self.m_output = self.GetOutputData( outInfoVec, 0 ) # type: ignore[no-untyped-call]
-
- assert input is not None, "Input object is null."
- assert self.m_output is not None, "Output pipeline is null."
-
- self.m_output.ShallowCopy( input )
- # conversion of vectorial attributes from Normal/Tangent basis to xyz basis
- assert ( self.convertLocalToXYZBasisAttributes() ), "Error while converting Local to XYZ basis attributes"
- # compute shear capacity utilization
- assert self.computeShearCapacityUtilization(), "Error while computing SCU."
- self.Modified()
- except AssertionError as e:
- mess: str = "Surface geomechanics attributes calculation failed due to:"
- self.m_logger.error( mess )
- self.m_logger.error( e, exc_info=True )
- return 0
- except Exception as e:
- mess1: str = "Surface geomechanics attributes calculation failed due to:"
- self.m_logger.critical( mess1 )
- self.m_logger.critical( e, exc_info=True )
- return 0
- mess2: str = "Surface geomechanics attributes were successfully computed."
- self.m_logger.info( mess2 )
- return 1
-
- def convertLocalToXYZBasisAttributes( self: Self ) -> bool:
- """Convert vectorial property coordinates from Local to canonic basis.
-
- Returns:
- bool: True if calculation successfully ended, False otherwise
- """
- # look for the list of attributes to convert
- attributesToConvert: set[ str ] = self.getAttributesToConvertFromLocalToXYZ()
- if len( attributesToConvert ) == 0:
- self.m_logger.info( "No attribute to convert from local to XYZ basis were found." )
- return False
-
- # get local coordinate vectors
- normalTangentVectors: npt.NDArray[ np.float64 ] = self.getNormalTangentsVectors()
- # do conversion for each attribute
- for attributName in attributesToConvert:
- # skip attribute if it is already in the object
- newAttrName: str = attributName + "_" + ComponentNameEnum.XYZ.name
- if isAttributeInObject( self.m_output, newAttrName, False ):
- continue
-
- attr: vtkDoubleArray = self.m_output.GetCellData().GetArray( attributName )
- assert attr is not None, "Attribute {attributName} is undefined."
- assert attr.GetNumberOfComponents() > 2, ( "Dimension of the attribute " +
- " must be equal or grater than 3." )
-
- attrArray: npt.NDArray[ np.float64 ] = vnp.vtk_to_numpy( attr ) # type: ignore[no-untyped-call]
- newAttrArray: npt.NDArray[ np.float64 ] = self.computeNewCoordinates( attrArray, normalTangentVectors,
- True )
-
- # create attribute
- createAttribute(
- self.m_output,
- newAttrArray,
- newAttrName,
- ComponentNameEnum.XYZ.value,
- False,
- )
- self.m_newAttributeNames.add( newAttrName )
- return True
-
- def convertXYZToLocalBasisAttributes( self: Self ) -> bool:
- """Convert vectorial property coordinates from canonic to local basis.
-
- Returns:
- bool: True if calculation successfully ended, False otherwise
- """
- # look for the list of attributes to convert
- # empty but to update if needed in the future
- attributesToConvert: set[ str ] = set()
- if len( attributesToConvert ) == 0:
- self.m_logger.info( "No attribute to convert from local to " + "canonic basis were found." )
- return False
-
- # get local coordinate vectors
- normalTangentVectors: npt.NDArray[ np.float64 ] = self.getNormalTangentsVectors()
- for attributName in attributesToConvert:
- # skip attribute if it is already in the object
- newAttrName: str = ( attributName + "_" + ComponentNameEnum.NORMAL_TANGENTS.name )
- if isAttributeInObject( self.m_output, newAttrName, False ):
- continue
- attr: vtkDoubleArray = self.m_output.GetCellData().GetArray( attributName )
- assert attr is not None, "Attribute {attributName} is undefined."
- assert attr.GetNumberOfComponents() > 2, ( "Dimension of the attribute " +
- " must be equal or grater than 3." )
-
- attrArray: npt.NDArray[ np.float64 ] = vnp.vtk_to_numpy( attr ) # type: ignore[no-untyped-call]
- newAttrArray: npt.NDArray[ np.float64 ] = self.computeNewCoordinates( attrArray, normalTangentVectors,
- False )
-
- # create attribute
- createAttribute(
- self.m_output,
- newAttrArray,
- newAttrName,
- ComponentNameEnum.NORMAL_TANGENTS.value,
- False,
- )
- self.m_newAttributeNames.add( newAttrName )
- return True
-
- def getAttributesToConvertFromLocalToXYZ( self: Self ) -> set[ str ]:
- """Get the list of attributes to convert from local to XYZ basis.
-
- Returns:
- set[str]: Set of the attribute names.
- """
- return self.filterAttributesToConvert( getAttributeToConvertFromLocalToXYZ() )
-
- def filterAttributesToConvert( self: Self, attributesToFilter0: set[ str ] ) -> set[ str ]:
- """Filter the set of attribute names if they are vectorial and present.
-
- Args:
- attributesToFilter0 (set[str]): set of attribute names to filter.
-
- Returns:
- set[str]: Set of the attribute names.
- """
- attributesFiltered: set[ str ] = set()
- attributeSet: set[ str ] = getAttributeSet( self.m_output, False )
- for attributName in attributesToFilter0:
- if attributName in attributeSet:
- attr: vtkDataArray = self.m_output.GetCellData().GetArray( attributName )
- if attr.GetNumberOfComponents() > 2:
- attributesFiltered.add( attributName )
- return attributesFiltered
-
- def computeNewCoordinates(
- self: Self,
- attrArray: npt.NDArray[ np.float64 ],
- normalTangentVectors: npt.NDArray[ np.float64 ],
- fromLocalToYXZ: bool,
- ) -> npt.NDArray[ np.float64 ]:
- """Compute the coordinates of a vectorial attribute.
-
- Args:
- attrArray (npt.NDArray[np.float64]): vector of attribute values
- normalTangentVectors (npt.NDArray[np.float64]): 3xNx3 local vector
- coordinates
- fromLocalToYXZ (bool): if True, conversion is done from local to XYZ
- basis, otherwise conversion is done from XZY to Local basis.
-
- Returns:
- npt.NDArray[np.float64]: Vector of new coordinates of the attribute.
- """
- attrArrayNew = np.full_like( attrArray, np.nan )
- # for each cell
- for i in range( attrArray.shape[ 0 ] ):
- # get the change of basis matrix
- localBasis: npt.NDArray[ np.float64 ] = normalTangentVectors[ :, i, : ]
- changeOfBasisMatrix = self.computeChangeOfBasisMatrix( localBasis, fromLocalToYXZ )
- if attrArray.shape[ 1 ] == 3:
- attrArrayNew[ i ] = self.computeNewCoordinatesVector3( attrArray[ i ], changeOfBasisMatrix )
- else:
- attrArrayNew[ i ] = self.computeNewCoordinatesVector6( attrArray[ i ], changeOfBasisMatrix )
-
- assert np.any( np.isfinite( attrArrayNew ) ), "Attribute new coordinate calculation failed."
- return attrArrayNew
-
- def computeNewCoordinatesVector3(
- self: Self,
- vector: npt.NDArray[ np.float64 ],
- changeOfBasisMatrix: npt.NDArray[ np.float64 ],
- ) -> npt.NDArray[ np.float64 ]:
- """Compute attribute new coordinates of vector of size 3.
-
- Args:
- vector (npt.NDArray[np.float64]): input coordinates.
- changeOfBasisMatrix (npt.NDArray[np.float64]): change of basis matrix
-
- Returns:
- npt.NDArray[np.float64]: new coordinates
- """
- return geom.computeCoordinatesInNewBasis( vector, changeOfBasisMatrix )
-
- def computeNewCoordinatesVector6(
- self: Self,
- vector: npt.NDArray[ np.float64 ],
- changeOfBasisMatrix: npt.NDArray[ np.float64 ],
- ) -> npt.NDArray[ np.float64 ]:
- """Compute attribute new coordinates of vector of size > 3.
-
- Args:
- vector (npt.NDArray[np.float64]): input coordinates.
- changeOfBasisMatrix (npt.NDArray[np.float64]): change of basis matrix
-
- Returns:
- npt.NDArray[np.float64]: new coordinates
- """
- attributeMatrix: npt.NDArray[ np.float64 ] = getAttributeMatrixFromVector( vector )
- attributeMatrixNew: npt.NDArray[ np.float64 ] = np.full_like( attributeMatrix, np.nan )
- # for each column of the matrix
- for j in range( attributeMatrix.shape[ 1 ] ):
- attributeMatrixNew[ :, j ] = geom.computeCoordinatesInNewBasis( attributeMatrix[ :, j ],
- changeOfBasisMatrix )
- return getAttributeVectorFromMatrix( attributeMatrixNew, vector.size )
-
- def computeChangeOfBasisMatrix( self: Self, localBasis: npt.NDArray[ np.float64 ],
- fromLocalToYXZ: bool ) -> npt.NDArray[ np.float64 ]:
- """Compute the change of basis matrix according to local coordinates.
-
- Args:
- localBasis (npt.NDArray[np.float64]): local coordinate vectors.
- fromLocalToYXZ (bool): if True, change of basis matrix is from local
- to XYZ bases, otherwise it is from XYZ to local bases.
-
- Returns:
- npt.NDArray[np.float64]: change of basis matrix.
- """
- P: npt.NDArray[ np.float64 ] = np.transpose( localBasis )
- if fromLocalToYXZ:
- return P
- # inverse the change of basis matrix
- return np.linalg.inv( P ).astype( np.float64 )
-
- def getNormalTangentsVectors( self: Self ) -> npt.NDArray[ np.float64 ]:
- """Compute the change of basis matrix from Local to XYZ bases.
-
- Returns:
- npt.NDArray[np.float64]: Nx3 matrix of local vector coordinates.
- """
- # get normal and first tangent components
- normals: npt.NDArray[ np.float64 ] = vnp.vtk_to_numpy(
- self.m_output.GetCellData().GetNormals() ) # type: ignore[no-untyped-call]
- assert normals is not None, "Normal attribute was not found."
- tangents1: npt.NDArray[ np.float64 ] = vnp.vtk_to_numpy(
- self.m_output.GetCellData().GetTangents() ) # type: ignore[no-untyped-call]
- assert tangents1 is not None, "Tangents attribute was not found."
-
- # compute second tangential component
- tangents2: npt.NDArray[ np.float64 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 )
- assert tangents2 is not None, "Local basis third axis was not computed."
-
- # put vectors as columns
- return np.array( ( normals, tangents1, tangents2 ) )
-
- def computeShearCapacityUtilization( self: Self ) -> bool:
- """Compute the shear capacity utilization on surface.
-
- Returns:
- bool: True if calculation successfully ended, False otherwise.
- """
- try:
- SCUAttributeName: str = PostProcessingOutputsEnum.SCU.attributeName
- if not isAttributeInObject( self.m_output, SCUAttributeName, self.m_attributeOnPoints ):
- tractionAttributeName: str = GeosMeshOutputsEnum.TRACTION.attributeName
- traction: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, tractionAttributeName,
- self.m_attributeOnPoints )
- assert traction is not None, ( f"{tractionAttributeName}" + " attribute is undefined." )
-
- scuAttribute: npt.NDArray[ np.float64 ] = fcts.shearCapacityUtilization(
- traction, self.m_rockCohesion, self.m_frictionAngle )
- createAttribute(
- self.m_output,
- scuAttribute,
- SCUAttributeName,
- (),
- self.m_attributeOnPoints,
- )
- self.m_newAttributeNames.add( SCUAttributeName )
-
- except AssertionError as e:
- self.m_logger.error( "Shear Capacity Utilization was not computed due to:" )
- self.m_logger.error( str( e ) )
- return False
- return True
diff --git a/geos-processing/pyproject.toml b/geos-processing/pyproject.toml
index 71fadf3e..3d947f67 100644
--- a/geos-processing/pyproject.toml
+++ b/geos-processing/pyproject.toml
@@ -67,4 +67,4 @@ python_files = "test*.py"
python_functions = "test*"
testpaths = ["tests"]
norecursedirs = "bin"
-filterwarnings = []
\ No newline at end of file
+filterwarnings = []
diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py
new file mode 100644
index 00000000..90320e82
--- /dev/null
+++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py
@@ -0,0 +1,411 @@
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Martin Lemay, Paloma Martinez
+# ruff: noqa: E402 # disable Module level import not at top of file
+import logging
+import numpy as np
+
+from typing_extensions import Self, Union
+import numpy.typing as npt
+
+from vtkmodules.vtkCommonCore import vtkDataArray
+from vtkmodules.vtkCommonDataModel import vtkPolyData
+
+from geos.mesh.utils.arrayModifiers import createAttribute
+from geos.mesh.utils.arrayHelpers import (
+ getArrayInObject,
+ getAttributeSet,
+ isAttributeInObject,
+)
+from geos.mesh.utils.genericHelpers import (
+ getLocalBasisVectors,
+ convertAttributeFromLocalToXYZForOneCell,
+)
+import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts
+from geos.utils.Logger import ( Logger, getLogger )
+from geos.utils.PhysicalConstants import (
+ DEFAULT_FRICTION_ANGLE_RAD,
+ DEFAULT_ROCK_COHESION,
+)
+from geos.utils.GeosOutputsConstants import (
+ ComponentNameEnum,
+ GeosMeshOutputsEnum,
+ PostProcessingOutputsEnum,
+)
+
+__doc__ = """
+SurfaceGeomechanics is a VTK filter that allows:
+ - Conversion of a set of attributes from local basis to XYZ basis
+ - Computation of the shear capacity utilization (SCU)
+
+.. Warning::
+ - The computation of the SCU requires the presence of a 'traction' attribute in the input mesh.
+ - Conversion from local to XYZ basis is currently only handled for cell attributes.
+
+.. Note::
+ Default values for physical constants used in this filter:
+ - rock cohesion: 0.0 Pa ( fractured case)
+ - friction angle: 10°
+
+Filter input and output types are vtkPolyData.
+
+To use the filter:
+
+.. code-block:: python
+
+ from geos.processing.post_processing.SurfaceGeomechanics import SurfaceGeomechanics
+
+ # filter inputs
+ inputMesh: vtkPolyData
+
+ # Optional inputs
+ rockCohesion: float # Pa
+ frictionAngle: float # angle in radian
+ speHandler: bool # defaults to false
+
+ # Instantiate the filter
+ sg: SurfaceGeomechanics = SurfaceGeomechanics( inputMesh, speHandler )
+
+ # Set the handler of yours (only if speHandler is True).
+ yourHandler: logging.Handler
+ sg.SetLoggerHandler( yourHandler )
+
+ # [optional] Set rock cohesion and friction angle
+ sg.SetRockCohesion( rockCohesion )
+ sg.SetFrictionAngle( frictionAngle )
+
+ # Do calculations
+ sg.applyFilter()
+
+ # Get output object
+ output: vtkPolyData = sg.GetOutputMesh()
+
+ # Get created attribute names
+ newAttributeNames: set[str] = sg.GetNewAttributeNames()
+
+
+.. Note::
+
+ By default, conversion of attributes from local to XYZ basis is performed for the following list: { 'displacementJump' }.
+ This list can be modified in different ways:
+
+ - Addition of one or several additional attributes to the set by using the filter function `AddAttributesToConvert`.
+ - Replace the list completely with the function `SetAttributesToConvert`.
+
+ Note that the dimension of the attributes to convert must be equal or greater than 3.
+"""
+loggerTitle: str = "Surface Geomechanics"
+
+
+class SurfaceGeomechanics:
+
+ def __init__( self: Self, surfacicMesh: vtkPolyData, speHandler: bool = False ) -> None:
+ """Vtk filter to compute geomechanical surfacic attributes.
+
+ Input and Output objects are a vtkPolydata with surfaces
+ objects with Normals and Tangential attributes.
+
+ Args:
+ surfacicMesh (vtkPolyData): The input surfacic mesh.
+ speHandler (bool, optional): True to use a specific handler, False to use the internal handler.
+ Defaults to False.
+ """
+ # Logger
+ self.logger: Logger
+ if not speHandler:
+ self.logger = getLogger( loggerTitle, True )
+ else:
+ self.logger = logging.getLogger( loggerTitle )
+ self.logger.setLevel( logging.INFO )
+
+ # Input surfacic mesh
+ if not surfacicMesh.IsA( "vtkPolyData" ):
+ self.logger.error( f"Input surface is expected to be a vtkPolyData, not a {type(surfacicMesh)}." )
+ self.inputMesh: vtkPolyData = surfacicMesh
+ # Identification of the input surface (logging purpose)
+ self.name: Union[ str, None ] = None
+ # Output surfacic mesh
+ self.outputMesh: vtkPolyData
+ # Default attributes to convert from local to XYZ
+ self.convertAttributesOn: bool = True
+ self.attributesToConvert: set[ str ] = {
+ GeosMeshOutputsEnum.DISPLACEMENT_JUMP.attributeName,
+ }
+
+ # Attributes are either on points or on cells
+ self.attributeOnPoints: bool = False
+ # Rock cohesion (Pa)
+ self.rockCohesion: float = DEFAULT_ROCK_COHESION
+ # Friction angle (rad)
+ self.frictionAngle: float = DEFAULT_FRICTION_ANGLE_RAD
+ # New created attributes names
+ self.newAttributeNames: set[ str ] = set()
+
+ def SetLoggerHandler( self: Self, handler: Logger ) -> None:
+ """Set a specific handler for the filter logger.
+
+ In this filter 4 log levels are use, .info, .error, .warning and .critical, be sure to have at least the same 4 levels.
+
+ Args:
+ handler (logging.Handler): The handler to add.
+ """
+ if not self.logger.hasHandlers():
+ self.logger.addHandler( handler )
+ else:
+ self.logger.warning(
+ "The logger already has an handler, to use yours set the argument 'speHandler' to True during the filter initialization."
+ )
+
+ def SetSurfaceName( self: Self, name: str ) -> None:
+ """Set a name for the input surface. For logging purpose only.
+
+ Args:
+ name (str): The identifier for the surface.
+ """
+ self.name = name
+
+ def SetRockCohesion( self: Self, rockCohesion: float ) -> None:
+ """Set rock cohesion value. Defaults to 0.0 Pa.
+
+ Args:
+ rockCohesion (float): The rock cohesion in Pascal.
+ """
+ self.rockCohesion = rockCohesion
+
+ def SetFrictionAngle( self: Self, frictionAngle: float ) -> None:
+ """Set friction angle value in radians. Defaults to 10 / 180 * pi rad.
+
+ Args:
+ frictionAngle (float): The friction angle in radians.
+ """
+ self.frictionAngle = frictionAngle
+
+ def ConvertAttributesOn( self: Self ) -> None:
+ """Activate the conversion of attributes from local to XYZ basis."""
+ self.convertAttributesOn = True
+
+ def ConvertAttributesOff( self: Self ) -> None:
+ """Deactivate the conversion of attributes from local to XYZ bais."""
+ self.convertAttributesOn = False
+
+ def GetConvertAttributes( self: Self ) -> bool:
+ """If convertAttributesOn is True, the set of attributes will be converted from local to XYZ during filter application. Default is True.
+
+ Returns:
+ bool: Current state of the attributes conversion request.
+ """
+ return self.convertAttributesOn
+
+ def GetAttributesToConvert( self: Self ) -> set[ str ]:
+ """Get the set of attributes that will be converted from local to XYZ basis.
+
+ Returns:
+ set[ str ]: The set of attributes that should be converted.
+ """
+ return self.attributesToConvert
+
+ def SetAttributesToConvert( self: Self, attributesToConvert: set[ str ] ) -> None:
+ """Set the list of attributes that will be converted from local to XYZ basis.
+
+ Args:
+ attributesToConvert (set[str]): The set of attributes names that will be converted from local to XYZ basis
+ """
+ self.attributesToConvert = attributesToConvert
+ if len( self.attributesToConvert ) != 0:
+ self.ConvertAttributesOn()
+ else:
+ self.ConvertAttributesOff()
+ self.logger.warning( "Empty set of attributes to convert." )
+
+ def AddAttributesToConvert( self: Self, attributeName: Union[ list[ str ], set[ str ] ] ) -> None:
+ """Add an attribute to the set of attributes to convert.
+
+ Args:
+ attributeName (Union[list[str],set[str]]): List of the attribute array names.
+ """
+ self.attributesToConvert.update( attributeName )
+ self.ConvertAttributesOn()
+
+ def GetNewAttributeNames( self: Self ) -> set[ str ]:
+ """Get the set of new attribute names that were created.
+
+ Returns:
+ set[str]: The set of new attribute names.
+ """
+ return self.newAttributeNames
+
+ def applyFilter( self: Self ) -> bool:
+ """Compute Geomechanical properties on input surface.
+
+ Returns:
+ int: 1 if calculation successfully ended, 0 otherwise.
+ """
+ msg = f"Applying filter {self.logger.name}"
+ if self.name is not None:
+ msg += f" on surface : {self.name}."
+ else:
+ msg += "."
+
+ self.logger.info( msg )
+
+ self.outputMesh = vtkPolyData()
+ self.outputMesh.ShallowCopy( self.inputMesh )
+
+ # Conversion of attributes from Normal/Tangent basis to xyz basis
+ if self.convertAttributesOn:
+ self.logger.info( "Conversion of attributes from local to XYZ basis." )
+ if not self.convertAttributesFromLocalToXYZBasis():
+ self.logger.error( "Error while converting attributes from local to XYZ basis." )
+ return False
+
+ # Compute shear capacity utilization
+ if not self.computeShearCapacityUtilization():
+ self.logger.error( "Error while computing SCU." )
+ return False
+
+ self.logger.info( f"Filter {self.logger.name} successfully applied on surface {self.name}." )
+ return True
+
+ def convertAttributesFromLocalToXYZBasis( self: Self ) -> bool:
+ """Convert attributes from local to XYZ basis.
+
+ Returns:
+ bool: True if calculation successfully ended or no attributes, False otherwise.
+ """
+ # Get the list of attributes to convert and filter
+ attributesToConvert: set[ str ] = self.__filterAttributesToConvert()
+
+ if len( attributesToConvert ) == 0:
+ self.logger.warning( "No attribute to convert from local to XYZ basis were found." )
+ return True
+
+ # Do conversion from local to XYZ for each attribute
+ for attrNameLocal in attributesToConvert:
+ attrNameXYZ: str = f"{attrNameLocal}_{ComponentNameEnum.XYZ.name}"
+
+ # Skip attribute if it is already in the object
+ if isAttributeInObject( self.outputMesh, attrNameXYZ, self.attributeOnPoints ):
+ continue
+
+ if self.attributeOnPoints:
+ self.logger.error(
+ "This filter can only convert cell attributes from local to XYZ basis, not point attributes." )
+ localArray: npt.NDArray[ np.float64 ] = getArrayInObject( self.outputMesh, attrNameLocal,
+ self.attributeOnPoints )
+
+ arrayXYZ: npt.NDArray[ np.float64 ] = self.__computeXYZCoordinates( localArray )
+
+ # Create converted attribute array in dataset
+ if createAttribute( self.outputMesh,
+ arrayXYZ,
+ attrNameXYZ,
+ ComponentNameEnum.XYZ.value,
+ onPoints=self.attributeOnPoints,
+ logger=self.logger ):
+ self.logger.info( f"Attribute {attrNameXYZ} added to the output mesh." )
+ self.newAttributeNames.add( attrNameXYZ )
+
+ return True
+
+ def __filterAttributesToConvert( self: Self ) -> set[ str ]:
+ """Filter the set of attribute names if they are vectorial and present.
+
+ Returns:
+ set[str]: Set of the attribute names.
+ """
+ attributesFiltered: set[ str ] = set()
+
+ if len( self.attributesToConvert ) != 0:
+ attributeSet: set[ str ] = getAttributeSet( self.outputMesh, False )
+ for attrName in self.attributesToConvert:
+ if attrName in attributeSet:
+ attr: vtkDataArray = self.outputMesh.GetCellData().GetArray( attrName )
+ if attr.GetNumberOfComponents() > 2:
+ attributesFiltered.add( attrName )
+ else:
+ self.logger.warning(
+ f"Attribute {attrName} filtered out. Dimension of the attribute must be equal or greater than 3."
+ )
+ else:
+ self.logger.warning( f"Attribute {attrName} not in the input mesh." )
+
+ if len( attributesFiltered ) == 0:
+ self.logger.warning( "All attributes filtered out." )
+ self.ConvertAttributesOff()
+
+ return attributesFiltered
+
+ # TODO: Adapt to handle point attributes.
+ def __computeXYZCoordinates(
+ self: Self,
+ attrArray: npt.NDArray[ np.float64 ],
+ ) -> npt.NDArray[ np.float64 ]:
+ """Compute the XYZ coordinates of a vectorial attribute.
+
+ Args:
+ attrArray (npt.NDArray[np.float64]): vector of attribute values
+
+ Returns:
+ npt.NDArray[np.float64]: Vector of new coordinates of the attribute.
+ """
+ attrXYZ: npt.NDArray[ np.float64 ] = np.full_like( attrArray, np.nan )
+
+ # Get all local basis vectors
+ localBasis: npt.NDArray[ np.float64 ] = getLocalBasisVectors( self.outputMesh )
+
+ for i, cellAttribute in enumerate( attrArray ):
+ if len( cellAttribute ) not in ( 3, 6, 9 ):
+ raise ValueError(
+ f"Inconsistent number of components for attribute. Expected 3, 6 or 9 but got { len( cellAttribute.shape ) }."
+ )
+
+ # Compute attribute XYZ components
+ cellLocalBasis: npt.NDArray[ np.float64 ] = localBasis[ :, i, : ]
+ attrXYZ[ i ] = convertAttributeFromLocalToXYZForOneCell( cellAttribute, cellLocalBasis )
+
+ if not np.any( np.isfinite( attrXYZ ) ):
+ self.logger.error( "Attribute new coordinate calculation failed." )
+
+ return attrXYZ
+
+ def computeShearCapacityUtilization( self: Self ) -> bool:
+ """Compute the shear capacity utilization (SCU) on surface.
+
+ Returns:
+ bool: True if calculation successfully ended, False otherwise.
+ """
+ SCUAttributeName: str = PostProcessingOutputsEnum.SCU.attributeName
+
+ if not isAttributeInObject( self.outputMesh, SCUAttributeName, self.attributeOnPoints ):
+ # Get the traction to compute the SCU
+ tractionAttributeName: str = GeosMeshOutputsEnum.TRACTION.attributeName
+ traction: npt.NDArray[ np.float64 ] = getArrayInObject( self.outputMesh, tractionAttributeName,
+ self.attributeOnPoints )
+
+ # Computation of the shear capacity utilization (SCU)
+ # TODO: better handling of errors in shearCapacityUtilization
+ try:
+ scuAttribute: npt.NDArray[ np.float64 ] = fcts.shearCapacityUtilization(
+ traction, self.rockCohesion, self.frictionAngle )
+ except AssertionError:
+ self.logger.error( f"Failed to compute {SCUAttributeName}." )
+ return False
+
+ # Create attribute
+ if not createAttribute(
+ self.outputMesh, scuAttribute, SCUAttributeName, (), self.attributeOnPoints, logger=self.logger ):
+ self.logger.error( f"Failed to create attribute {SCUAttributeName}." )
+ return False
+ else:
+ self.logger.info( "SCU computed and added to the output mesh." )
+ self.newAttributeNames.add( SCUAttributeName )
+
+ return True
+
+ def GetOutputMesh( self: Self ) -> vtkPolyData:
+ """Get the output mesh with computed attributes.
+
+ Returns:
+ vtkPolyData: The surfacic output mesh.
+ """
+ return self.outputMesh
diff --git a/geos-processing/tests/test_SurfaceGeomechanics.py b/geos-processing/tests/test_SurfaceGeomechanics.py
new file mode 100644
index 00000000..7623e802
--- /dev/null
+++ b/geos-processing/tests/test_SurfaceGeomechanics.py
@@ -0,0 +1,84 @@
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Paloma Martinez
+# SPDX-License-Identifier: Apache 2.0
+# ruff: noqa: E402 # disable Module level import not at top of file
+
+import pytest
+from typing_extensions import Self, Union
+import numpy.typing as npt
+import numpy as np
+
+from dataclasses import dataclass, field
+from vtkmodules.vtkCommonDataModel import vtkPolyData, vtkTriangle, vtkCellArray
+from vtkmodules.vtkCommonCore import vtkPoints, vtkFloatArray
+import vtkmodules.util.numpy_support as vnp
+from geos.processing.post_processing.SurfaceGeomechanics import SurfaceGeomechanics
+
+
+@dataclass
+class TriangulatedSurfaceTestCase:
+ pointsCoords: npt.NDArray[ np.float64 ]
+ triangles: list[ tuple[ int, int, int ] ]
+ mesh: vtkPolyData = field( init=False )
+ attributes: Union[ None, dict[ str, npt.NDArray ] ]
+
+ def __post_init__( self: Self ) -> None:
+ """Generate the mesh."""
+ vtk_points = vtkPoints()
+ for coord in self.pointsCoords:
+ vtk_points.InsertNextPoint( coord )
+
+ vtk_triangles = vtkCellArray()
+ for t in self.triangles:
+ triangle = vtkTriangle()
+ triangle.GetPointIds().SetId( 0, t[ 0 ] )
+ triangle.GetPointIds().SetId( 1, t[ 1 ] )
+ triangle.GetPointIds().SetId( 2, t[ 2 ] )
+ vtk_triangles.InsertNextCell( triangle )
+
+ polydata = vtkPolyData()
+ polydata.SetPoints( vtk_points )
+ polydata.SetPolys( vtk_triangles )
+
+ self.mesh = polydata
+
+ if self.attributes is not None:
+ for attrName, attrValue in self.attributes.items():
+ newAttribute: vtkFloatArray = vnp.numpy_to_vtk( attrValue, deep=True )
+ newAttribute.SetName( attrName )
+
+ self.mesh.GetCellData().AddArray( newAttribute )
+ self.mesh.Modified()
+
+ assert self.mesh.GetCellData().HasArray( attrName )
+
+
+# yapf: disable
+pointsCoords: npt.NDArray[ np.float64 ] = np.array( [ [ 0, 0, 0 ], [ 0, 1, 0 ], [ 0, 2, 0 ], [ 0, 2, 1 ], [ 1, 1, 1.5 ], [ 0, 0, 1 ] ] )
+triangles: list[ tuple[ int, int, int ] ] = [ ( 0, 1, 5 ), ( 1, 4, 5 ), ( 1, 2, 3 ), ( 1, 4, 3 ) ]
+attributes: dict[ str, npt.NDArray ] = { "traction": np.array( [ [ -9e-10, 3e-15, 0.0 ], [ -9e-10, 3e-15, 0.0 ], [ 0.0, 3e-15, -3e-17 ], [ 0.0, 3e-15, -3e-17 ], ] ),
+ "displacementJump": np.array( [ [ 4e-02, 8e-07, 3e-08 ], [ 4e-02, 8e-07, 3e-08 ], [ 2e-02, 4e-07, -8e-08 ], [ 2e-02, 4e-07, -8e-08 ], ] )
+}
+# yapf: enable
+
+
+def test_SurfaceGeomechanics() -> None:
+ """Test SurfaceGeomechanics vtk filter."""
+ testCase: TriangulatedSurfaceTestCase = TriangulatedSurfaceTestCase( pointsCoords, triangles, attributes )
+
+ sgFilter: SurfaceGeomechanics = SurfaceGeomechanics( testCase.mesh )
+
+ assert sgFilter.applyFilter()
+
+ mesh: vtkPolyData = sgFilter.GetOutputMesh()
+ assert mesh.GetCellData().HasArray( "SCU" )
+ assert mesh.GetCellData().HasArray( "displacementJump_XYZ" )
+
+
+def test_failingSurfaceGeomechanics() -> None:
+ """Test failing of SurfaceGeomechanics due to absence of attributes in the mesh."""
+ failingCase: TriangulatedSurfaceTestCase = TriangulatedSurfaceTestCase( pointsCoords, triangles, None )
+ sgFilter: SurfaceGeomechanics = SurfaceGeomechanics( failingCase.mesh )
+ with pytest.raises( AssertionError ):
+ assert sgFilter.applyFilter()
diff --git a/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py
new file mode 100644
index 00000000..2858766e
--- /dev/null
+++ b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py
@@ -0,0 +1,167 @@
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Martin Lemay, Paloma Martinez
+# ruff: noqa: E402 # disable Module level import not at top of file
+import sys
+from pathlib import Path
+import numpy as np
+from typing_extensions import Self
+
+from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found]
+ VTKPythonAlgorithmBase, smdomain, smproperty,
+)
+from paraview.detail.loghandler import ( # type: ignore[import-not-found]
+ VTKHandler, )
+
+# update sys.path to load all GEOS Python Package dependencies
+geos_pv_path: Path = Path( __file__ ).parent.parent.parent.parent.parent
+sys.path.insert( 0, str( geos_pv_path / "src" ) )
+from geos.pv.utils.config import update_paths
+from geos.pv.utils.details import ( SISOFilter, FilterCategory )
+
+update_paths()
+
+from geos.utils.PhysicalConstants import (
+ DEFAULT_FRICTION_ANGLE_DEG,
+ DEFAULT_ROCK_COHESION,
+)
+from geos.processing.post_processing.SurfaceGeomechanics import SurfaceGeomechanics
+from geos.mesh.utils.multiblockHelpers import (
+ getBlockElementIndexesFlatten,
+ getBlockFromFlatIndex,
+)
+from vtkmodules.vtkCommonCore import (
+ vtkDataArray, )
+from vtkmodules.vtkCommonDataModel import (
+ vtkMultiBlockDataSet,
+ vtkPolyData,
+)
+
+__doc__ = """
+PVSurfaceGeomechanics is a Paraview plugin that allows to compute
+additional geomechanical attributes from the input surfaces, such as shear capacity utilization (SCU).
+
+Input and output are vtkMultiBlockDataSet.
+.. Important::
+ - Please refer to the GeosExtractMergeBlockVolumeSurface* filters to provide the correct input.
+ - This filter only works on triangles at the moment. Please apply a triangulation algorithm beforehand if required.
+
+
+To use it:
+
+* Load the module in Paraview: Tools>Manage Plugins...>Load new>PVSurfaceGeomechanics.
+* Select any pipeline child of the second ouput from
+ GeosExtractMergeBlocksVolumeSurface* filter.
+* Select Filters > 3- Geos Geomechanics > Geos Surface Geomechanics.
+* (Optional) Set rock cohesion and/or friction angle.
+* Apply.
+
+"""
+
+
+@SISOFilter( category=FilterCategory.GEOS_GEOMECHANICS,
+ decoratedLabel="Geos Surface Geomechanics",
+ decoratedType="vtkMultiBlockDataSet" )
+class PVSurfaceGeomechanics( VTKPythonAlgorithmBase ):
+
+ def __init__( self: Self ) -> None:
+ """Compute additional geomechanical surface outputs.
+
+ Input is a vtkMultiBlockDataSet containing surfaces.
+ """
+ # rock cohesion (Pa)
+ self.rockCohesion: float = DEFAULT_ROCK_COHESION
+ # friction angle (°)
+ self.frictionAngle: float = DEFAULT_FRICTION_ANGLE_DEG
+
+ @smproperty.doublevector(
+ name="RockCohesion",
+ label="Rock Cohesion (Pa)",
+ default_values=DEFAULT_ROCK_COHESION,
+ panel_visibility="default",
+ )
+ @smdomain.xml( """
+
+ Reference rock cohesion to compute critical pore pressure.
+ The unit is Pa. Default is fractured case (i.e., 0. Pa).
+
+ """ )
+ def a01SetRockCohesion( self: Self, value: float ) -> None:
+ """Set rock cohesion.
+
+ Args:
+ value (float): rock cohesion (Pa)
+ """
+ self.rockCohesion = value
+ self.Modified()
+
+ @smproperty.doublevector(
+ name="FrictionAngle",
+ label="Friction Angle (°)",
+ default_values=DEFAULT_FRICTION_ANGLE_DEG,
+ panel_visibility="default",
+ )
+ @smdomain.xml( """
+
+ Reference friction angle to compute critical pore pressure.
+ The unit is °. Default is no friction case (i.e., 0.°).
+
+ """ )
+ def a02SetFrictionAngle( self: Self, value: float ) -> None:
+ """Set friction angle.
+
+ Args:
+ value (float): friction angle (°)
+ """
+ self.frictionAngle = value
+ self.Modified()
+
+ def Filter( self: Self, inputMesh: vtkMultiBlockDataSet, outputMesh: vtkMultiBlockDataSet ) -> None:
+ """Apply SurfaceGeomechanics filter to the mesh.
+
+ Args:
+ inputMesh (vtkMultiBlockDataSet): The input multiblock mesh with surfaces.
+ outputMesh (vtkMultiBlockDataSet): The output multiblock mesh with converted attributes and SCU.
+ """
+ outputMesh.ShallowCopy( inputMesh )
+
+ surfaceBlockIndexes: list[ int ] = getBlockElementIndexesFlatten( inputMesh )
+ for blockIndex in surfaceBlockIndexes:
+ surfaceBlock: vtkPolyData = vtkPolyData.SafeDownCast( getBlockFromFlatIndex( outputMesh, blockIndex ) )
+
+ sgFilter: SurfaceGeomechanics = SurfaceGeomechanics( surfaceBlock, True )
+ sgFilter.SetSurfaceName( f"blockIndex {blockIndex}" )
+ if not sgFilter.logger.hasHandlers():
+ sgFilter.SetLoggerHandler( VTKHandler() )
+
+ sgFilter.SetRockCohesion( self._getRockCohesion() )
+ sgFilter.SetFrictionAngle( self._getFrictionAngle() )
+ sgFilter.applyFilter()
+
+ outputSurface: vtkPolyData = sgFilter.GetOutputMesh()
+
+ # add attributes to output surface mesh
+ for attributeName in sgFilter.GetNewAttributeNames():
+ attr: vtkDataArray = outputSurface.GetCellData().GetArray( attributeName )
+ surfaceBlock.GetCellData().AddArray( attr )
+ surfaceBlock.GetCellData().Modified()
+ surfaceBlock.Modified()
+
+ outputMesh.Modified()
+ return
+
+ def _getFrictionAngle( self: Self ) -> float:
+ """Get friction angle in radian.
+
+ Returns:
+ float: The friction angle.
+ """
+ return self.frictionAngle * np.pi / 180.0
+
+ def _getRockCohesion( self: Self ) -> float:
+ """Get rock cohesion.
+
+ Returns:
+ float: rock cohesion.
+ """
+ return self.rockCohesion
diff --git a/geos-utils/src/geos/utils/GeosOutputsConstants.py b/geos-utils/src/geos/utils/GeosOutputsConstants.py
index 8a53dca1..a2fc72cb 100644
--- a/geos-utils/src/geos/utils/GeosOutputsConstants.py
+++ b/geos-utils/src/geos/utils/GeosOutputsConstants.py
@@ -300,14 +300,3 @@ def getAttributeToTransferFromInitialTime() -> dict[ str, str ]:
PostProcessingOutputsEnum.POISSON_RATIO.attributeName:
PostProcessingOutputsEnum.POISSON_RATIO_INITIAL.attributeName,
}
-
-
-def getAttributeToConvertFromLocalToXYZ() -> set[ str ]:
- """Get the list of attribute names to convert from local to xyz basis.
-
- Returns:
- list[str]: list of attributes to convert
- """
- return {
- GeosMeshOutputsEnum.DISPLACEMENT_JUMP.attributeName,
- }
diff --git a/geos-utils/src/geos/utils/algebraFunctions.py b/geos-utils/src/geos/utils/algebraFunctions.py
index eeeb88e4..876f8887 100644
--- a/geos-utils/src/geos/utils/algebraFunctions.py
+++ b/geos-utils/src/geos/utils/algebraFunctions.py
@@ -1,12 +1,12 @@
# SPDX-License-Identifier: Apache-2.0
# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
-# SPDX-FileContributor: Martin Lemay
+# SPDX-FileContributor: Martin Lemay, Paloma Martinez
+# ruff: noqa: E402 # disable Module level import not at top of file
import numpy as np
import numpy.typing as npt
__doc__ = """
-GeosUtils module defines usefull methods to process GEOS outputs according to
-GEOS conventions.
+The algebraFunctions module of `geos-utils` package defines methods to switch from vector <> matrix and conversely, following GEOS outputs conventions.
"""
@@ -18,44 +18,69 @@ def getAttributeMatrixFromVector( attrArray: npt.NDArray[ np.float64 ], ) -> npt
* if input vector size is 3:
.. math::
- (v1, v2, v3) => \begin{bmatrix}
+ (v1, v2, v3) => \\begin{bmatrix}
v0 & 0 & 0 \\
0 & v1 & 0 \\
0 & 0 & v2
- \end{bmatrix}
+ \\end{bmatrix}
- * if input vector size is 6:
+ * if input vector size is 9:
.. math::
- (v1, v2, v3, v4, v5, v6) => \begin{bmatrix}
+ (v1, v2, v3, v4, v5, v6, v7, v8, v9) => \\begin{bmatrix}
+ v1 & v6 & v5 \\
+ v9 & v2 & v4 \\
+ v8 & v7 & v3
+ \\end{bmatrix}
+
+ * if input vector size is 6 (symmetrical tensor):
+ .. math::
+
+ (v1, v2, v3, v4, v5, v6) => \\begin{bmatrix}
v1 & v6 & v5 \\
v6 & v2 & v4 \\
v5 & v4 & v3
- \end{bmatrix}
-
+ \\end{bmatrix}
- .. WARNING:: Input vector must be of size 3 or 6.
+ .. Note::
+ In the case of 3 x 3 tensors, GEOS is currently only implemented for symmetrical cases.
Args:
attrArray (npt.NDArray[np.float64]): Vector of attribute values.
- Returns:
- npt.NDArray[np.float64]: matrix of attribute values
+ Raises:
+ ValueError: The input vector must be of size 3, 9 or 6 (symmetrical case).
+ Returns:
+ npt.NDArray[np.float64]: The matrix of attribute values.
"""
- assert attrArray.size > 2, ( "Vectorial attribute must contains at least " + "3 components." )
+ if attrArray.size not in ( 3, 6, 9 ):
+ raise ValueError( "Vectorial attribute must contain 3, 6 or 9 components." )
+
# diagonal terms
matrix: npt.NDArray[ np.float64 ] = np.diagflat( attrArray[ :3 ] )
+
# shear stress components
if attrArray.size == 6:
- matrix[ 0, 1 ] = attrArray[ 5 ]
+ matrix[ 0, 1 ] = attrArray[ 5 ] # v1
matrix[ 1, 0 ] = attrArray[ 5 ]
- matrix[ 0, 2 ] = attrArray[ 4 ]
+ matrix[ 0, 2 ] = attrArray[ 4 ] # v5
matrix[ 2, 0 ] = attrArray[ 4 ]
- matrix[ 1, 2 ] = attrArray[ 3 ]
+ matrix[ 1, 2 ] = attrArray[ 3 ] # v4
matrix[ 2, 1 ] = attrArray[ 3 ]
+
+ elif attrArray.size == 9:
+ matrix[ 0, 1 ] = attrArray[ 5 ] # v1
+ matrix[ 1, 0 ] = attrArray[ 8 ] # v9
+
+ matrix[ 0, 2 ] = attrArray[ 4 ] # v5
+ matrix[ 2, 0 ] = attrArray[ 7 ] # v8
+
+ matrix[ 1, 2 ] = attrArray[ 3 ] # v4
+ matrix[ 2, 1 ] = attrArray[ 6 ] # v7
+
return matrix
@@ -67,37 +92,68 @@ def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], size: i
* 3x3 diagonal matrix:
.. math::
- \begin{bmatrix}
+ \\begin{bmatrix}
M00 & 0 & 0 \\
0 & M11 & 0 \\
0 & 0 & M22
- \end{bmatrix}
+ \\end{bmatrix}
=> (M00, M11, M22)
- * otherwise:
+ * 3x3 Generic matrix:
.. math::
- \begin{bmatrix}
+ \\begin{bmatrix}
+ M00 & M01 & M02 \\
+ M10 & M11 & M12 \\
+ M20 & M21 & M22
+ \\end{matrix}
+ => (M00, M11, M22, M12, M02, M01, M21, M20, M10)
+
+ * Symmetric case
+ .. math::
+
+ \\begin{bmatrix}
M00 & M01 & M02 \\
M01 & M11 & M12 \\
M02 & M12 & M22
- \end{bmatrix}
+ \\end{bmatrix}
=> (M00, M11, M22, M12, M02, M01)
+ .. Note::
+ In the case of 3 x 3 tensors, GEOS is currently only implemented for symmetrical cases.
+
Args:
attrMatrix (npt.NDArray[np.float64]): Matrix of attribute values.
- size (int): Size of the final vector.
+ size (int): Size of the final vector. Values accepted are 3, 9 or 6.
- Returns:
- npt.NDArray[np.float64]: vector of attribute values
+ Raises:
+ ValueError: The output vector size requested can only be 3, 9 or 6.
+ ValueError: Input matrix shape must be (3,3).
+ Returns:
+ npt.NDArray[np.float64]: Vector of attribute values.
"""
attrArray: npt.NDArray[ np.float64 ] = np.full( size, np.nan )
- # diagonal terms
+ if size not in ( 3, 6, 9 ):
+ raise ValueError( "Requested output size can only be 3, 9 or 6 (symmetrical case)." )
+ if attrMatrix.shape != ( 3, 3 ):
+ raise ValueError( "Input matrix shape must be (3,3)." )
+
+ # Diagonal terms
attrArray[ :3 ] = np.diag( attrMatrix )
+
# shear stress components
if attrArray.size == 6:
attrArray[ 3 ] = attrMatrix[ 1, 2 ]
attrArray[ 4 ] = attrMatrix[ 0, 2 ]
attrArray[ 5 ] = attrMatrix[ 0, 1 ]
+
+ elif attrArray.size == 9:
+ attrArray[ 3 ] = attrMatrix[ 1, 2 ]
+ attrArray[ 4 ] = attrMatrix[ 0, 2 ]
+ attrArray[ 5 ] = attrMatrix[ 0, 1 ]
+ attrArray[ 6 ] = attrMatrix[ 2, 1 ]
+ attrArray[ 7 ] = attrMatrix[ 2, 0 ]
+ attrArray[ 8 ] = attrMatrix[ 1, 0 ]
+
return attrArray
diff --git a/geos-utils/src/geos/utils/geometryFunctions.py b/geos-utils/src/geos/utils/geometryFunctions.py
index 6898b5e5..eb184375 100644
--- a/geos-utils/src/geos/utils/geometryFunctions.py
+++ b/geos-utils/src/geos/utils/geometryFunctions.py
@@ -9,6 +9,10 @@
__doc__ = """Functions to permform geometry calculations."""
+CANONICAL_BASIS_3D: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ],
+ npt.NDArray[ np.float64 ] ] = ( np.array( [ 1.0, 0.0, 0.0 ] ), np.array( [ 0.0, 1.0, 0.0 ] ),
+ np.array( [ 0.0, 0.0, 1.0 ] ) )
+
def getChangeOfBasisMatrix(
basisFrom: tuple[ npt.NDArray[ np.floating[ Any ] ], npt.NDArray[ np.floating[ Any ] ],
diff --git a/geos-utils/tests/test_algebraFunctions.py b/geos-utils/tests/test_algebraFunctions.py
new file mode 100644
index 00000000..104a8136
--- /dev/null
+++ b/geos-utils/tests/test_algebraFunctions.py
@@ -0,0 +1,73 @@
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Paloma Martinez
+
+import numpy as np
+import numpy.typing as npt
+from unittest import TestCase
+from typing_extensions import Self
+
+from geos.utils.algebraFunctions import (
+ getAttributeMatrixFromVector,
+ getAttributeVectorFromMatrix,
+)
+
+
+class TestAttributeMatrixFromVector( TestCase ):
+ """Tests for conversion from matrix to vector with GEOS convention."""
+
+ def test_wrongInputVectorSize( self: Self ) -> None:
+ """Test failure on incorrect input vector size."""
+ emptyVector: npt.NDArray[ np.float64 ] = np.array( [] )
+ with self.assertRaises( ValueError ):
+ getAttributeMatrixFromVector( emptyVector )
+
+ def test_vector3size( self: Self ) -> None:
+ """Test for an input vector size of 3."""
+ vector3 = np.arange( 1, 4 )
+ expectedMatrix: npt.NDArray[ np.float64 ] = np.array( [ [ 1, 0, 0 ], [ 0, 2, 0 ], [ 0, 0, 3 ] ] )
+
+ self.assertTrue( np.array_equal( expectedMatrix, getAttributeMatrixFromVector( vector3 ) ) )
+
+ def test_vector6( self: Self ) -> None:
+ """Test for an input vector size of 6."""
+ vector6 = np.arange( 1, 7 )
+ expectedMatrix = np.array( [ [ 1, 6, 5 ], [ 6, 2, 4 ], [ 5, 4, 3 ] ] )
+
+ self.assertTrue( np.array_equal( expectedMatrix, getAttributeMatrixFromVector( vector6 ) ) )
+
+ def test_vector9( self: Self ) -> None:
+ """Test for an input vector size of 9."""
+ vector9 = np.arange( 1, 10 )
+ expectedMatrix = np.array( [ [ 1, 6, 5 ], [ 9, 2, 4 ], [ 8, 7, 3 ] ] )
+
+ self.assertTrue( np.array_equal( expectedMatrix, getAttributeMatrixFromVector( vector9 ) ) )
+
+
+class TestAttributeVectorFromMatrix( TestCase ):
+ """Tests for conversion from vector to matrix with GEOS convention."""
+
+ def setUp( self: Self ) -> None:
+ """Set up parameters."""
+ self.rdMatrix = np.arange( 1, 10 ).reshape( 3, 3 )
+ self.expected: npt.NDArray[ np.float64 ] = np.array( [ 1, 5, 9, 6, 3, 2, 8, 7, 4 ] )
+
+ def test_wrongInputMatrixShape( self ) -> None:
+ """Test failure on empty input matrix."""
+ emptyMatrix = np.array( [] )
+ with self.assertRaises( ValueError ):
+ getAttributeVectorFromMatrix( emptyMatrix, 3 )
+
+ def test_wrongOutputSize( self: Self ) -> None:
+ """Test failure on incorrect output size requested."""
+ size = 4
+ with self.assertRaises( ValueError ):
+ getAttributeVectorFromMatrix( self.rdMatrix, size )
+
+ def test_vecOutput( self: Self ) -> None:
+ """Test correct output for requested size."""
+ listSize = ( 3, 6, 9 )
+ for size in listSize:
+ with self.subTest( size ):
+ expectedVector = np.array( self.expected[ :size ] )
+ self.assertTrue( np.array_equal( expectedVector, getAttributeVectorFromMatrix( self.rdMatrix, size ) ) )
diff --git a/geos-utils/tests/testsFunctionsGeosUtils.py b/geos-utils/tests/testsFunctionsGeosUtils.py
deleted file mode 100644
index 673553f8..00000000
--- a/geos-utils/tests/testsFunctionsGeosUtils.py
+++ /dev/null
@@ -1,26 +0,0 @@
-# SPDX-License-Identifier: Apache-2.0
-# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
-# SPDX-FileContributor: Alexandre Benedicto, Martin Lemay
-# ruff: noqa: E402 # disable Module level import not at top of file
-import unittest
-
-import geos.utils.algebraFunctions as fcts
-import numpy as np
-import numpy.typing as npt
-from typing_extensions import Self
-
-matrix: npt.NDArray[ np.float64 ] = np.array( [ [ 11, 21, 31 ], [ 21, 22, 23 ], [ 31, 23, 33 ] ] )
-vector: npt.NDArray[ np.float64 ] = np.array( [ 11, 22, 33, 23, 31, 21 ] )
-
-
-class TestsFunctionsalgebraFunctions( unittest.TestCase ):
-
- def test_getAttributeMatrixFromVector( self: Self ) -> None:
- """Test conversion from Matrix to Vector for Geos stress."""
- obtained: npt.NDArray[ np.float64 ] = fcts.getAttributeMatrixFromVector( vector )
- self.assertTrue( np.all( np.equal( obtained, matrix ) ) )
-
- def test_getAttributeVectorFromMatrix( self: Self ) -> None:
- """Test conversion from Vector to Matrix for Geos stress."""
- obtained: npt.NDArray[ np.float64 ] = fcts.getAttributeVectorFromMatrix( matrix, 6 )
- self.assertTrue( np.all( np.equal( obtained, vector ) ) )
diff --git a/install_packages.sh b/install_packages.sh
index e7160228..b6726ef7 100755
--- a/install_packages.sh
+++ b/install_packages.sh
@@ -2,8 +2,8 @@
python -m pip install --upgrade ./geos-utils
python -m pip install --upgrade ./geos-geomechanics
python -m pip install --upgrade ./geos-mesh
-python -m pip install --upgrade ./geos-posp
python -m pip install --upgrade ./geos-processing
+python -m pip install --upgrade ./geos-posp
python -m pip install --upgrade ./geos-xml-tools
python -m pip install --upgrade ./geos-xml-viewer
python -m pip install --upgrade ./hdf5-wrapper