From bb2911efe62dc31f9ef4771252731b96b9cf5f73 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Tue, 9 Sep 2025 14:36:04 +0200 Subject: [PATCH 01/23] Move SurfaceGeomechanicsFilter and refactor functions --- .../geos_posp/filters/SurfaceGeomechanics.py | 457 ----------------- geos-processing/pyproject.toml | 70 +++ .../post_processing/SurfaceGeomechanics.py | 473 ++++++++++++++++++ .../src/geos/utils/GeosOutputsConstants.py | 11 +- 4 files changed, 553 insertions(+), 458 deletions(-) delete mode 100644 geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py create mode 100644 geos-processing/pyproject.toml create mode 100644 geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py 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 new file mode 100644 index 00000000..3d947f67 --- /dev/null +++ b/geos-processing/pyproject.toml @@ -0,0 +1,70 @@ +[build-system] +requires = ["setuptools>=61.2", "wheel >= 0.37.1"] +build-backend = "setuptools.build_meta" + +[tool.setuptools] +include-package-data = true + +[tool.setuptools.packages.find] +where = ["src"] +include = ["geos.processing*"] +exclude = ['tests*'] + +[project] +name = "geos-processing" +version = "0.1.0" +description = "GEOS post- and pre-processing tools" +authors = [{name = "GEOS Contributors" }] +maintainers = [ + {name = "Alexandre Benedicto", email = "alexandre.benedicto@external.totalenergies.com" }, + {name = "Romain Baville", email = "romain.baville@external.totalenergies.com" }, + {name = "Paloma Martinez", email = "paloma.martinez@external.totalenergies.com" }, +] +license = {text = "Apache-2.0"} +classifiers = [ + "Development Status :: 4 - Beta", + "Programming Language :: Python" +] + +requires-python = ">=3.10" + +dependencies = [ + "geos-geomechanics", + "geos-utils", + "geos-mesh", + "vtk >= 9.3", + "numpy >= 2.2", + "typing_extensions >= 4.12", +] + + +[project.urls] +Homepage = "https://github.com/GEOS-DEV/geosPythonPackages" +Documentation = "https://geosx-geosx.readthedocs-hosted.com/projects/geosx-geospythonpackages/en/latest/" +Repository = "https://github.com/GEOS-DEV/geosPythonPackages.git" +"Bug Tracker" = "https://github.com/GEOS-DEV/geosPythonPackages/issues" + +[project.optional-dependencies] +build = [ + "build ~= 1.2" +] +dev = [ + "mypy", + "ruff", + "yapf", +] +test = [ + "pytest-cov", + "pytest" +] + +[tool.pytest.ini_options] +addopts = "--import-mode=importlib" +console_output_style = "count" +pythonpath = ["src"] +python_classes = "Test" +python_files = "test*.py" +python_functions = "test*" +testpaths = ["tests"] +norecursedirs = "bin" +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..f5e83be6 --- /dev/null +++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py @@ -0,0 +1,473 @@ +# 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 +from enum import Enum +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, + getAttributeToConvertFromXYZToLocal, +) +from geos.utils.Logger import logging, 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, +) +from vtkmodules.vtkCommonDataModel import ( + vtkPolyData, ) + +from geos.mesh.utils.arrayModifiers import createAttribute +from geos.mesh.utils.arrayHelpers import ( + getArrayInObject, + getAttributeSet, + isAttributeInObject, +) + + +__doc__ = """ +SurfaceGeomechanics vtk filter allows to compute Geomechanical +properties on surfaces. + +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 + filter: SurfaceGeomechanics = SurfaceGeomechanics( inputMesh, speHandler ) + + # Set the handler of yours (only if speHandler is True). + yourHandler: logging.Handler + filter.SetLoggerHandler( yourHandler ) + + # Set rock cohesion and friction angle (optional) + filter.SetRockCohesion( rockCohesion ) + filter.SetFrictionAngle( frictionAngle ) + + # Do calculations + filter.applyFilter() + + # Get output object + output: vtkPolyData = filter.GetOutputMesh() + + # Get created attribute names + newAttributeNames: set[str] = filter.GetNewAttributeNames() +""" +loggerTitle: str = "Surface Geomechanics" + +class SurfaceGeomechanics( VTKPythonAlgorithmBase ): + + 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. + """ + super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkPolyData" ) # type: ignore[no-untyped-call] + + # 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" surfacicMesh parameter is expected to be a vtkPolyData, not a {type(surfacicMesh)}." ) + self.inputMesh: vtkPolyData = surfacicMesh + # Output surfacic mesh + self.outputMesh: vtkPolyData + + # 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 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. Defaults to 10 / 180 * pi rad. + + Args: + frictionAngle (float): The friction angle in radians. + """ + self.frictionAngle = frictionAngle + + 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. + """ + self.logger.info( f"Applying filter { self.logger.name }." ) + + self.outputMesh = vtkPolyData() + self.outputMesh.ShallowCopy( self.inputMesh ) + + # Conversion of vectorial attributes from Normal/Tangent basis to xyz basis + if not self.convertLocalToXYZBasisAttributes(): + self.logger.error( "Error while converting Local to XYZ basis attributes." ) + # Compute shear capacity utilization + if not self.computeShearCapacityUtilization(): + self.logger.error( "Error while computing SCU." ) + + return True + + 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 + basisFrom = "local" + basisTo = "XYZ" + attributesToConvert: set[ str ] = self.__getAttributesToConvert( basisFrom, basisTo ) + if len( attributesToConvert ) == 0: + self.logger.error( f"No attribute to convert from {basisFrom} to {basisTo} basis were found." ) + return False + + # Get local coordinate vectors + normalTangentVectors: npt.NDArray[ np.float64 ] = self.getNormalTangentsVectors() + # Do conversion from local to canonic for each attribute + for attrName in attributesToConvert: + # Skip attribute if it is already in the object + newAttrName: str = attrName + "_" + ComponentNameEnum.XYZ.name + if isAttributeInObject( self.outputMesh, newAttrName, False ): + continue + + attr: vtkDoubleArray = self.outputMesh.GetCellData().GetArray( attrName ) + if attr is None: + self.logger.error( f"Attribute {attrName} is undefined." ) + if not attr.GetNumberOfComponents() > 2: + self.logger.error( f"Dimension of the attribute {attrName} must be equal or greater 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 + if not createAttribute( + self.outputMesh, + newAttrArray, + newAttrName, + ComponentNameEnum.XYZ.value, + False, + logger = self.logger + ): + self.logger.error( f"Failed to create attribute {newAttrName}." ) + else: + self.logger.info( f"Attribute {newAttrName} added to the output mesh." ) + self.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 + basisFrom = "canonical" + basisTo = "local" + attributesToConvert: set[ str ] = self.__getAttributesToConvert( basisFrom, basisTo ) + if len( attributesToConvert ) == 0: + self.logger.error( "No attribute to convert from canonic to local basis were found." ) + return False + + # get local coordinate vectors + normalTangentVectors: npt.NDArray[ np.float64 ] = self.getNormalTangentsVectors() + # Do conversion from canonic to local for each attribute + for attrName in attributesToConvert: + # Skip attribute if it is already in the object + newAttrName: str = attrName + "_" + ComponentNameEnum.NORMAL_TANGENTS.name + if isAttributeInObject( self.outputMesh, newAttrName, False ): + continue + + attr: vtkDoubleArray = self.outputMesh.GetCellData().GetArray( attrName ) + if attr is None: + self.logger.error( f"Attribute {attrName} is undefined." ) + if not attr.GetNumberOfComponents() > 2: + self.logger.error( f"Dimension of the attribute {attrName} must be equal or greater 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 + if not createAttribute( + self.outputMesh, + newAttrArray, + newAttrName, + ComponentNameEnum.NORMAL_TANGENTS.value, + False, + logger = self.logger + ): + self.logger.error( f"Failed to create attribute {newAttrName}." ) + return False + else: + self.logger.info( f"Attribute {newAttrName} added to the output mesh." ) + self.newAttributeNames.add( newAttrName ) + return True + + + def __getAttributesToConvert( self: Self, basisFrom: str, basisTo: str ): + """Get the list of attributes to convert from one basis to another. + Choices of basis are "local" or "XYZ". + + Returns: + set[ str ]: Set of the attributes names + """ + if basisFrom == "local" and basisTo == "XYZ": + return getAttributeToConvertFromLocalToXYZ() + elif basisFrom == "XYZ" and basisTo == "local": + return getAttributeToConvertFromXYZToLocal() + else: + self.logger.error( f"Unknow basis named {basisFrom} or {basisTo}." ) + + + def filterAttributesToConvert( self: Self, attributesToFilter: set[ str ] ) -> set[ str ]: + """Filter the set of attribute names if they are vectorial and present. + + Args: + attributesToFilter (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.outputMesh, False ) + for attrName in attributesToFilter: + if attrName in attributeSet: + attr: vtkDataArray = self.outputMesh.GetCellData().GetArray( attrName ) + if attr.GetNumberOfComponents() > 2: + attributesFiltered.add( attrName ) + 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 ) + + if not np.any( np.isfinite( attrArrayNew ) ): + self.logger.error( "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.outputMesh.GetCellData().GetNormals() ) # type: ignore[no-untyped-call] + if normals is None: + self.logger.error( "Normal attribute was not found." ) + tangents1: npt.NDArray[ np.float64 ] = vnp.vtk_to_numpy( + self.outputMesh.GetCellData().GetTangents() ) # type: ignore[no-untyped-call] + if tangents1 is None: + self.logger.error( "Tangents attribute was not found." ) + + # Compute second tangential component + tangents2: npt.NDArray[ np.float64 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 ) + if tangents2 is None: + self.logger.error( "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 (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 ): + tractionAttributeName: str = GeosMeshOutputsEnum.TRACTION.attributeName + traction: npt.NDArray[ np.float64 ] = getArrayInObject( self.outputMesh, tractionAttributeName, + self.attributeOnPoints ) + if traction is None: + self.logger.error( f"{tractionAttributeName} attribute is undefined." ) + self.logger.error( f"Failed to compute {SCUAttributeName}." ) + + scuAttribute: npt.NDArray[ np.float64 ] = fcts.shearCapacityUtilization( + traction, self.rockCohesion, self.frictionAngle ) + + # 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( f"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-utils/src/geos/utils/GeosOutputsConstants.py b/geos-utils/src/geos/utils/GeosOutputsConstants.py index 3e684aaa..c98b3fe8 100644 --- a/geos-utils/src/geos/utils/GeosOutputsConstants.py +++ b/geos-utils/src/geos/utils/GeosOutputsConstants.py @@ -306,8 +306,17 @@ 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 + list[str]: set of attributes to convert """ return { GeosMeshOutputsEnum.DISPLACEMENT_JUMP.attributeName, } + +def getAttributeToConvertFromXYZToLocal() -> set[ str ]: + """Get the list of attribute names to convert from canonical (XYZ) to local basis. + This list is empty at the moment. + + Returns: + set[ str ]: set of attributes to convert + """ + return set() \ No newline at end of file From 391aeea9d3671118a7f21189d64b05f3ddc784a1 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Tue, 9 Sep 2025 16:49:14 +0200 Subject: [PATCH 02/23] Move and refactor plugin --- .../src/PVplugins/PVSurfaceGeomechanics.py | 250 ------------------ .../post_processing/SurfaceGeomechanics.py | 6 + geos-pv/requirements.txt | 3 +- .../geos/pv/plugins/PVSurfaceGeomechanics.py | 197 ++++++++++++++ 4 files changed, 205 insertions(+), 251 deletions(-) delete mode 100644 geos-posp/src/PVplugins/PVSurfaceGeomechanics.py create mode 100644 geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py 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-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py index f5e83be6..30d3d6ec 100644 --- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py +++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py @@ -109,6 +109,8 @@ def __init__( self: Self, self.logger.setLevel( logging.INFO ) # Input surfacic mesh + if surfacicMesh is None: + self.logger.error( "Input surface is undefined." ) if not surfacicMesh.IsA( "vtkPolyData" ): self.logger.error( f" surfacicMesh parameter is expected to be a vtkPolyData, not a {type(surfacicMesh)}." ) self.inputMesh: vtkPolyData = surfacicMesh @@ -178,10 +180,14 @@ def applyFilter( self: Self ) -> bool: # Conversion of vectorial attributes from Normal/Tangent basis to xyz basis if not self.convertLocalToXYZBasisAttributes(): self.logger.error( "Error while converting Local to XYZ basis attributes." ) + 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} succeeded." ) return True def convertLocalToXYZBasisAttributes( self: Self ) -> bool: diff --git a/geos-pv/requirements.txt b/geos-pv/requirements.txt index edb1046c..6ec1b5a9 100644 --- a/geos-pv/requirements.txt +++ b/geos-pv/requirements.txt @@ -1,4 +1,5 @@ geos-geomechanics geos-mesh geos-posp -geos-utils \ No newline at end of file +geos-utils +geos-processing \ No newline at end of file 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..756f2afb --- /dev/null +++ b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py @@ -0,0 +1,197 @@ +# 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, smhint, smproperty, smproxy, +) +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 + +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 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 ( + 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. + + +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. + +""" + +@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: + """Compute additional geomechanical surface outputs. + + Input is a vtkMultiBlockDataSet that contains surfaces with + Normals and Tangential attributes. + """ + super().__init__( + nInputPorts=1, + nOutputPorts=1, + inputType="vtkMultiBlockDataSet", + outputType="vtkMultiBlockDataSet", + ) + # 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 frition angle. + + Args: + value (float): friction angle (°) + """ + self.frictionAngle = value + self.Modified() + + 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. + """ + inputMesh: vtkMultiBlockDataSet = vtkMultiBlockDataSet.GetData( inInfoVec[ 0 ] ) + output: vtkMultiBlockDataSet = self.GetOutputData( outInfoVec, 0 ) + + assert inputMesh is not None, "Input surface is null." + assert output is not None, "Output pipeline is null." + + output.ShallowCopy( inputMesh ) + + surfaceBlockIndexes: list[ int ] = getBlockElementIndexesFlatten( inputMesh ) + for blockIndex in surfaceBlockIndexes: + surfaceBlock: vtkPolyData = vtkPolyData.SafeDownCast( getBlockFromFlatIndex( output, blockIndex ) ) + + filter: SurfaceGeomechanics = SurfaceGeomechanics( surfaceBlock, True ) + + if not filter.logger.hasHandlers(): + filter.SetLoggerHandler( VTKHandler() ) + + filter.SetRockCohesion( self._getRockCohesion() ) + filter.SetFrictionAngle( self._getFrictionAngle() ) + filter.applyFilter() + + outputSurface: vtkPolyData = filter.GetOutputMesh() + + # 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() + + # TODO : modify logger + return 1 + + 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 \ No newline at end of file From 176b508791a494d8facee083a01b4f1842364dcb Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Wed, 10 Sep 2025 14:48:47 +0200 Subject: [PATCH 03/23] add geos-processing to install script --- .../src/geos/processing/post_processing/SurfaceGeomechanics.py | 2 +- install_packages.sh | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py index 30d3d6ec..e9856bec 100644 --- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py +++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py @@ -112,7 +112,7 @@ def __init__( self: Self, if surfacicMesh is None: self.logger.error( "Input surface is undefined." ) if not surfacicMesh.IsA( "vtkPolyData" ): - self.logger.error( f" surfacicMesh parameter is expected to be a vtkPolyData, not a {type(surfacicMesh)}." ) + self.logger.error( f"Input surface is expected to be a vtkPolyData, not a {type(surfacicMesh)}." ) self.inputMesh: vtkPolyData = surfacicMesh # Output surfacic mesh self.outputMesh: vtkPolyData diff --git a/install_packages.sh b/install_packages.sh index 74238daf..e7160228 100755 --- a/install_packages.sh +++ b/install_packages.sh @@ -3,6 +3,7 @@ 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-xml-tools python -m pip install --upgrade ./geos-xml-viewer python -m pip install --upgrade ./hdf5-wrapper From 8b72f3f79b4ed1ef76059d1c7d731010c5a61fcc Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Fri, 12 Sep 2025 11:31:51 +0200 Subject: [PATCH 04/23] Add more log --- .../post_processing/SurfaceGeomechanics.py | 25 ++++++++++++++++--- .../geos/pv/plugins/PVSurfaceGeomechanics.py | 7 ++---- 2 files changed, 23 insertions(+), 9 deletions(-) diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py index e9856bec..3f2f9f43 100644 --- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py +++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py @@ -2,7 +2,7 @@ # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. # SPDX-FileContributor: Martin Lemay # ruff: noqa: E402 # disable Module level import not at top of file -from enum import Enum +from typing_extensions import Self import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts import geos.utils.geometryFunctions as geom import numpy as np @@ -24,7 +24,7 @@ DEFAULT_FRICTION_ANGLE_RAD, DEFAULT_ROCK_COHESION, ) -from typing_extensions import Self + from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase from vtkmodules.vtkCommonCore import ( vtkDataArray, @@ -114,6 +114,8 @@ def __init__( self: Self, 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 = None # Output surfacic mesh self.outputMesh: vtkPolyData @@ -142,6 +144,14 @@ def SetLoggerHandler( self: Self, handler: Logger ) -> None: "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 ): + """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. @@ -158,6 +168,7 @@ def SetFrictionAngle( self: Self, frictionAngle: float ) -> None: """ self.frictionAngle = frictionAngle + def GetNewAttributeNames( self: Self ) -> set[ str ]: """Get the set of new attribute names that were created. @@ -172,7 +183,13 @@ def applyFilter( self: Self ) -> bool: Returns: int: 1 if calculation successfully ended, 0 otherwise. """ - self.logger.info( f"Applying filter { self.logger.name }." ) + 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 ) @@ -242,7 +259,7 @@ def convertXYZToLocalBasisAttributes( self: Self ) -> bool: """Convert vectorial property coordinates from canonic to local basis. Returns: - bool: True if calculation successfully ended, False otherwise + bool: True if calculation successfully ended, False otherwise. """ # Look for the list of attributes to convert basisFrom = "canonical" diff --git a/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py index 756f2afb..ff525676 100644 --- a/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py +++ b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py @@ -14,7 +14,6 @@ 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" ) ) @@ -158,7 +157,7 @@ def RequestData( surfaceBlock: vtkPolyData = vtkPolyData.SafeDownCast( getBlockFromFlatIndex( output, blockIndex ) ) filter: SurfaceGeomechanics = SurfaceGeomechanics( surfaceBlock, True ) - + filter.SetSurfaceName( f"blockIndex {blockIndex}" ) if not filter.logger.hasHandlers(): filter.SetLoggerHandler( VTKHandler() ) @@ -176,8 +175,6 @@ def RequestData( surfaceBlock.Modified() output.Modified() - - # TODO : modify logger return 1 def _getFrictionAngle( self: Self ) -> float: @@ -194,4 +191,4 @@ def _getRockCohesion( self: Self ) -> float: Returns: float: rock cohesion. """ - return self.rockCohesion \ No newline at end of file + return self.rockCohesion From 5bd267c78d054757820e2dd9118d92a0358db960 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Thu, 18 Sep 2025 16:46:25 +0200 Subject: [PATCH 05/23] Add documentation --- docs/geos-processing.rst | 14 ++++++++++++++ .../generic_processing_tool.rst | 4 ++++ docs/geos_processing_docs/post_processing.rst | 15 +++++++++++++++ docs/geos_processing_docs/pre_processing.rst | 4 ++++ docs/index.rst | 6 ++++-- 5 files changed, 41 insertions(+), 2 deletions(-) create mode 100644 docs/geos-processing.rst create mode 100644 docs/geos_processing_docs/generic_processing_tool.rst create mode 100644 docs/geos_processing_docs/post_processing.rst create mode 100644 docs/geos_processing_docs/pre_processing.rst diff --git a/docs/geos-processing.rst b/docs/geos-processing.rst new file mode 100644 index 00000000..80f56429 --- /dev/null +++ b/docs/geos-processing.rst @@ -0,0 +1,14 @@ +GEOS Processing tools +======================== + + +.. toctree:: + :maxdepth: 5 + :caption: Contents: + + ./geos_processing_docs/post_processing.rst + + ./geos_processing_docs/pre_processing.rst + + ./geos_processing_docs/generic_processing_tools.rst + diff --git a/docs/geos_processing_docs/generic_processing_tool.rst b/docs/geos_processing_docs/generic_processing_tool.rst new file mode 100644 index 00000000..4fdb32fc --- /dev/null +++ b/docs/geos_processing_docs/generic_processing_tool.rst @@ -0,0 +1,4 @@ +Generic processing filters +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In progress.. \ No newline at end of file diff --git a/docs/geos_processing_docs/post_processing.rst b/docs/geos_processing_docs/post_processing.rst new file mode 100644 index 00000000..f23e4fea --- /dev/null +++ b/docs/geos_processing_docs/post_processing.rst @@ -0,0 +1,15 @@ +Post-processing filters +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This module contains GEOS post-processing tools. + +Geomechanics post-processing tools +===================================== + +geos.processing.post_processing.SurfaceGeomechanics +------------------------------------------------------- + +.. automodule:: geos.processing.post_processing.SurfaceGeomechanics + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/geos_processing_docs/pre_processing.rst b/docs/geos_processing_docs/pre_processing.rst new file mode 100644 index 00000000..85973093 --- /dev/null +++ b/docs/geos_processing_docs/pre_processing.rst @@ -0,0 +1,4 @@ +Pre-processing filters +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In progress.. \ No newline at end of file diff --git a/docs/index.rst b/docs/index.rst index 536616d1..034d0e45 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -47,7 +47,7 @@ To do this, you can clone a copy of the geosPythonPackages repository and instal cd geosPythonPackages/ python -m pip install --upgrade geos-ats - + .. note:: To upgrade an existing installation, the python executable in the above command should correspond to the version you indicated in your host config. If you have previously built the tools, this version will be linked in the build directory: `build_dir/bin/python`. @@ -87,7 +87,9 @@ Packages geos-mesh geos-posp - + + geos-processing + geos-pv geos-timehistory From 4277ebb1348005d2ea8e357893887d2e44c469e1 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Fri, 3 Oct 2025 11:15:54 +0200 Subject: [PATCH 06/23] Filter cleaning --- docs/geos_processing_docs/post_processing.rst | 2 +- .../post_processing/SurfaceGeomechanics.py | 281 +++++++++--------- .../geos/pv/plugins/PVSurfaceGeomechanics.py | 5 +- .../src/geos/utils/GeosOutputsConstants.py | 22 +- 4 files changed, 152 insertions(+), 158 deletions(-) diff --git a/docs/geos_processing_docs/post_processing.rst b/docs/geos_processing_docs/post_processing.rst index f23e4fea..45ea4857 100644 --- a/docs/geos_processing_docs/post_processing.rst +++ b/docs/geos_processing_docs/post_processing.rst @@ -12,4 +12,4 @@ geos.processing.post_processing.SurfaceGeomechanics .. automodule:: geos.processing.post_processing.SurfaceGeomechanics :members: :undoc-members: - :show-inheritance: \ No newline at end of file + :show-inheritance: diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py index 3f2f9f43..0662ce24 100644 --- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py +++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py @@ -2,29 +2,13 @@ # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. # SPDX-FileContributor: Martin Lemay # ruff: noqa: E402 # disable Module level import not at top of file -from typing_extensions import Self -import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts -import geos.utils.geometryFunctions as geom +import logging import numpy as np + +from typing_extensions import Self, Union 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, - getAttributeToConvertFromXYZToLocal, -) -from geos.utils.Logger import logging, Logger, getLogger -from geos.utils.PhysicalConstants import ( - DEFAULT_FRICTION_ANGLE_RAD, - DEFAULT_ROCK_COHESION, -) +import vtkmodules.util.numpy_support as vnp from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase from vtkmodules.vtkCommonCore import ( vtkDataArray, @@ -39,11 +23,35 @@ getAttributeSet, isAttributeInObject, ) - +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.algebraFunctions import ( + getAttributeMatrixFromVector, + getAttributeVectorFromMatrix, +) +from geos.utils.GeosOutputsConstants import ( + ComponentNameEnum, + GeosMeshOutputsEnum, + PostProcessingOutputsEnum, +) +import geos.utils.geometryFunctions as geom __doc__ = """ -SurfaceGeomechanics vtk filter allows to compute Geomechanical -properties on surfaces. +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. + +.. 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. @@ -68,7 +76,7 @@ yourHandler: logging.Handler filter.SetLoggerHandler( yourHandler ) - # Set rock cohesion and friction angle (optional) + # [optional] Set rock cohesion and friction angle filter.SetRockCohesion( rockCohesion ) filter.SetFrictionAngle( frictionAngle ) @@ -80,6 +88,14 @@ # Get created attribute names newAttributeNames: set[str] = filter.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" @@ -118,6 +134,9 @@ def __init__( self: Self, self.name = 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 @@ -161,7 +180,7 @@ def SetRockCohesion( self: Self, rockCohesion: float ) -> None: self.rockCohesion = rockCohesion def SetFrictionAngle( self: Self, frictionAngle: float ) -> None: - """Set friction angle value. Defaults to 10 / 180 * pi rad. + """Set friction angle value in radians. Defaults to 10 / 180 * pi rad. Args: frictionAngle (float): The friction angle in radians. @@ -169,6 +188,58 @@ def SetFrictionAngle( self: Self, frictionAngle: float ) -> None: 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.add( attributeName ) + self.ConvertAttributesOn() + + def GetNewAttributeNames( self: Self ) -> set[ str ]: """Get the set of new attribute names that were created. @@ -177,6 +248,7 @@ def GetNewAttributeNames( self: Self ) -> set[ str ]: """ return self.newAttributeNames + def applyFilter( self: Self ) -> bool: """Compute Geomechanical properties on input surface. @@ -194,10 +266,12 @@ def applyFilter( self: Self ) -> bool: self.outputMesh = vtkPolyData() self.outputMesh.ShallowCopy( self.inputMesh ) - # Conversion of vectorial attributes from Normal/Tangent basis to xyz basis - if not self.convertLocalToXYZBasisAttributes(): - self.logger.error( "Error while converting Local to XYZ basis attributes." ) - return False + # 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(): @@ -207,138 +281,75 @@ def applyFilter( self: Self ) -> bool: self.logger.info( f"Filter {self.logger.name} succeeded." ) return True - def convertLocalToXYZBasisAttributes( self: Self ) -> bool: - """Convert vectorial property coordinates from local to canonic basis. + + def _convertAttributesFromLocalToXYZBasis( self: Self ) -> bool: + """Convert attributes from local to XYZ basis. Returns: - bool: True if calculation successfully ended, False otherwise + bool: True if calculation successfully ended or no attributes, False otherwise """ - # Look for the list of attributes to convert - basisFrom = "local" - basisTo = "XYZ" - attributesToConvert: set[ str ] = self.__getAttributesToConvert( basisFrom, basisTo ) + # Get the list of attributes to convert and filter + attributesToConvert: set[ str ] = self.__filterAttributesToConvert( self.attributesToConvert ) + if len( attributesToConvert ) == 0: - self.logger.error( f"No attribute to convert from {basisFrom} to {basisTo} basis were found." ) - return False + self.logger.warning( f"No attribute to convert from local to XYZ basis were found." ) + return True # Get local coordinate vectors - normalTangentVectors: npt.NDArray[ np.float64 ] = self.getNormalTangentsVectors() - # Do conversion from local to canonic for each attribute - for attrName in attributesToConvert: + normalTangentVectors: npt.NDArray[ np.float64 ] = self.__getNormalTangentsVectors() + # Do conversion from local to XYZ for each attribute + for attrNameLocal in attributesToConvert: # Skip attribute if it is already in the object - newAttrName: str = attrName + "_" + ComponentNameEnum.XYZ.name - if isAttributeInObject( self.outputMesh, newAttrName, False ): + attrNameXYZ: str = f" {attrNameLocal}_{ComponentNameEnum.XYZ.name}" + if isAttributeInObject( self.outputMesh, attrNameXYZ, self.attributeOnPoints ): continue - attr: vtkDoubleArray = self.outputMesh.GetCellData().GetArray( attrName ) - if attr is None: - self.logger.error( f"Attribute {attrName} is undefined." ) - if not attr.GetNumberOfComponents() > 2: - self.logger.error( f"Dimension of the attribute {attrName} must be equal or greater than 3." ) + attrArray: vtkDoubleArray = getArrayInObject( self.outputMesh, attrNameLocal, self.attributeOnPoints) - attrArray: npt.NDArray[ np.float64 ] = vnp.vtk_to_numpy( attr ) # type: ignore[no-untyped-call] - newAttrArray: npt.NDArray[ np.float64 ] = self.computeNewCoordinates( attrArray, normalTangentVectors, + newAttrArray: npt.NDArray[ np.float64 ] = self.__computeNewCoordinates( attrArray, normalTangentVectors, True ) # create attribute - if not createAttribute( + if createAttribute( self.outputMesh, newAttrArray, - newAttrName, + attrNameXYZ, ComponentNameEnum.XYZ.value, - False, - logger = self.logger - ): - self.logger.error( f"Failed to create attribute {newAttrName}." ) - else: - self.logger.info( f"Attribute {newAttrName} added to the output mesh." ) - self.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 - basisFrom = "canonical" - basisTo = "local" - attributesToConvert: set[ str ] = self.__getAttributesToConvert( basisFrom, basisTo ) - if len( attributesToConvert ) == 0: - self.logger.error( "No attribute to convert from canonic to local basis were found." ) - return False - - # get local coordinate vectors - normalTangentVectors: npt.NDArray[ np.float64 ] = self.getNormalTangentsVectors() - # Do conversion from canonic to local for each attribute - for attrName in attributesToConvert: - # Skip attribute if it is already in the object - newAttrName: str = attrName + "_" + ComponentNameEnum.NORMAL_TANGENTS.name - if isAttributeInObject( self.outputMesh, newAttrName, False ): - continue + onPoints = self.attributeOnPoints, + logger = self.logger): + self.logger.info( f"Attribute {attrNameXYZ} added to the output mesh." ) + self.newAttributeNames.add( attrNameXYZ ) - attr: vtkDoubleArray = self.outputMesh.GetCellData().GetArray( attrName ) - if attr is None: - self.logger.error( f"Attribute {attrName} is undefined." ) - if not attr.GetNumberOfComponents() > 2: - self.logger.error( f"Dimension of the attribute {attrName} must be equal or greater 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 - if not createAttribute( - self.outputMesh, - newAttrArray, - newAttrName, - ComponentNameEnum.NORMAL_TANGENTS.value, - False, - logger = self.logger - ): - self.logger.error( f"Failed to create attribute {newAttrName}." ) - return False - else: - self.logger.info( f"Attribute {newAttrName} added to the output mesh." ) - self.newAttributeNames.add( newAttrName ) return True - def __getAttributesToConvert( self: Self, basisFrom: str, basisTo: str ): - """Get the list of attributes to convert from one basis to another. - Choices of basis are "local" or "XYZ". - - Returns: - set[ str ]: Set of the attributes names - """ - if basisFrom == "local" and basisTo == "XYZ": - return getAttributeToConvertFromLocalToXYZ() - elif basisFrom == "XYZ" and basisTo == "local": - return getAttributeToConvertFromXYZToLocal() - else: - self.logger.error( f"Unknow basis named {basisFrom} or {basisTo}." ) - - - def filterAttributesToConvert( self: Self, attributesToFilter: set[ str ] ) -> set[ str ]: + def __filterAttributesToConvert( self: Self ) -> set[ str ]: """Filter the set of attribute names if they are vectorial and present. - Args: - attributesToFilter (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.outputMesh, False ) - for attrName in attributesToFilter: - if attrName in attributeSet: - attr: vtkDataArray = self.outputMesh.GetCellData().GetArray( attrName ) - if attr.GetNumberOfComponents() > 2: - attributesFiltered.add( attrName ) + + 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 - def computeNewCoordinates( + def __computeNewCoordinates( self: Self, attrArray: npt.NDArray[ np.float64 ], normalTangentVectors: npt.NDArray[ np.float64 ], @@ -427,7 +438,7 @@ def __computeChangeOfBasisMatrix( self: Self, localBasis: npt.NDArray[ np.float6 # Inverse the change of basis matrix return np.linalg.inv( P ).astype( np.float64 ) - def getNormalTangentsVectors( self: Self ) -> npt.NDArray[ np.float64 ]: + def __getNormalTangentsVectors( self: Self ) -> npt.NDArray[ np.float64 ]: """Compute the change of basis matrix from Local to XYZ bases. Returns: @@ -460,6 +471,7 @@ def computeShearCapacityUtilization( self: Self ) -> bool: 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 ) @@ -467,6 +479,7 @@ def computeShearCapacityUtilization( self: Self ) -> bool: self.logger.error( f"{tractionAttributeName} attribute is undefined." ) self.logger.error( f"Failed to compute {SCUAttributeName}." ) + # Computation of the shear capacity utilization (SCU) scuAttribute: npt.NDArray[ np.float64 ] = fcts.shearCapacityUtilization( traction, self.rockCohesion, self.frictionAngle ) diff --git a/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py index ff525676..ff198ba2 100644 --- a/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py +++ b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py @@ -49,7 +49,8 @@ Input and output are vtkMultiBlockDataSet. .. Important:: - Please refer to the GeosExtractMergeBlockVolumeSurface* filters to provide the correct input. + - 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: @@ -57,7 +58,7 @@ * 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. +* Select Filters > 3- Geos Geomechanics > Geos Surface Geomechanics. * (Optional) Set rock cohesion and/or friction angle. * Apply. diff --git a/geos-utils/src/geos/utils/GeosOutputsConstants.py b/geos-utils/src/geos/utils/GeosOutputsConstants.py index c98b3fe8..1ad99fbf 100644 --- a/geos-utils/src/geos/utils/GeosOutputsConstants.py +++ b/geos-utils/src/geos/utils/GeosOutputsConstants.py @@ -299,24 +299,4 @@ def getAttributeToTransferFromInitialTime() -> dict[ str, str ]: PostProcessingOutputsEnum.YOUNG_MODULUS_INITIAL.attributeName, 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]: set of attributes to convert - """ - return { - GeosMeshOutputsEnum.DISPLACEMENT_JUMP.attributeName, - } - -def getAttributeToConvertFromXYZToLocal() -> set[ str ]: - """Get the list of attribute names to convert from canonical (XYZ) to local basis. - This list is empty at the moment. - - Returns: - set[ str ]: set of attributes to convert - """ - return set() \ No newline at end of file + } \ No newline at end of file From 2e715381b0dcfde2c5ba2a05c9a3c95e947e9493 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Fri, 3 Oct 2025 11:36:52 +0200 Subject: [PATCH 07/23] Documentation --- docs/geos_processing_docs/post_processing.rst | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/docs/geos_processing_docs/post_processing.rst b/docs/geos_processing_docs/post_processing.rst index 45ea4857..839f4a7e 100644 --- a/docs/geos_processing_docs/post_processing.rst +++ b/docs/geos_processing_docs/post_processing.rst @@ -6,6 +6,17 @@ 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: + - rock cohesion: 0.0 Pa $( fractured case ) + - friction angle: 10° + + geos.processing.post_processing.SurfaceGeomechanics ------------------------------------------------------- From b6f9d859f0b713c469ce3f52a5801ab51d8c5615 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Wed, 15 Oct 2025 16:50:04 +0200 Subject: [PATCH 08/23] Refactor the basis change computation --- .../src/geos/mesh/utils/genericHelpers.py | 83 +++++++++- .../post_processing/SurfaceGeomechanics.py | 148 ++++-------------- geos-utils/src/geos/utils/algebraFunctions.py | 98 ++++++++++-- 3 files changed, 201 insertions(+), 128 deletions(-) diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py index 481b391f..57ee6b22 100644 --- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py +++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py @@ -4,12 +4,22 @@ import numpy as np import numpy.typing as npt from typing import Iterator, List, Sequence, Any, Union -from vtkmodules.util.numpy_support import numpy_to_vtk +from vtkmodules.util.numpy_support import ( + numpy_to_vtk, + vtk_to_numpy +) 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 geos.mesh.utils.multiblockHelpers import ( getBlockElementIndexesFlatten, getBlockFromFlatIndex ) +from geos.utils.algebraFunctions import ( + getAttributeMatrixFromVector, + getAttributeVectorFromMatrix, + getLocalToXYZTransformMatrix, +) + + __doc__ = """ Generic VTK utilities. @@ -267,3 +277,74 @@ def createVertices( cellPtsCoord: list[ npt.NDArray[ np.float64 ] ], cellVertexMap += [ ptId.get() ] # type: ignore cellVertexMapAll += [ tuple( cellVertexMap ) ] return points, cellVertexMapAll + + +def convertAttributeFromLocalToXYZForOneCell( vector, localBasisVectors ) -> 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 (np.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 ] = getLocalToXYZTransformMatrix( localBasisVectors ) + + # Apply transformation + arrayXYZ: npt.NDArray[ np.float64 ] = transformMatrix @ matrix3x3 @ transformMatrix.T + + # Convert back to GEOS type attribute + vectorXYZ: npt.NDArray[ np.float64 ] = getAttributeVectorFromMatrix( arrayXYZ, vector.size ) + + return vectorXYZ + + +def getNormalVectors( dataset, ) -> npt.NDArray[ np.float64 ]: + normals: npt.NDArray[ np.float64 ] = vtk_to_numpy( + dataset.GetCellData().GetNormals() ) # type: ignore[no-untyped-call] + # TODO: if normals not defined, compute it on the fly + if normals is None: + raise ValueError( "Normal attribute was not found." ) + + return normals + + +def getTangentsVectors( dataset, normals = None ) -> npt.NDArray[ np.float64 ]: + # Get first tangential component + tangents1: npt.NDArray[ np.float64 ] = vtk_to_numpy( + dataset.GetCellData().GetTangents() ) # type: ignore[no-untyped-call] + if tangents1 is None: + raise ValueError( "Tangents attribute was not found." ) + + # Compute second tangential component + if normals is None: + normals = getNormalVectors( dataset ) + + tangents2: npt.NDArray[ np.float64 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 ) + if tangents2 is None: + raise ValueError( "Local basis third axis was not computed." ) + + return ( tangents1, tangents2 ) + + +def getLocalBasisVectors( dataset: vtkPolyData ) -> npt.NDArray[ np.float64 ]: + """Return the local basis vectors for all cells of the input dataset. + + Args: + dataset(vtkPolydata): The input dataset + + """ + normals = getNormalVectors( dataset ) + tangents = getTangentsVectors( dataset, normals=normals ) + + return np.array( ( normals, *tangents ) ) diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py index 0662ce24..c5dba4d9 100644 --- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py +++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py @@ -23,16 +23,16 @@ 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.algebraFunctions import ( - getAttributeMatrixFromVector, - getAttributeVectorFromMatrix, -) from geos.utils.GeosOutputsConstants import ( ComponentNameEnum, GeosMeshOutputsEnum, @@ -46,7 +46,8 @@ - Computation of the shear capacity utilization (SCU) .. Warning:: - The computation of the SCU requires the presence of a 'traction' attribute in the input mesh. + - 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: @@ -269,7 +270,7 @@ def applyFilter( self: Self ) -> bool: # 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(): + if not self.convertAttributesFromLocalToXYZBasis(): self.logger.error( "Error while converting attributes from local to XYZ basis." ) return False @@ -278,41 +279,41 @@ def applyFilter( self: Self ) -> bool: self.logger.error( "Error while computing SCU." ) return False - self.logger.info( f"Filter {self.logger.name} succeeded." ) + self.logger.info( f"Filter {self.logger.name} successfully applied on surface {self.name}." ) return True - def _convertAttributesFromLocalToXYZBasis( self: Self ) -> bool: + def convertAttributesFromLocalToXYZBasis( self: Self ) -> bool: """Convert attributes from local to XYZ basis. Returns: - bool: True if calculation successfully ended or no attributes, False otherwise + 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( self.attributesToConvert ) + attributesToConvert: set[ str ] = self.__filterAttributesToConvert() if len( attributesToConvert ) == 0: self.logger.warning( f"No attribute to convert from local to XYZ basis were found." ) return True - # Get local coordinate vectors - normalTangentVectors: npt.NDArray[ np.float64 ] = self.__getNormalTangentsVectors() # 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 - attrNameXYZ: str = f" {attrNameLocal}_{ComponentNameEnum.XYZ.name}" if isAttributeInObject( self.outputMesh, attrNameXYZ, self.attributeOnPoints ): continue - attrArray: vtkDoubleArray = getArrayInObject( self.outputMesh, attrNameLocal, self.attributeOnPoints) + if self.attributeOnPoints: + self.logger.error( "This filter can only convert cell attributes from local to XYZ basis, not point attributes." ) + localArray: vtkDoubleArray = getArrayInObject( self.outputMesh, attrNameLocal, self.attributeOnPoints) - newAttrArray: npt.NDArray[ np.float64 ] = self.__computeNewCoordinates( attrArray, normalTangentVectors, - True ) + arrayXYZ: npt.NDArray[ np.float64 ] = self.__computeXYZCoordinates( localArray ) - # create attribute + # Create converted attribute array in dataset if createAttribute( self.outputMesh, - newAttrArray, + arrayXYZ, attrNameXYZ, ComponentNameEnum.XYZ.value, onPoints = self.attributeOnPoints, @@ -349,118 +350,37 @@ def __filterAttributesToConvert( self: Self ) -> set[ str ]: return attributesFiltered - def __computeNewCoordinates( + # TODO: Adapt to handle point attributes. + def __computeXYZCoordinates( 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. + """Compute the XYZ 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 ) - - if not np.any( np.isfinite( attrArrayNew ) ): - self.logger.error( "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 ) + attrXYZ: npt.NDArray[ np.float64 ] = np.full_like( attrArray, np.nan ) - 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. + # Get all local basis vectors + localBasis: npt.NDArray[ np.float64 ] = getLocalBasisVectors( self.outputMesh ) - Args: - vector (npt.NDArray[np.float64]): input coordinates. - changeOfBasisMatrix (npt.NDArray[np.float64]): change of basis matrix + 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 ) }." ) - 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. + # Compute attribute XYZ components + cellLocalBasis: npt.NDArray[ np.float64 ] = localBasis[ :, i, : ] + attrXYZ[ i ] = convertAttributeFromLocalToXYZForOneCell( cellAttribute, cellLocalBasis ) - 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 ) + if not np.any( np.isfinite( attrXYZ ) ): + self.logger.error( "Attribute new coordinate calculation failed." ) - def __getNormalTangentsVectors( self: Self ) -> npt.NDArray[ np.float64 ]: - """Compute the change of basis matrix from Local to XYZ bases. + return attrXYZ - 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.outputMesh.GetCellData().GetNormals() ) # type: ignore[no-untyped-call] - if normals is None: - self.logger.error( "Normal attribute was not found." ) - tangents1: npt.NDArray[ np.float64 ] = vnp.vtk_to_numpy( - self.outputMesh.GetCellData().GetTangents() ) # type: ignore[no-untyped-call] - if tangents1 is None: - self.logger.error( "Tangents attribute was not found." ) - - # Compute second tangential component - tangents2: npt.NDArray[ np.float64 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 ) - if tangents2 is None: - self.logger.error( "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 (SCU) on surface. diff --git a/geos-utils/src/geos/utils/algebraFunctions.py b/geos-utils/src/geos/utils/algebraFunctions.py index eeeb88e4..0afb84cd 100644 --- a/geos-utils/src/geos/utils/algebraFunctions.py +++ b/geos-utils/src/geos/utils/algebraFunctions.py @@ -11,8 +11,7 @@ def getAttributeMatrixFromVector( attrArray: npt.NDArray[ np.float64 ], ) -> npt.NDArray[ np.float64 ]: - r"""Get the matrix of attribute values from the vector. - + """Get the matrix of attribute values from the vector. Matrix to vector conversion is the following: * if input vector size is 3: @@ -24,7 +23,16 @@ def getAttributeMatrixFromVector( attrArray: npt.NDArray[ np.float64 ], ) -> npt 0 & 0 & v2 \end{bmatrix} - * if input vector size is 6: + * if input vector size is 9: + .. math:: + + (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} @@ -33,8 +41,8 @@ def getAttributeMatrixFromVector( attrArray: npt.NDArray[ np.float64 ], ) -> npt v5 & v4 & v3 \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. @@ -42,25 +50,42 @@ def getAttributeMatrixFromVector( attrArray: npt.NDArray[ np.float64 ], ) -> npt Returns: npt.NDArray[np.float64]: matrix of attribute values + Raises: + ValueError: The input vector must be of size 3, 9 or 6 (symmetrical case). + """ - 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 def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], size: int ) -> npt.NDArray[ np.float64 ]: - r"""Get the vector of attribute values from the matrix. + """Get the vector of attribute values from the matrix. Matrix to vector conversion is the following: @@ -74,7 +99,17 @@ def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], size: i \end{bmatrix} => (M00, M11, M22) - * otherwise: + * 3x3 Generic matrix: + .. math:: + + \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} @@ -84,20 +119,57 @@ def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], size: i \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. + """ 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)." ) + + # 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 + + +def getLocalToXYZTransformMatrix( localBasis ) -> npt.NDArray[ np.float64 ]: + """Compute and return the transformation matrix to convert from local basis to XYZ basis. + + Args: + localBasis (npt.NDArray[np.float64]): Local basis coordinates vectors in fonction of XYZ basis. + + Returns: + npt.NDArray[ np.float64 ]: Local to XYZ transformation matrix. + """ + if localBasis.shape != ( 3, 3 ): + raise ValueError( "Expected a 3x3 matrix." ) + + P: npt.NDArray[ np.float64 ] = np.transpose( localBasis ) + + return P \ No newline at end of file From 9c10aeabb44768bcabceabf0b03210b40618cb0d Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Tue, 21 Oct 2025 15:06:28 +0200 Subject: [PATCH 09/23] Refactor of normal and tangential vectors computation --- .../src/geos/mesh/utils/genericHelpers.py | 171 ++++++++++++++---- .../src/geos_posp/filters/GeosBlockMerge.py | 69 +------ 2 files changed, 147 insertions(+), 93 deletions(-) diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py index 57ee6b22..cb7483c2 100644 --- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py +++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py @@ -3,14 +3,22 @@ # SPDX-FileContributor: Martin Lemay, Paloma Martinez import numpy as np import numpy.typing as npt -from typing import Iterator, List, Sequence, Any, Union -from vtkmodules.util.numpy_support import ( - numpy_to_vtk, - vtk_to_numpy +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 +from vtkmodules.vtkCommonDataModel import ( + vtkUnstructuredGrid, + vtkMultiBlockDataSet, + vtkPolyData, + vtkDataSet, + vtkDataObject, + vtkPlane, + vtkCellTypes, + vtkIncrementalOctreePointLocator, ) -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 vtkmodules.vtkFiltersCore import ( vtk3DLinearGridPlaneCutter, vtkPolyDataNormals, vtkPolyDataTangents ) +from vtkmodules.vtkFiltersTexture import vtkTextureMapToPlane + from geos.mesh.utils.multiblockHelpers import ( getBlockElementIndexesFlatten, getBlockFromFlatIndex ) from geos.utils.algebraFunctions import ( @@ -19,7 +27,6 @@ getLocalToXYZTransformMatrix, ) - __doc__ = """ Generic VTK utilities. @@ -279,9 +286,9 @@ def createVertices( cellPtsCoord: list[ npt.NDArray[ np.float64 ] ], return points, cellVertexMapAll -def convertAttributeFromLocalToXYZForOneCell( vector, localBasisVectors ) -> npt.NDArray[ np.float64 ]: - """ - Convert one cell attribute from local to XYZ basis. +def convertAttributeFromLocalToXYZForOneCell( + vector: npt.NDArray[ np.float64 ], localBasisVectors: 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 \ @@ -303,48 +310,150 @@ def convertAttributeFromLocalToXYZForOneCell( vector, localBasisVectors ) -> npt # Apply transformation arrayXYZ: npt.NDArray[ np.float64 ] = transformMatrix @ matrix3x3 @ transformMatrix.T - # Convert back to GEOS type attribute - vectorXYZ: npt.NDArray[ np.float64 ] = getAttributeVectorFromMatrix( arrayXYZ, vector.size ) + # Convert back to GEOS type attribute and return + return getAttributeVectorFromMatrix( arrayXYZ, vector.size ) - return vectorXYZ +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. -def getNormalVectors( dataset, ) -> npt.NDArray[ np.float64 ]: - normals: npt.NDArray[ np.float64 ] = vtk_to_numpy( - dataset.GetCellData().GetNormals() ) # type: ignore[no-untyped-call] - # TODO: if normals not defined, compute it on the fly + Returns: + npt.NDArray[ np.float64 ]: The normal vectors of the input surface. + """ + normals: Union[ npt.NDArray[ np.float64 ], + None ] = vtk_to_numpy( surface.GetCellData().GetNormals() ) # type: ignore[no-untyped-call] if normals is None: - raise ValueError( "Normal attribute was not found." ) + raise ValueError( "No normal attribute found in the mesh. Use the computeNormals function beforehand." ) return normals -def getTangentsVectors( dataset, normals = None ) -> npt.NDArray[ np.float64 ]: +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 - tangents1: npt.NDArray[ np.float64 ] = vtk_to_numpy( - dataset.GetCellData().GetTangents() ) # type: ignore[no-untyped-call] + tangents1: Union[ npt.NDArray[ np.float64 ], + None ] = vtk_to_numpy( surface.GetCellData().GetTangents() ) # type: ignore[no-untyped-call] if tangents1 is None: - raise ValueError( "Tangents attribute was not found." ) + raise ValueError( "No tangential attribute found in the mesh. Use the computeTangents function beforehand." ) # Compute second tangential component - if normals is None: - normals = getNormalVectors( dataset ) + normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface ) - tangents2: npt.NDArray[ np.float64 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 ) + tangents2: Union[ npt.NDArray[ np.float64 ], None ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 ) if tangents2 is None: raise ValueError( "Local basis third axis was not computed." ) return ( tangents1, tangents2 ) -def getLocalBasisVectors( dataset: vtkPolyData ) -> npt.NDArray[ np.float64 ]: - """Return the local basis vectors for all cells of the input dataset. +def getLocalBasisVectors( surface: vtkPolyData ) -> npt.NDArray[ np.float64 ]: + """Return the local basis vectors for all cells of the input surface. Args: - dataset(vtkPolydata): The input dataset + surface(vtkPolydata): The input surface. + Returns: + Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]]: Array with normal, tangential 1 and tangential 2 vectors. """ - normals = getNormalVectors( dataset ) - tangents = getTangentsVectors( dataset, normals=normals ) + try: + normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface ) + except ValueError: + surfaceWithNormals: vtkPolyData = computeNormals( surface ) + normals = getNormalVectors( surfaceWithNormals ) + + try: + tangents: Tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ] = getTangentsVectors( surface ) + except ValueError: + surfaceWithTangents: vtkPolyData = computeTangents( surface ) + tangents = getTangentsVectors( surfaceWithTangents ) return np.array( ( normals, *tangents ) ) + + +def computeNormals( surface: vtkPolyData ) -> vtkPolyData: + """Compute and set the normals of a given surface. + + Args: + surface (vtkPolyData): The input surface. + + Returns: + vtkPolyData: The surface with normal attribute. + """ + normalFilter: vtkPolyDataNormals = vtkPolyDataNormals() + normalFilter.SetInputData( surface ) + normalFilter.ComputeCellNormalsOn() + normalFilter.ComputePointNormalsOff() + normalFilter.Update() + + return normalFilter.GetOutput() + + +def computeTangents( triangulatedSurface: vtkPolyData ) -> 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. + + 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 + 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 + + +def computeSurfaceTextureCoordinates( surface: vtkPolyData ) -> 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. + + Returns: + vtkPolyData: The input surface with generated texture map. + """ + # Need to compute texture coordinates required for tangent calculation + textureFilter: vtkTextureMapToPlane = vtkTextureMapToPlane() + textureFilter.SetInputData( surface ) + textureFilter.AutomaticPlaneGenerationOn() + textureFilter.Update() + + return textureFilter.GetOutput() diff --git a/geos-posp/src/geos_posp/filters/GeosBlockMerge.py b/geos-posp/src/geos_posp/filters/GeosBlockMerge.py index 8d24b593..154d8cac 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 @@ -397,10 +395,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 @@ -437,56 +435,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 From 4643add850a84c2762a30b979334aacee0f9c2d0 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Wed, 29 Oct 2025 09:09:16 +0100 Subject: [PATCH 10/23] Replace and add tests for attribute to vector functions --- .../src/geos/utils/GeosOutputsConstants.py | 2 +- geos-utils/src/geos/utils/algebraFunctions.py | 84 ++++++++----------- .../src/geos/utils/geometryFunctions.py | 4 + geos-utils/tests/test_algebraFunctions.py | 73 ++++++++++++++++ geos-utils/tests/testsFunctionsGeosUtils.py | 26 ------ 5 files changed, 112 insertions(+), 77 deletions(-) create mode 100644 geos-utils/tests/test_algebraFunctions.py delete mode 100644 geos-utils/tests/testsFunctionsGeosUtils.py diff --git a/geos-utils/src/geos/utils/GeosOutputsConstants.py b/geos-utils/src/geos/utils/GeosOutputsConstants.py index 1ad99fbf..4b10574a 100644 --- a/geos-utils/src/geos/utils/GeosOutputsConstants.py +++ b/geos-utils/src/geos/utils/GeosOutputsConstants.py @@ -299,4 +299,4 @@ def getAttributeToTransferFromInitialTime() -> dict[ str, str ]: PostProcessingOutputsEnum.YOUNG_MODULUS_INITIAL.attributeName, PostProcessingOutputsEnum.POISSON_RATIO.attributeName: PostProcessingOutputsEnum.POISSON_RATIO_INITIAL.attributeName, - } \ No newline at end of file + } diff --git a/geos-utils/src/geos/utils/algebraFunctions.py b/geos-utils/src/geos/utils/algebraFunctions.py index 0afb84cd..876f8887 100644 --- a/geos-utils/src/geos/utils/algebraFunctions.py +++ b/geos-utils/src/geos/utils/algebraFunctions.py @@ -1,45 +1,46 @@ # 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. """ def getAttributeMatrixFromVector( attrArray: npt.NDArray[ np.float64 ], ) -> npt.NDArray[ np.float64 ]: - """Get the matrix of attribute values from the vector. + r"""Get the matrix of attribute values from the vector. + Matrix to vector conversion is the following: * 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 9: .. math:: - (v1, v2, v3, v4, v5, v6, v7, v8, v9) => \begin{bmatrix} + (v1, v2, v3, v4, v5, v6, v7, v8, v9) => \\begin{bmatrix} v1 & v6 & v5 \\ v9 & v2 & v4 \\ v8 & v7 & v3 - \end{bmatrix} + \\end{bmatrix} * if input vector size is 6 (symmetrical tensor): .. math:: - (v1, v2, v3, v4, v5, v6) => \begin{bmatrix} + (v1, v2, v3, v4, v5, v6) => \\begin{bmatrix} v1 & v6 & v5 \\ v6 & v2 & v4 \\ v5 & v4 & v3 - \end{bmatrix} + \\end{bmatrix} .. Note:: In the case of 3 x 3 tensors, GEOS is currently only implemented for symmetrical cases. @@ -47,12 +48,11 @@ def getAttributeMatrixFromVector( attrArray: npt.NDArray[ np.float64 ], ) -> npt 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. """ if attrArray.size not in ( 3, 6, 9 ): raise ValueError( "Vectorial attribute must contain 3, 6 or 9 components." ) @@ -62,64 +62,63 @@ def getAttributeMatrixFromVector( attrArray: npt.NDArray[ np.float64 ], ) -> npt # shear stress components if attrArray.size == 6: - matrix[ 0, 1 ] = attrArray[ 5 ] # v1 + matrix[ 0, 1 ] = attrArray[ 5 ] # v1 matrix[ 1, 0 ] = attrArray[ 5 ] - matrix[ 0, 2 ] = attrArray[ 4 ] # v5 + matrix[ 0, 2 ] = attrArray[ 4 ] # v5 matrix[ 2, 0 ] = attrArray[ 4 ] - matrix[ 1, 2 ] = attrArray[ 3 ] # v4 + 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, 1 ] = attrArray[ 5 ] # v1 + matrix[ 1, 0 ] = attrArray[ 8 ] # v9 - matrix[ 0, 2 ] = attrArray[ 4 ] # v5 - matrix[ 2, 0 ] = attrArray[ 7 ] # v8 + matrix[ 0, 2 ] = attrArray[ 4 ] # v5 + matrix[ 2, 0 ] = attrArray[ 7 ] # v8 - matrix[ 1, 2 ] = attrArray[ 3 ] # v4 - matrix[ 2, 1 ] = attrArray[ 6 ] # v7 + matrix[ 1, 2 ] = attrArray[ 3 ] # v4 + matrix[ 2, 1 ] = attrArray[ 6 ] # v7 return matrix def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], size: int ) -> npt.NDArray[ np.float64 ]: - """Get the vector of attribute values from the matrix. + r"""Get the vector of attribute values from the matrix. Matrix to vector conversion is the following: * 3x3 diagonal matrix: .. math:: - \begin{bmatrix} + \\begin{bmatrix} M00 & 0 & 0 \\ 0 & M11 & 0 \\ 0 & 0 & M22 - \end{bmatrix} + \\end{bmatrix} => (M00, M11, M22) * 3x3 Generic matrix: .. math:: - \begin{bmatrix} + \\begin{bmatrix} M00 & M01 & M02 \\ M10 & M11 & M12 \\ M20 & M21 & M22 - \end{matrix} + \\end{matrix} => (M00, M11, M22, M12, M02, M01, M21, M20, M10) * Symmetric case .. math:: - \begin{bmatrix} + \\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. @@ -127,16 +126,18 @@ def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], size: i attrMatrix (npt.NDArray[np.float64]): Matrix of attribute values. 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 ) 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 ) @@ -156,20 +157,3 @@ def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], size: i attrArray[ 8 ] = attrMatrix[ 1, 0 ] return attrArray - - -def getLocalToXYZTransformMatrix( localBasis ) -> npt.NDArray[ np.float64 ]: - """Compute and return the transformation matrix to convert from local basis to XYZ basis. - - Args: - localBasis (npt.NDArray[np.float64]): Local basis coordinates vectors in fonction of XYZ basis. - - Returns: - npt.NDArray[ np.float64 ]: Local to XYZ transformation matrix. - """ - if localBasis.shape != ( 3, 3 ): - raise ValueError( "Expected a 3x3 matrix." ) - - P: npt.NDArray[ np.float64 ] = np.transpose( localBasis ) - - return P \ No newline at end of file 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..7d3effba --- /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 ) ) ) From 64494c8197ccd4d129d14dd7373d638921a5330c Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Wed, 29 Oct 2025 09:17:40 +0100 Subject: [PATCH 11/23] Add tests and better error handling --- .../src/geos/mesh/utils/genericHelpers.py | 153 +++++++-- .../test_genericHelpers_normalAndTangents.py | 302 ++++++++++++++++++ 2 files changed, 423 insertions(+), 32 deletions(-) create mode 100644 geos-mesh/tests/test_genericHelpers_normalAndTangents.py diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py index 9d9a61fa..5297c960 100644 --- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py +++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py @@ -2,10 +2,11 @@ # 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, Tuple from vtkmodules.util.numpy_support import ( numpy_to_vtk, vtk_to_numpy ) -from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints, reference, vtkDataArray +from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints, reference, vtkDataArray, vtkLogger from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkMultiBlockDataSet, @@ -24,8 +25,10 @@ from geos.utils.algebraFunctions import ( getAttributeMatrixFromVector, getAttributeVectorFromMatrix, - getLocalToXYZTransformMatrix, ) +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. @@ -287,7 +290,9 @@ def createVertices( cellPtsCoord: list[ npt.NDArray[ np.float64 ] ], def convertAttributeFromLocalToXYZForOneCell( - vector: npt.NDArray[ np.float64 ], localBasisVectors: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]: + 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:: @@ -296,7 +301,7 @@ def convertAttributeFromLocalToXYZForOneCell( 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 (np.NDArray[np.float64]): The local basis vectors expressed in the XYZ basis. + 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. @@ -305,7 +310,7 @@ def convertAttributeFromLocalToXYZForOneCell( matrix3x3: npt.NDArray[ np.float64 ] = getAttributeMatrixFromVector( vector ) # Get transform matrix - transformMatrix: npt.NDArray[ np.float64 ] = getLocalToXYZTransformMatrix( localBasisVectors ) + transformMatrix: npt.NDArray[ np.float64 ] = getChangeOfBasisMatrix( localBasisVectors, CANONICAL_BASIS_3D ) # Apply transformation arrayXYZ: npt.NDArray[ np.float64 ] = transformMatrix @ matrix3x3 @ transformMatrix.T @@ -327,11 +332,11 @@ def getNormalVectors( surface: vtkPolyData ) -> npt.NDArray[ np.float64 ]: npt.NDArray[ np.float64 ]: The normal vectors of the input surface. """ normals: Union[ npt.NDArray[ np.float64 ], - None ] = vtk_to_numpy( surface.GetCellData().GetNormals() ) # type: ignore[no-untyped-call] + 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 normals + return vtk_to_numpy( normals ) def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ]: @@ -349,9 +354,11 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64 """ # Get first tangential component tangents1: Union[ npt.NDArray[ np.float64 ], - None ] = vtk_to_numpy( surface.GetCellData().GetTangents() ) # type: ignore[no-untyped-call] + None ] = surface.GetCellData().GetTangents() # type: ignore[no-untyped-call] if tangents1 is None: raise ValueError( "No tangential attribute found in the mesh. Use the computeTangents function beforehand." ) + else: + tangents1 = vtk_to_numpy( tangents1 ) # Compute second tangential component normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface ) @@ -374,44 +381,82 @@ def getLocalBasisVectors( surface: vtkPolyData ) -> npt.NDArray[ np.float64 ]: """ try: normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface ) + surfaceWithNormals: vtkPolyData = surface + # ValueError raised if no normals found in the mesh except ValueError: - surfaceWithNormals: vtkPolyData = computeNormals( surface ) + # 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( surface ) + 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 ValueError: - surfaceWithTangents: vtkPolyData = computeTangents( surface ) + surfaceWithTangents: vtkPolyData = computeTangents( surfaceWithNormals ) tangents = getTangentsVectors( surfaceWithTangents ) return np.array( ( normals, *tangents ) ) -def computeNormals( surface: vtkPolyData ) -> vtkPolyData: +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. """ - normalFilter: vtkPolyDataNormals = vtkPolyDataNormals() - normalFilter.SetInputData( surface ) - normalFilter.ComputeCellNormalsOn() - normalFilter.ComputePointNormalsOff() - normalFilter.Update() + 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() ) - return normalFilter.GetOutput() + outputSurface = normalFilter.GetOutput() + if outputSurface.GetCellData().GetNormals() is None: + raise VTKError( "Something went wrong with VTK calculation of the normals." ) -def computeTangents( triangulatedSurface: vtkPolyData ) -> vtkPolyData: + 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 @@ -422,18 +467,39 @@ def computeTangents( triangulatedSurface: vtkPolyData ) -> vtkPolyData: # TODO: triangulate the surface before computation of the tangents if needed. # 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." + 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() - assert array is not None, "Attribute Tangents is not in the mesh." + if array is None: + raise VTKError( "Attribute Tangents is not in the mesh." ) surface1.GetCellData().SetTangents( array ) surface1.GetCellData().Modified() @@ -441,19 +507,42 @@ def computeTangents( triangulatedSurface: vtkPolyData ) -> vtkPolyData: return surface1 -def computeSurfaceTextureCoordinates( surface: vtkPolyData ) -> vtkPolyData: +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 - textureFilter: vtkTextureMapToPlane = vtkTextureMapToPlane() - textureFilter.SetInputData( surface ) - textureFilter.AutomaticPlaneGenerationOn() - textureFilter.Update() + 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/test_genericHelpers_normalAndTangents.py b/geos-mesh/tests/test_genericHelpers_normalAndTangents.py new file mode 100644 index 00000000..22acccc4 --- /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( ValueError ): + 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 ) From 98a825b6f1f9ae1f28e7cfebf558722463a66f9e Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Wed, 29 Oct 2025 09:18:40 +0100 Subject: [PATCH 12/23] Adding tests for SurfaceGeomechanics filter --- geos-mesh/tests/conftest.py | 2 +- .../post_processing/SurfaceGeomechanics.py | 110 ++++++++---------- .../tests/test_SurfaceGeomechanics.py | 84 +++++++++++++ .../geos/pv/plugins/PVSurfaceGeomechanics.py | 7 +- 4 files changed, 133 insertions(+), 70 deletions(-) create mode 100644 geos-processing/tests/test_SurfaceGeomechanics.py diff --git a/geos-mesh/tests/conftest.py b/geos-mesh/tests/conftest.py index 156d37f4..05fd5221 100644 --- a/geos-mesh/tests/conftest.py +++ b/geos-mesh/tests/conftest.py @@ -10,7 +10,7 @@ import numpy.typing as npt from vtkmodules.vtkCommonDataModel import vtkDataSet, vtkMultiBlockDataSet, vtkPolyData -from vtkmodules.vtkIOXML import vtkXMLGenericDataObjectReader, vtkXMLMultiBlockDataReader +from vtkmodules.vtkIOXML import vtkXMLGenericDataObjectReader @pytest.fixture diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py index c5dba4d9..23b68d02 100644 --- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py +++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py @@ -8,14 +8,9 @@ from typing_extensions import Self, Union import numpy.typing as npt -import vtkmodules.util.numpy_support as vnp from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase -from vtkmodules.vtkCommonCore import ( - vtkDataArray, - vtkDoubleArray, -) -from vtkmodules.vtkCommonDataModel import ( - vtkPolyData, ) +from vtkmodules.vtkCommonCore import vtkDataArray +from vtkmodules.vtkCommonDataModel import vtkPolyData from geos.mesh.utils.arrayModifiers import createAttribute from geos.mesh.utils.arrayHelpers import ( @@ -38,7 +33,6 @@ GeosMeshOutputsEnum, PostProcessingOutputsEnum, ) -import geos.utils.geometryFunctions as geom __doc__ = """ SurfaceGeomechanics is a VTK filter that allows: @@ -71,24 +65,24 @@ speHandler: bool # defaults to false # Instantiate the filter - filter: SurfaceGeomechanics = SurfaceGeomechanics( inputMesh, speHandler ) + sg: SurfaceGeomechanics = SurfaceGeomechanics( inputMesh, speHandler ) # Set the handler of yours (only if speHandler is True). yourHandler: logging.Handler - filter.SetLoggerHandler( yourHandler ) + sg.SetLoggerHandler( yourHandler ) # [optional] Set rock cohesion and friction angle - filter.SetRockCohesion( rockCohesion ) - filter.SetFrictionAngle( frictionAngle ) + sg.SetRockCohesion( rockCohesion ) + sg.SetFrictionAngle( frictionAngle ) # Do calculations - filter.applyFilter() + sg.applyFilter() # Get output object - output: vtkPolyData = filter.GetOutputMesh() + output: vtkPolyData = sg.GetOutputMesh() # Get created attribute names - newAttributeNames: set[str] = filter.GetNewAttributeNames() + newAttributeNames: set[str] = sg.GetNewAttributeNames() .. Note:: @@ -100,11 +94,10 @@ """ loggerTitle: str = "Surface Geomechanics" + class SurfaceGeomechanics( VTKPythonAlgorithmBase ): - def __init__( self: Self, - surfacicMesh: vtkPolyData, - speHandler: bool = False ) -> None: + 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 @@ -126,18 +119,18 @@ def __init__( self: Self, self.logger.setLevel( logging.INFO ) # Input surfacic mesh - if surfacicMesh is None: - self.logger.error( "Input surface is undefined." ) 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 = None + 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, } + self.attributesToConvert: set[ str ] = { + GeosMeshOutputsEnum.DISPLACEMENT_JUMP.attributeName, + } # Attributes are either on points or on cells self.attributeOnPoints: bool = False @@ -148,7 +141,6 @@ def __init__( self: Self, # New created attributes names self.newAttributeNames: set[ str ] = set() - def SetLoggerHandler( self: Self, handler: Logger ) -> None: """Set a specific handler for the filter logger. @@ -164,7 +156,7 @@ def SetLoggerHandler( self: Self, handler: Logger ) -> None: "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 ): + def SetSurfaceName( self: Self, name: str ) -> None: """Set a name for the input surface. For logging purpose only. Args: @@ -188,17 +180,14 @@ def SetFrictionAngle( self: Self, frictionAngle: float ) -> None: """ 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. @@ -207,7 +196,6 @@ def GetConvertAttributes( self: Self ) -> bool: """ return self.convertAttributesOn - def GetAttributesToConvert( self: Self ) -> set[ str ]: """Get the set of attributes that will be converted from local to XYZ basis. @@ -216,7 +204,6 @@ def GetAttributesToConvert( self: Self ) -> set[ str ]: """ 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. @@ -230,17 +217,15 @@ def SetAttributesToConvert( self: Self, attributesToConvert: set[ str ] ) -> Non 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. + attributeName (Union[list[str],set[str]]): List of the attribute array names. """ - self.attributesToConvert.add( attributeName ) + self.attributesToConvert.update( attributeName ) self.ConvertAttributesOn() - def GetNewAttributeNames( self: Self ) -> set[ str ]: """Get the set of new attribute names that were created. @@ -249,7 +234,6 @@ def GetNewAttributeNames( self: Self ) -> set[ str ]: """ return self.newAttributeNames - def applyFilter( self: Self ) -> bool: """Compute Geomechanical properties on input surface. @@ -269,7 +253,7 @@ def applyFilter( self: Self ) -> bool: # Conversion of attributes from Normal/Tangent basis to xyz basis if self.convertAttributesOn: - self.logger.info( "Conversion of attributes from local to XYZ basis.") + 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 @@ -282,7 +266,6 @@ def applyFilter( self: Self ) -> bool: 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. @@ -293,7 +276,7 @@ def convertAttributesFromLocalToXYZBasis( self: Self ) -> bool: attributesToConvert: set[ str ] = self.__filterAttributesToConvert() if len( attributesToConvert ) == 0: - self.logger.warning( f"No attribute to convert from local to XYZ basis were found." ) + 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 @@ -305,25 +288,25 @@ def convertAttributesFromLocalToXYZBasis( self: Self ) -> bool: continue if self.attributeOnPoints: - self.logger.error( "This filter can only convert cell attributes from local to XYZ basis, not point attributes." ) - localArray: vtkDoubleArray = getArrayInObject( self.outputMesh, attrNameLocal, 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): + 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. @@ -340,9 +323,11 @@ def __filterAttributesToConvert( self: Self ) -> set[ str ]: 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." ) + 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.") + self.logger.warning( f"Attribute {attrName} not in the input mesh." ) if len( attributesFiltered ) == 0: self.logger.warning( "All attributes filtered out." ) @@ -370,7 +355,9 @@ def __computeXYZCoordinates( 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 ) }." ) + 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, : ] @@ -381,7 +368,6 @@ def __computeXYZCoordinates( return attrXYZ - def computeShearCapacityUtilization( self: Self ) -> bool: """Compute the shear capacity utilization (SCU) on surface. @@ -395,27 +381,23 @@ def computeShearCapacityUtilization( self: Self ) -> bool: tractionAttributeName: str = GeosMeshOutputsEnum.TRACTION.attributeName traction: npt.NDArray[ np.float64 ] = getArrayInObject( self.outputMesh, tractionAttributeName, self.attributeOnPoints ) - if traction is None: - self.logger.error( f"{tractionAttributeName} attribute is undefined." ) - self.logger.error( f"Failed to compute {SCUAttributeName}." ) # Computation of the shear capacity utilization (SCU) - scuAttribute: npt.NDArray[ np.float64 ] = fcts.shearCapacityUtilization( - traction, self.rockCohesion, self.frictionAngle ) + # 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.outputMesh, scuAttribute, SCUAttributeName, (), self.attributeOnPoints, logger=self.logger ): self.logger.error( f"Failed to create attribute {SCUAttributeName}." ) return False else: - self.logger.info( f"SCU computed and added to the output mesh." ) + self.logger.info( "SCU computed and added to the output mesh." ) self.newAttributeNames.add( SCUAttributeName ) return True diff --git a/geos-processing/tests/test_SurfaceGeomechanics.py b/geos-processing/tests/test_SurfaceGeomechanics.py new file mode 100644 index 00000000..a940f903 --- /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 ) + + filter: SurfaceGeomechanics = SurfaceGeomechanics( testCase.mesh ) + + assert filter.applyFilter() + + mesh: vtkPolyData = filter.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 ) + filter: SurfaceGeomechanics = SurfaceGeomechanics( failingCase.mesh ) + with pytest.raises( AssertionError ): + assert filter.applyFilter() diff --git a/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py index ff198ba2..a5a2a26c 100644 --- a/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py +++ b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py @@ -11,8 +11,7 @@ VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, ) from paraview.detail.loghandler import ( # type: ignore[import-not-found] - VTKHandler, -) + VTKHandler, ) # update sys.path to load all GEOS Python Package dependencies geos_pv_path: Path = Path( __file__ ).parent.parent.parent.parent.parent @@ -30,9 +29,6 @@ getBlockElementIndexesFlatten, getBlockFromFlatIndex, ) -from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] - VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, -) from vtkmodules.vtkCommonCore import ( vtkDataArray, vtkInformation, @@ -64,6 +60,7 @@ """ + @smproxy.filter( name="PVSurfaceGeomechanics", label="Geos Surface Geomechanics" ) @smhint.xml( '' ) @smproperty.input( name="Input", port_index=0 ) From 52bcd566c4de16537a849b9141293d6f106f4f5d Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Wed, 29 Oct 2025 10:12:38 +0100 Subject: [PATCH 13/23] Fix --- docs/geos_posp_docs/PVplugins.rst | 6 ------ docs/geos_posp_docs/filters.rst | 8 -------- docs/geos_processing_docs/generic_processing_tool.rst | 4 ---- .../src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py | 2 +- .../PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py | 2 +- install_packages.sh | 2 +- 6 files changed, 3 insertions(+), 21 deletions(-) delete mode 100644 docs/geos_processing_docs/generic_processing_tool.rst 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 b3c3cd89..f0c00c1c 100644 --- a/docs/geos_posp_docs/filters.rst +++ b/docs/geos_posp_docs/filters.rst @@ -18,11 +18,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/generic_processing_tool.rst b/docs/geos_processing_docs/generic_processing_tool.rst deleted file mode 100644 index 4fdb32fc..00000000 --- a/docs/geos_processing_docs/generic_processing_tool.rst +++ /dev/null @@ -1,4 +0,0 @@ -Generic processing filters -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -In progress.. \ No newline at end of file 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/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 From d3e4fc6fbab6df3ff66c5aa979429aa0438b1986 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Wed, 29 Oct 2025 11:33:26 +0100 Subject: [PATCH 14/23] typo --- geos-utils/tests/test_algebraFunctions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geos-utils/tests/test_algebraFunctions.py b/geos-utils/tests/test_algebraFunctions.py index 7d3effba..104a8136 100644 --- a/geos-utils/tests/test_algebraFunctions.py +++ b/geos-utils/tests/test_algebraFunctions.py @@ -5,7 +5,7 @@ import numpy as np import numpy.typing as npt from unittest import TestCase -from typing.extensions import Self +from typing_extensions import Self from geos.utils.algebraFunctions import ( getAttributeMatrixFromVector, From badc7cd7227d206f8f779c9bc67cafa22ef519c7 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Thu, 30 Oct 2025 10:15:55 +0100 Subject: [PATCH 15/23] Small modifications from Romain review --- docs/geos_processing_docs/post_processing.rst | 2 ++ .../post_processing/SurfaceGeomechanics.py | 11 +++++++---- .../tests/test_SurfaceGeomechanics.py | 10 +++++----- .../geos/pv/plugins/PVSurfaceGeomechanics.py | 18 +++++++++--------- 4 files changed, 23 insertions(+), 18 deletions(-) diff --git a/docs/geos_processing_docs/post_processing.rst b/docs/geos_processing_docs/post_processing.rst index 8698539a..d22b46a0 100644 --- a/docs/geos_processing_docs/post_processing.rst +++ b/docs/geos_processing_docs/post_processing.rst @@ -13,6 +13,8 @@ GEOS computes many outputs including flow and geomechanic properties if coupled 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° diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py index 23b68d02..2f527efd 100644 --- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py +++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py @@ -1,6 +1,6 @@ # 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 logging import numpy as np @@ -86,16 +86,19 @@ .. 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`. + + - 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( VTKPythonAlgorithmBase ): +class SurfaceGeomechanics: def __init__( self: Self, surfacicMesh: vtkPolyData, speHandler: bool = False ) -> None: """Vtk filter to compute geomechanical surfacic attributes. diff --git a/geos-processing/tests/test_SurfaceGeomechanics.py b/geos-processing/tests/test_SurfaceGeomechanics.py index a940f903..7623e802 100644 --- a/geos-processing/tests/test_SurfaceGeomechanics.py +++ b/geos-processing/tests/test_SurfaceGeomechanics.py @@ -67,11 +67,11 @@ def test_SurfaceGeomechanics() -> None: """Test SurfaceGeomechanics vtk filter.""" testCase: TriangulatedSurfaceTestCase = TriangulatedSurfaceTestCase( pointsCoords, triangles, attributes ) - filter: SurfaceGeomechanics = SurfaceGeomechanics( testCase.mesh ) + sgFilter: SurfaceGeomechanics = SurfaceGeomechanics( testCase.mesh ) - assert filter.applyFilter() + assert sgFilter.applyFilter() - mesh: vtkPolyData = filter.GetOutputMesh() + mesh: vtkPolyData = sgFilter.GetOutputMesh() assert mesh.GetCellData().HasArray( "SCU" ) assert mesh.GetCellData().HasArray( "displacementJump_XYZ" ) @@ -79,6 +79,6 @@ def test_SurfaceGeomechanics() -> None: def test_failingSurfaceGeomechanics() -> None: """Test failing of SurfaceGeomechanics due to absence of attributes in the mesh.""" failingCase: TriangulatedSurfaceTestCase = TriangulatedSurfaceTestCase( pointsCoords, triangles, None ) - filter: SurfaceGeomechanics = SurfaceGeomechanics( failingCase.mesh ) + sgFilter: SurfaceGeomechanics = SurfaceGeomechanics( failingCase.mesh ) with pytest.raises( AssertionError ): - assert filter.applyFilter() + assert sgFilter.applyFilter() diff --git a/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py index a5a2a26c..50e7dc34 100644 --- a/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py +++ b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py @@ -154,19 +154,19 @@ def RequestData( for blockIndex in surfaceBlockIndexes: surfaceBlock: vtkPolyData = vtkPolyData.SafeDownCast( getBlockFromFlatIndex( output, blockIndex ) ) - filter: SurfaceGeomechanics = SurfaceGeomechanics( surfaceBlock, True ) - filter.SetSurfaceName( f"blockIndex {blockIndex}" ) - if not filter.logger.hasHandlers(): - filter.SetLoggerHandler( VTKHandler() ) + sgFilter: SurfaceGeomechanics = SurfaceGeomechanics( surfaceBlock, True ) + sgFilter.SetSurfaceName( f"blockIndex {blockIndex}" ) + if not sgFilter.logger.hasHandlers(): + sgFilter.SetLoggerHandler( VTKHandler() ) - filter.SetRockCohesion( self._getRockCohesion() ) - filter.SetFrictionAngle( self._getFrictionAngle() ) - filter.applyFilter() + sgFilter.SetRockCohesion( self._getRockCohesion() ) + sgFilter.SetFrictionAngle( self._getFrictionAngle() ) + sgFilter.applyFilter() - outputSurface: vtkPolyData = filter.GetOutputMesh() + outputSurface: vtkPolyData = sgFilter.GetOutputMesh() # add attributes to output surface mesh - for attributeName in filter.GetNewAttributeNames(): + for attributeName in sgFilter.GetNewAttributeNames(): attr: vtkDataArray = outputSurface.GetCellData().GetArray( attributeName ) surfaceBlock.GetCellData().AddArray( attr ) surfaceBlock.GetCellData().Modified() From 335782c2898585fc43f7526d0e9a44475c332e4b Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Fri, 31 Oct 2025 11:15:10 +0100 Subject: [PATCH 16/23] Fix typing and error in getTangents function --- .../src/geos/mesh/utils/genericHelpers.py | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py index 5297c960..e4d66506 100644 --- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py +++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py @@ -6,7 +6,7 @@ import numpy.typing as npt 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 +from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints, reference, vtkDataArray, vtkLogger, vtkFloatArray from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkMultiBlockDataSet, @@ -353,19 +353,19 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64 Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: The tangents vectors of the input surface. """ # Get first tangential component - tangents1: Union[ npt.NDArray[ np.float64 ], - None ] = surface.GetCellData().GetTangents() # type: ignore[no-untyped-call] - if tangents1 is None: - raise ValueError( "No tangential attribute found in the mesh. Use the computeTangents function beforehand." ) - else: - tangents1 = vtk_to_numpy( tangents1 ) + vtkTangents: Union[ vtkFloatArray, None ] = surface.GetCellData().GetTangents() # type: ignore[no-untyped-call] + tangents1: npt.NDArray[ np.float64 ] - # Compute second tangential component - normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface ) + 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: Union[ npt.NDArray[ np.float64 ], None ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 ) - if tangents2 is None: - raise ValueError( "Local basis third axis was not computed." ) + tangents2: npt.NDArray[ np.float64 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 ) return ( tangents1, tangents2 ) From 359ceb024e17f450454f2dc6ed5f8e0bd01ec389 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Fri, 31 Oct 2025 13:16:40 +0100 Subject: [PATCH 17/23] Fix tests following previous commit modifs --- geos-mesh/src/geos/mesh/utils/genericHelpers.py | 2 +- geos-mesh/tests/test_genericHelpers_normalAndTangents.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py index e4d66506..fd592984 100644 --- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py +++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py @@ -393,7 +393,7 @@ def getLocalBasisVectors( surface: vtkPolyData ) -> npt.NDArray[ np.float64 ]: 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 ValueError: + except VTKError: surfaceWithTangents: vtkPolyData = computeTangents( surfaceWithNormals ) tangents = getTangentsVectors( surfaceWithTangents ) diff --git a/geos-mesh/tests/test_genericHelpers_normalAndTangents.py b/geos-mesh/tests/test_genericHelpers_normalAndTangents.py index 22acccc4..5cfc3f27 100644 --- a/geos-mesh/tests/test_genericHelpers_normalAndTangents.py +++ b/geos-mesh/tests/test_genericHelpers_normalAndTangents.py @@ -254,7 +254,7 @@ def test_failingComputeNormals( emptySurface: vtkPolyData ) -> None: def test_failingGetTangents( emptySurface: vtkPolyData ) -> None: """Test error raising when getting the surface tangents.""" - with pytest.raises( ValueError ): + with pytest.raises( VTKError ): getTangentsVectors( emptySurface ) From 27ae09d9a624e45b95b7a18354dc622672a5b58e Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Fri, 31 Oct 2025 14:03:49 +0100 Subject: [PATCH 18/23] SISO filter --- .../post_processing/SurfaceGeomechanics.py | 1 - .../geos/pv/plugins/PVSurfaceGeomechanics.py | 57 ++++++------------- 2 files changed, 16 insertions(+), 42 deletions(-) diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py index 2f527efd..d49db1b9 100644 --- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py +++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py @@ -111,7 +111,6 @@ def __init__( self: Self, surfacicMesh: vtkPolyData, speHandler: bool = False ) speHandler (bool, optional): True to use a specific handler, False to use the internal handler. Defaults to False. """ - super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkPolyData" ) # type: ignore[no-untyped-call] # Logger self.logger: Logger diff --git a/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py index 50e7dc34..2858766e 100644 --- a/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py +++ b/geos-pv/src/geos/pv/plugins/PVSurfaceGeomechanics.py @@ -8,7 +8,7 @@ from typing_extensions import Self from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] - VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, + VTKPythonAlgorithmBase, smdomain, smproperty, ) from paraview.detail.loghandler import ( # type: ignore[import-not-found] VTKHandler, ) @@ -17,6 +17,7 @@ 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() @@ -30,10 +31,7 @@ getBlockFromFlatIndex, ) from vtkmodules.vtkCommonCore import ( - vtkDataArray, - vtkInformation, - vtkInformationVector, -) + vtkDataArray, ) from vtkmodules.vtkCommonDataModel import ( vtkMultiBlockDataSet, vtkPolyData, @@ -61,24 +59,16 @@ """ -@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 ) +@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 that contains surfaces with - Normals and Tangential attributes. + Input is a vtkMultiBlockDataSet containing surfaces. """ - super().__init__( - nInputPorts=1, - nOutputPorts=1, - inputType="vtkMultiBlockDataSet", - outputType="vtkMultiBlockDataSet", - ) # rock cohesion (Pa) self.rockCohesion: float = DEFAULT_ROCK_COHESION # friction angle (°) @@ -118,7 +108,7 @@ def a01SetRockCohesion( self: Self, value: float ) -> None: """ ) def a02SetFrictionAngle( self: Self, value: float ) -> None: - """Set frition angle. + """Set friction angle. Args: value (float): friction angle (°) @@ -126,33 +116,18 @@ def a02SetFrictionAngle( self: Self, value: float ) -> None: self.frictionAngle = value self.Modified() - def RequestData( - self: Self, - request: vtkInformation, # noqa: F841 - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestData. + def Filter( self: Self, inputMesh: vtkMultiBlockDataSet, outputMesh: vtkMultiBlockDataSet ) -> None: + """Apply SurfaceGeomechanics filter to the mesh. Args: - request (vtkInformation): Request - inInfoVec (list[vtkInformationVector]): Input objects - outInfoVec (vtkInformationVector): Output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. + inputMesh (vtkMultiBlockDataSet): The input multiblock mesh with surfaces. + outputMesh (vtkMultiBlockDataSet): The output multiblock mesh with converted attributes and SCU. """ - inputMesh: vtkMultiBlockDataSet = vtkMultiBlockDataSet.GetData( inInfoVec[ 0 ] ) - output: vtkMultiBlockDataSet = self.GetOutputData( outInfoVec, 0 ) - - assert inputMesh is not None, "Input surface is null." - assert output is not None, "Output pipeline is null." - - output.ShallowCopy( inputMesh ) + outputMesh.ShallowCopy( inputMesh ) surfaceBlockIndexes: list[ int ] = getBlockElementIndexesFlatten( inputMesh ) for blockIndex in surfaceBlockIndexes: - surfaceBlock: vtkPolyData = vtkPolyData.SafeDownCast( getBlockFromFlatIndex( output, blockIndex ) ) + surfaceBlock: vtkPolyData = vtkPolyData.SafeDownCast( getBlockFromFlatIndex( outputMesh, blockIndex ) ) sgFilter: SurfaceGeomechanics = SurfaceGeomechanics( surfaceBlock, True ) sgFilter.SetSurfaceName( f"blockIndex {blockIndex}" ) @@ -172,8 +147,8 @@ def RequestData( surfaceBlock.GetCellData().Modified() surfaceBlock.Modified() - output.Modified() - return 1 + outputMesh.Modified() + return def _getFrictionAngle( self: Self ) -> float: """Get friction angle in radian. From 58e30e332ee8e3ba9ffd447200b678611dd81222 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Fri, 31 Oct 2025 14:10:00 +0100 Subject: [PATCH 19/23] fix docstring --- geos-mesh/src/geos/mesh/utils/genericHelpers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py index fd592984..0e2ebb39 100644 --- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py +++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py @@ -377,7 +377,7 @@ def getLocalBasisVectors( surface: vtkPolyData ) -> npt.NDArray[ np.float64 ]: surface(vtkPolydata): The input surface. Returns: - Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]]: Array with normal, tangential 1 and tangential 2 vectors. + npt.NDArray[np.float64]: Array with normal, tangential 1 and tangential 2 vectors. """ try: normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface ) From c56b2de25c2c8a4b8cf4d52eba36a879094ce028 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Fri, 31 Oct 2025 14:44:47 +0100 Subject: [PATCH 20/23] fix the fix --- .../src/geos/processing/post_processing/SurfaceGeomechanics.py | 2 -- geos-processing/tests/test_SurfaceGeomechanics.py | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py index d49db1b9..90320e82 100644 --- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py +++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py @@ -8,7 +8,6 @@ from typing_extensions import Self, Union import numpy.typing as npt -from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase from vtkmodules.vtkCommonCore import vtkDataArray from vtkmodules.vtkCommonDataModel import vtkPolyData @@ -111,7 +110,6 @@ def __init__( self: Self, surfacicMesh: vtkPolyData, speHandler: bool = False ) 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: diff --git a/geos-processing/tests/test_SurfaceGeomechanics.py b/geos-processing/tests/test_SurfaceGeomechanics.py index 7623e802..dfdb030a 100644 --- a/geos-processing/tests/test_SurfaceGeomechanics.py +++ b/geos-processing/tests/test_SurfaceGeomechanics.py @@ -14,7 +14,7 @@ from vtkmodules.vtkCommonCore import vtkPoints, vtkFloatArray import vtkmodules.util.numpy_support as vnp from geos.processing.post_processing.SurfaceGeomechanics import SurfaceGeomechanics - +from geos.utils.Errors import VTKError @dataclass class TriangulatedSurfaceTestCase: From cdd923cd3e918608fd5e0bdde56f4cb8553904f9 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Fri, 31 Oct 2025 14:51:14 +0100 Subject: [PATCH 21/23] typing --- geos-processing/tests/test_SurfaceGeomechanics.py | 1 - 1 file changed, 1 deletion(-) diff --git a/geos-processing/tests/test_SurfaceGeomechanics.py b/geos-processing/tests/test_SurfaceGeomechanics.py index dfdb030a..8ebda0e7 100644 --- a/geos-processing/tests/test_SurfaceGeomechanics.py +++ b/geos-processing/tests/test_SurfaceGeomechanics.py @@ -14,7 +14,6 @@ from vtkmodules.vtkCommonCore import vtkPoints, vtkFloatArray import vtkmodules.util.numpy_support as vnp from geos.processing.post_processing.SurfaceGeomechanics import SurfaceGeomechanics -from geos.utils.Errors import VTKError @dataclass class TriangulatedSurfaceTestCase: From 24e4c5b9655b0f1b91ba5b09f60711c80013b3ae Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Fri, 31 Oct 2025 15:03:24 +0100 Subject: [PATCH 22/23] typing --- geos-processing/tests/test_SurfaceGeomechanics.py | 1 + 1 file changed, 1 insertion(+) diff --git a/geos-processing/tests/test_SurfaceGeomechanics.py b/geos-processing/tests/test_SurfaceGeomechanics.py index 8ebda0e7..7623e802 100644 --- a/geos-processing/tests/test_SurfaceGeomechanics.py +++ b/geos-processing/tests/test_SurfaceGeomechanics.py @@ -15,6 +15,7 @@ import vtkmodules.util.numpy_support as vnp from geos.processing.post_processing.SurfaceGeomechanics import SurfaceGeomechanics + @dataclass class TriangulatedSurfaceTestCase: pointsCoords: npt.NDArray[ np.float64 ] From 7e7eca119559adbbcf528000bbadf02fea83d7a9 Mon Sep 17 00:00:00 2001 From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com> Date: Fri, 31 Oct 2025 15:53:19 +0100 Subject: [PATCH 23/23] . --- geos-mesh/tests/conftest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geos-mesh/tests/conftest.py b/geos-mesh/tests/conftest.py index 05fd5221..40f5dd46 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