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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] 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/34] . --- 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 From dab09d01969ebd642ae8e2251052492b45390250 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Fri, 31 Oct 2025 16:47:50 +0100 Subject: [PATCH 24/34] Move GeosBlockMerge from geos-posp to geos-processing --- docs/geos_posp_docs/filters.rst | 12 ------------ docs/geos_posp_docs/modules.rst | 2 -- docs/geos_processing_docs/post_processing.rst | 9 +++++++++ geos-posp/src/geos_posp/filters/__init__.py | 0 .../processing/post_processing}/GeosBlockMerge.py | 0 5 files changed, 9 insertions(+), 14 deletions(-) delete mode 100644 docs/geos_posp_docs/filters.rst delete mode 100644 geos-posp/src/geos_posp/filters/__init__.py rename {geos-posp/src/geos_posp/filters => geos-processing/src/geos/processing/post_processing}/GeosBlockMerge.py (100%) diff --git a/docs/geos_posp_docs/filters.rst b/docs/geos_posp_docs/filters.rst deleted file mode 100644 index 7bd6a88d..00000000 --- a/docs/geos_posp_docs/filters.rst +++ /dev/null @@ -1,12 +0,0 @@ -vtk Filters -=========== - -This package defines vtk filters that allows to process Geos outputs. - -geos_posp.filters.GeosBlockMerge module -------------------------------------------- - -.. automodule:: geos_posp.filters.GeosBlockMerge - :members: - :undoc-members: - :show-inheritance: diff --git a/docs/geos_posp_docs/modules.rst b/docs/geos_posp_docs/modules.rst index cf77644b..a87e24eb 100644 --- a/docs/geos_posp_docs/modules.rst +++ b/docs/geos_posp_docs/modules.rst @@ -4,6 +4,4 @@ Processing .. toctree:: :maxdepth: 5 - filters - pyvistaTools diff --git a/docs/geos_processing_docs/post_processing.rst b/docs/geos_processing_docs/post_processing.rst index 58224d9f..19b2906a 100644 --- a/docs/geos_processing_docs/post_processing.rst +++ b/docs/geos_processing_docs/post_processing.rst @@ -44,3 +44,12 @@ geos.processing.post_processing.GeosBlockExtractor module :members: :undoc-members: :show-inheritance: + + +geos.processing.post_processing.GeosBlockMerge module +----------------------------------------------------- + +.. automodule:: geos.processing.post_processing.GeosBlockMerge + :members: + :undoc-members: + :show-inheritance: diff --git a/geos-posp/src/geos_posp/filters/__init__.py b/geos-posp/src/geos_posp/filters/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/geos-posp/src/geos_posp/filters/GeosBlockMerge.py b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py similarity index 100% rename from geos-posp/src/geos_posp/filters/GeosBlockMerge.py rename to geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py From 7677e88c7bc6a12d20c40c2ea2b3f19c0612f7fd Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Fri, 31 Oct 2025 17:06:47 +0100 Subject: [PATCH 25/34] Remove PVPythonAlgorythmBase --- .../post_processing/GeosBlockMerge.py | 56 ++----------------- 1 file changed, 5 insertions(+), 51 deletions(-) diff --git a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py index 4b76d9a2..e4816b43 100644 --- a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py +++ b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py @@ -69,65 +69,24 @@ """ -class GeosBlockMerge( VTKPythonAlgorithmBase ): +class GeosBlockMerge(): - def __init__( self: Self ) -> None: + def __init__( self: Self, inputMesh: vtkMultiBlockDataSet ) -> None: """VTK Filter that perform GEOS rank merge. The filter returns a multiblock mesh composed of elementary blocks. """ - super().__init__( nInputPorts=1, nOutputPorts=1, - outputType="vtkMultiBlockDataSet" ) # type: ignore[no-untyped-call] - self.m_input: vtkMultiBlockDataSet - self.m_output: vtkMultiBlockDataSet + self.m_input: vtkMultiBlockDataSet = inputMesh + self.m_output: vtkMultiBlockDataSet = vtkMultiBlockDataSet() self.m_convertFaultToSurface: bool = True # set logger self.m_logger: Logger = getLogger( "Geos Block Merge Filter" ) - 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, - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: + def applyFilter( self: Self ) -> int: """Inherited from VTKPythonAlgorithmBase::RequestData. Args: @@ -139,11 +98,6 @@ def RequestData( int: 1 if calculation successfully ended, 0 otherwise. """ try: - self.m_input = vtkMultiBlockDataSet.GetData( inInfoVec[ 0 ] ) - - # initialize output objects -- TODO: separate it as soon as we are sure not to alter behavior - self.m_output = self.GetOutputData( outInfoVec, 0 ) # type: ignore[no-untyped-call] - self.doMerge() except ( ValueError, TypeError, RuntimeError ) as e: self.m_logger.critical( "Geos block merge failed due to:" ) From 95072db3e606c2c4187054fbf8fbae9b3be5b587 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Mon, 3 Nov 2025 16:49:54 +0100 Subject: [PATCH 26/34] clean & refactor the code --- .../post_processing/GeosBlockMerge.py | 315 ++++-------------- 1 file changed, 72 insertions(+), 243 deletions(-) diff --git a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py index e4816b43..f20d815d 100644 --- a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py +++ b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py @@ -1,36 +1,27 @@ # SPDX-License-Identifier: Apache-2.0 # # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -# SPDX-FileContributor: Martin Lemay +# SPDX-FileContributor: Martin Lemay, Romain Baville # ruff: noqa: E402 # disable Module level import not at top of file from geos.utils.GeosOutputsConstants import ( PHASE_SEP, - FluidPrefixEnum, PhaseTypeEnum, PostProcessingOutputsEnum, getRockSuffixRenaming, ) from geos.utils.Logger import Logger, getLogger from typing_extensions import Self -from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase -from vtkmodules.vtkCommonCore import ( - vtkInformation, - vtkInformationVector, -) from vtkmodules.vtkCommonDataModel import ( vtkCompositeDataSet, - vtkDataObjectTreeIterator, vtkMultiBlockDataSet, vtkPolyData, vtkUnstructuredGrid, ) -from vtkmodules.vtkFiltersCore import vtkArrayRename from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter 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.arrayModifiers import createConstantAttribute, renameAttribute from geos.mesh.utils.multiblockModifiers import mergeBlocks from geos.mesh.utils.genericHelpers import ( computeNormals, @@ -53,7 +44,7 @@ logger :Logger input :vtkMultiBlockDataSet - # instanciate the filter + # instantiate the filter mergeBlockFilter :GeosBlockMerge = GeosBlockMerge() # set the logger mergeBlockFilter.SetLogger(logger) @@ -77,35 +68,18 @@ def __init__( self: Self, inputMesh: vtkMultiBlockDataSet ) -> None: The filter returns a multiblock mesh composed of elementary blocks. """ + self.m_inputMesh: vtkMultiBlockDataSet = inputMesh + self.m_outputMesh: vtkMultiBlockDataSet = vtkMultiBlockDataSet() - self.m_input: vtkMultiBlockDataSet = inputMesh - self.m_output: vtkMultiBlockDataSet = vtkMultiBlockDataSet() - - self.m_convertFaultToSurface: bool = True + self.m_convertFaultToSurface: bool = False + self.phaseNameDict: dict[ str, set[ str ] ] = { + PhaseTypeEnum.ROCK.type: {}, + PhaseTypeEnum.FLUID.type: {}, + } # set logger self.m_logger: Logger = getLogger( "Geos Block Merge Filter" ) - def applyFilter( self: Self ) -> 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: - self.doMerge() - except ( ValueError, TypeError, RuntimeError ) as e: - self.m_logger.critical( "Geos block merge failed due to:" ) - self.m_logger.critical( e, exc_info=True ) - return 0 - else: - return 1 - def SetLogger( self: Self, logger: Logger ) -> None: """Set the logger. @@ -115,235 +89,90 @@ def SetLogger( self: Self, logger: Logger ) -> None: self.m_logger = logger def ConvertSurfaceMeshOn( self: Self ) -> None: - """Activate surface conversion from vtkUnstructredGrid to vtkPolyData.""" + """Activate surface conversion from vtkUnstructuredGrid to vtkPolyData.""" self.m_convertFaultToSurface = True def ConvertSurfaceMeshOff( self: Self ) -> None: - """Deactivate surface conversion from vtkUnstructredGrid to vtkPolyData.""" + """Deactivate surface conversion from vtkUnstructuredGrid to vtkPolyData.""" self.m_convertFaultToSurface = False - def doMerge( self: Self ) -> int: - """Apply block merge. - - Returns: - bool: True if block merge successfully ended, False otherwise. - """ - self.mergeRankBlocks() - if self.m_convertFaultToSurface: - self.convertFaultsToSurfaces() - - self.m_output.ShallowCopy( self.m_outputMesh ) - return 1 - - def mergeRankBlocks( self: Self ) -> bool: - """Merge all elementary node that belong to a same parent node. - - Returns: - bool: True if calculation successfully ended, False otherwise - """ - # display phase names + def applyFilter( self: Self ) -> None: + """Merge all elementary node that belong to a same parent node.""" try: - phaseClassification: dict[ str, PhaseTypeEnum ] = self.getPhases( True ) - if phaseClassification is not None: - for phaseTypeRef in list( PhaseTypeEnum ): - phases = [ - phaseName for phaseName, phaseType in phaseClassification.items() if phaseType is phaseTypeRef - ] - if len( phases ) > 0: - self.m_logger.info( f"Identified {phaseTypeRef.type} phase(s) are: {phases}" ) - except AssertionError as e: - self.m_logger.warning( "Phases were not identified due to:" ) - self.m_logger.warning( e ) - - compositeBlockIndexesToMerge: dict[ str, int ] = ( getElementaryCompositeBlockIndexes( self.m_input ) ) - nbBlocks: int = len( compositeBlockIndexesToMerge ) - self.m_outputMesh = vtkMultiBlockDataSet() - self.m_outputMesh.SetNumberOfBlocks( nbBlocks ) - for newIndex, ( blockName, blockIndex ) in enumerate( compositeBlockIndexesToMerge.items() ): - # extract composite block - blockToMerge1: vtkMultiBlockDataSet = extractBlock( self.m_input, blockIndex ) - assert blockToMerge1 is not None, "Extracted block to merge is null." - - # rename attributes - blockToMerge2: vtkMultiBlockDataSet = self.renameAttributes( blockToMerge1, phaseClassification ) - assert blockToMerge2 is not None, "Attribute renaming failed." - - # merge all its children - mergedBlock: vtkUnstructuredGrid = self.mergeChildBlocks( blockToMerge2 ) - - # create index attribute keeping the index in intial mesh - if not createConstantAttribute( - mergedBlock, - [ - blockIndex, - ], - PostProcessingOutputsEnum.BLOCK_INDEX.attributeName, - (), - False, - ): - self.m_logger.warning( "BlockIndex attribute was not created." ) - - # set this composite block into the output - self.m_outputMesh.SetBlock( newIndex, mergedBlock ) - self.m_outputMesh.GetMetaData( newIndex ).Set( vtkCompositeDataSet.NAME(), blockName ) - - assert ( self.m_outputMesh.GetNumberOfBlocks() == nbBlocks ), "Final number of merged blocks is wrong." - return True + # Display phase names + self.commutePhaseNames() + for phase, phaseNames in self.phaseNameDict.items(): + if len( phaseNames ) > 0: + self.m_logger.info( f"Identified { phase } phase(s) are: { phaseNames }." ) + else: + self.m_logger.info( f"No { phase } phase has been identified." ) + + # Merge all the composite blocks + compositeBlockIndexesToMerge: dict[ str, int ] = getElementaryCompositeBlockIndexes( self.m_inputMesh ) + nbBlocks: int = len( compositeBlockIndexesToMerge ) + self.m_outputMesh.SetNumberOfBlocks( nbBlocks ) + for newIndex, ( blockName, blockIndex ) in enumerate( compositeBlockIndexesToMerge.items() ): + # Set the name of the merged block + self.m_outputMesh.GetMetaData( newIndex ).Set( vtkCompositeDataSet.NAME(), blockName ) + + # Merge blocks + blockToMerge: vtkMultiBlockDataSet = vtkMultiBlockDataSet.SafeDownCast( self.m_inputMesh.GetBlock( blockIndex ) ) + volumeMesh: vtkUnstructuredGrid = mergeBlocks( blockToMerge, keepPartialAttributes=True, logger=self.m_logger ) + + # Create index attribute keeping the index in initial mesh + if not createConstantAttribute( volumeMesh, [ blockIndex ], PostProcessingOutputsEnum.BLOCK_INDEX.attributeName, onPoints=False, logger=self.m_logger ): + self.m_logger.warning( "BlockIndex attribute was not created." ) + + # Rename attributes + self.renameAttributes( volumeMesh ) + + # Convert the volume mesh to a surface mesh + if self.m_convertFaultToSurface: + surfaceMesh: vtkPolyData = self.convertBlockToSurface( volumeMesh ) + assert surfaceMesh is not None, "Surface extraction from block failed." + surfaceMesh.ShallowCopy( computeNormals( surfaceMesh, logger=self.m_logger ) ) + assert surfaceMesh is not None, "Normal calculation failed." + surfaceMesh.ShallowCopy( computeTangents( surfaceMesh, logger=self.m_logger ) ) + assert surfaceMesh is not None, "Tangent calculation failed." + # Add the merged block to the output + self.m_outputMesh.SetBlock( newIndex, surfaceMesh ) + else: + self.m_outputMesh.SetBlock( newIndex, volumeMesh ) + except ( ValueError, TypeError, RuntimeError ) as e: + self.m_logger.critical( "Geos block merge failed due to:" ) + self.m_logger.critical( e, exc_info=True ) + + return def renameAttributes( self: Self, - mesh: vtkMultiBlockDataSet, - phaseClassification: dict[ str, PhaseTypeEnum ], - ) -> vtkMultiBlockDataSet: + mesh: vtkUnstructuredGrid, + ) -> None: """Rename attributes to harmonize throughout the mesh. Args: mesh (vtkMultiBlockDataSet): input mesh phaseClassification (dict[str, PhaseTypeEnum]): phase classification detected from attributes - - Returns: - vtkMultiBlockDataSet: output mesh with renamed attributes """ - assert phaseClassification is not None, "Phases were not correctly identified." - - renameFilter: vtkArrayRename = vtkArrayRename() - renameFilter.SetInputData( mesh ) - rockPhase: list[ str ] = [ - phaseName for phaseName, phaseType in phaseClassification.items() if phaseType is PhaseTypeEnum.ROCK - ] for attributeName in getAttributeSet( mesh, False ): - for phaseName in rockPhase: - if phaseName in attributeName: - for suffix, newName in getRockSuffixRenaming().items(): - if suffix in attributeName: - renameFilter.SetCellArrayName( attributeName, newName ) - - renameFilter.Update() - output: vtkMultiBlockDataSet = renameFilter.GetOutput() - return output - - def getPhaseNames( self: Self, attributeSet: set[ str ] ) -> set[ str ]: - """Get the names of the phases in the mesh from Point/Cell attributes. + for suffix, newName in getRockSuffixRenaming().items(): + if suffix in attributeName: + renameAttribute( mesh, attributeName, newName, False ) - Args: - attributeSet (set[str]): list of attributes where to find phases - - Returns: - set[str]: the list of phase names that appear at least twice. - """ - phaseNameDict: dict[ str, int ] = {} - for name in attributeSet: + def commutePhaseNames( self: Self ) -> None: + """Get the names of the phases in the mesh from Cell attributes.""" + for name in getAttributeSet( self.m_inputMesh, False ): if PHASE_SEP in name: - # use the last occurence of PHASE_SEP to split phase name from - # property name index = name.rindex( PHASE_SEP ) phaseName: str = name[ :index ] - if phaseName in phaseNameDict: - phaseNameDict[ phaseName ] += 1 - else: - phaseNameDict[ phaseName ] = 1 - - # remove names that appear only once - return set( phaseNameDict.keys() ) - - def getPhases( self: Self, onCells: bool = True ) -> dict[ str, PhaseTypeEnum ]: - """Get the dictionnary of phases classified according to PhaseTypeEnum. - - Args: - onCells (bool, optional): Attributes are on Cells (Default) or on - Points. + suffixName: str = name[ index: ] + if suffixName in PhaseTypeEnum.ROCK.attributes: + self.phaseNameDict[ PhaseTypeEnum.ROCK.type ].add( phaseName ) + elif suffixName in PhaseTypeEnum.FLUID.attributes: + self.phaseNameDict[ PhaseTypeEnum.FLUID.type ].add( phaseName ) - Defaults to True. - - Returns: - dict[str, PhaseTypeEnum]: a dictionnary with phase names as keys and - phase type as value. - """ - attributeSet: set[ str ] = getAttributeSet( self.m_input, not onCells ) - assert len( attributeSet ) > 0, "Input object does not have any attribute." - - phaseClassification: dict[ str, PhaseTypeEnum ] = {} - phaseNames: set[ str ] = self.getPhaseNames( attributeSet ) - for phaseName in phaseNames: - # check for fluid phase names (most often the same names: fluid, water, or gas) - if any( fluidPrefix.value.lower() in phaseName.lower() for fluidPrefix in list( FluidPrefixEnum ) ): - phaseClassification[ phaseName ] = PhaseTypeEnum.FLUID - continue - - for attributeName in attributeSet: - if phaseName in attributeName: - index = attributeName.index( phaseName ) + len( phaseName ) - suffix = attributeName[ index: ] - for phaseType in list( PhaseTypeEnum ): - if suffix in phaseType.attributes: - if ( phaseName in phaseClassification ) and ( phaseType - is not phaseClassification[ phaseName ] ): - self.m_logger.warning( f"The phase {phaseName} may be misclassified " + - "since the same name is used for both " + - "{phaseType.type} and " + - "{phaseClassification[phaseName].type} types" ) - phaseClassification[ phaseName ] = phaseType - - return phaseClassification - - def mergeChildBlocks( self: Self, compositeBlock: vtkMultiBlockDataSet ) -> vtkUnstructuredGrid: - """Merge all children of the input composite block. - - Args: - compositeBlock (vtkMultiBlockDataSet): composite block - - Returns: - vtkUnstructuredGrid: merged block - """ - # fill partial attributes in all children blocks - if not fillAllPartialAttributes( compositeBlock ): - self.m_logger.warning( "Some partial attributes may not have been propagated to the whole mesh." ) - - # merge blocks - mergedBlocks = mergeBlocks( compositeBlock ) - return mergedBlocks - - def convertFaultsToSurfaces( self: Self ) -> bool: - """Convert blocks corresponding to faults to surface. - - Returns: - bool: True if calculation succesfully ended, False otherwise - """ - assert self.m_outputMesh is not None, "Output surface is null." - - transientMesh: vtkMultiBlockDataSet = vtkMultiBlockDataSet() - nbSurfaces: int = self.m_outputMesh.GetNumberOfBlocks() - transientMesh.SetNumberOfBlocks( nbSurfaces ) - - # initialize data object tree iterator - iter: vtkDataObjectTreeIterator = vtkDataObjectTreeIterator() - iter.SetDataSet( self.m_outputMesh ) - iter.VisitOnlyLeavesOn() - iter.GoToFirstItem() - surfaceIndex: int = 0 - while iter.GetCurrentDataObject() is not None: - surfaceName: str = iter.GetCurrentMetaData().Get( vtkCompositeDataSet.NAME() ) - # convert block to surface - surface0: vtkUnstructuredGrid = vtkUnstructuredGrid.SafeDownCast( iter.GetCurrentDataObject() ) - surface1: vtkPolyData = self.convertBlockToSurface( surface0 ) - assert surface1 is not None, "Surface extraction from block failed." - # compute normals - surface2: vtkPolyData = computeNormals( surface1 ) - assert surface2 is not None, "Normal calculation failed." - # compute tangents - surface3: vtkPolyData = computeTangents( surface2 ) - assert surface3 is not None, "Tangent calculation failed." - - # set surface to output multiBlockDataSet - transientMesh.SetBlock( surfaceIndex, surface3 ) - transientMesh.GetMetaData( surfaceIndex ).Set( vtkCompositeDataSet.NAME(), surfaceName ) - - iter.GoToNextItem() - surfaceIndex += 1 - - self.m_outputMesh.ShallowCopy( transientMesh ) - return True + return def convertBlockToSurface( self: Self, block: vtkUnstructuredGrid ) -> vtkPolyData: """Convert vtkUnstructuredGrid to a surface vtkPolyData. From c2ab8fb9ef4116bfd0ccc97fac7169581fd1357e Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Mon, 3 Nov 2025 16:56:46 +0100 Subject: [PATCH 27/34] Update import --- geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py | 2 +- .../src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurface.py | 2 +- .../geos/pv/plugins/PVExtractMergeBlocksVolumeSurfaceWell.py | 2 +- geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeWell.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py index ef97c659..b495e747 100644 --- a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py +++ b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py @@ -23,7 +23,7 @@ ) from geos.utils.Logger import ERROR, INFO, Logger, getLogger from geos.processing.post_processing.GeosBlockExtractor import GeosBlockExtractor -from geos_posp.filters.GeosBlockMerge import GeosBlockMerge +from geos.processing.post_processing.GeosBlockMerge import GeosBlockMerge from geos.mesh.utils.arrayModifiers import ( copyAttribute, createCellCenterAttribute, diff --git a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurface.py b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurface.py index c9a27e53..33ac4f29 100644 --- a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurface.py +++ b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurface.py @@ -24,7 +24,7 @@ ) from geos.utils.Logger import ERROR, INFO, Logger, getLogger from geos.processing.post_processing.GeosBlockExtractor import GeosBlockExtractor -from geos_posp.filters.GeosBlockMerge import GeosBlockMerge +from geos.processing.post_processing.GeosBlockMerge import GeosBlockMerge from geos.mesh.utils.arrayModifiers import ( copyAttribute, createCellCenterAttribute, diff --git a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurfaceWell.py b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurfaceWell.py index 74710a34..ec388e48 100644 --- a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurfaceWell.py +++ b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurfaceWell.py @@ -24,7 +24,7 @@ ) from geos.utils.Logger import ERROR, INFO, Logger, getLogger from geos.processing.post_processing.GeosBlockExtractor import GeosBlockExtractor -from geos_posp.filters.GeosBlockMerge import GeosBlockMerge +from geos.processing.post_processing.GeosBlockMerge import GeosBlockMerge from geos.mesh.utils.arrayModifiers import ( copyAttribute, createCellCenterAttribute, diff --git a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeWell.py b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeWell.py index 4b7e4389..9ae83c7f 100644 --- a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeWell.py +++ b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeWell.py @@ -27,7 +27,7 @@ ) from geos.utils.Logger import ERROR, INFO, Logger, getLogger from geos.processing.post_processing.GeosBlockExtractor import GeosBlockExtractor -from geos_posp.filters.GeosBlockMerge import GeosBlockMerge +from geos.processing.post_processing.GeosBlockMerge import GeosBlockMerge from geos.mesh.utils.arrayModifiers import ( copyAttribute, createCellCenterAttribute, From c2df04e03aa2c85f03972c94680356d5df1f9543 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 4 Nov 2025 08:41:50 +0100 Subject: [PATCH 28/34] Fix density name for fluid phase --- .../post_processing/GeosBlockMerge.py | 27 ++++++++++++++----- .../pv/plugins/PVExtractMergeBlocksVolume.py | 11 +++----- .../src/geos/utils/GeosOutputsConstants.py | 4 +-- 3 files changed, 27 insertions(+), 15 deletions(-) diff --git a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py index f20d815d..bdc34e3f 100644 --- a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py +++ b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py @@ -6,6 +6,7 @@ from geos.utils.GeosOutputsConstants import ( PHASE_SEP, PhaseTypeEnum, + FluidPrefixEnum, PostProcessingOutputsEnum, getRockSuffixRenaming, ) @@ -19,7 +20,7 @@ ) from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter -from geos.mesh.utils.multiblockHelpers import getElementaryCompositeBlockIndexes +from geos.mesh.utils.multiblockHelpers import getElementaryCompositeBlockIndexes, extractBlock from geos.mesh.utils.arrayHelpers import getAttributeSet from geos.mesh.utils.arrayModifiers import createConstantAttribute, renameAttribute from geos.mesh.utils.multiblockModifiers import mergeBlocks @@ -73,8 +74,8 @@ def __init__( self: Self, inputMesh: vtkMultiBlockDataSet ) -> None: self.m_convertFaultToSurface: bool = False self.phaseNameDict: dict[ str, set[ str ] ] = { - PhaseTypeEnum.ROCK.type: {}, - PhaseTypeEnum.FLUID.type: {}, + PhaseTypeEnum.ROCK.type: set(), + PhaseTypeEnum.FLUID.type: set(), } # set logger @@ -96,6 +97,10 @@ def ConvertSurfaceMeshOff( self: Self ) -> None: """Deactivate surface conversion from vtkUnstructuredGrid to vtkPolyData.""" self.m_convertFaultToSurface = False + def getOutput ( self: Self ) -> vtkMultiBlockDataSet: + """Get the mesh with the composite blocks merged.""" + return self.m_outputMesh + def applyFilter( self: Self ) -> None: """Merge all elementary node that belong to a same parent node.""" try: @@ -116,7 +121,7 @@ def applyFilter( self: Self ) -> None: self.m_outputMesh.GetMetaData( newIndex ).Set( vtkCompositeDataSet.NAME(), blockName ) # Merge blocks - blockToMerge: vtkMultiBlockDataSet = vtkMultiBlockDataSet.SafeDownCast( self.m_inputMesh.GetBlock( blockIndex ) ) + blockToMerge: vtkMultiBlockDataSet = extractBlock( self.m_inputMesh, blockIndex ) volumeMesh: vtkUnstructuredGrid = mergeBlocks( blockToMerge, keepPartialAttributes=True, logger=self.m_logger ) # Create index attribute keeping the index in initial mesh @@ -158,7 +163,12 @@ def renameAttributes( for attributeName in getAttributeSet( mesh, False ): for suffix, newName in getRockSuffixRenaming().items(): if suffix in attributeName: - renameAttribute( mesh, attributeName, newName, False ) + if suffix == "_density": + for phaseName in self.phaseNameDict[ PhaseTypeEnum.ROCK.type ]: + if phaseName in attributeName: + renameAttribute( mesh, attributeName, newName, False ) + elif suffix != "_density": + renameAttribute( mesh, attributeName, newName, False ) def commutePhaseNames( self: Self ) -> None: """Get the names of the phases in the mesh from Cell attributes.""" @@ -167,7 +177,12 @@ def commutePhaseNames( self: Self ) -> None: index = name.rindex( PHASE_SEP ) phaseName: str = name[ :index ] suffixName: str = name[ index: ] - if suffixName in PhaseTypeEnum.ROCK.attributes: + if suffixName == "_density": + if any( phaseName in fluidPrefix.value for fluidPrefix in list( FluidPrefixEnum ) ): + self.phaseNameDict[ PhaseTypeEnum.FLUID.type ].add( phaseName ) + else: + self.phaseNameDict[ PhaseTypeEnum.ROCK.type ].add( phaseName ) + elif suffixName in PhaseTypeEnum.ROCK.attributes: self.phaseNameDict[ PhaseTypeEnum.ROCK.type ].add( phaseName ) elif suffixName in PhaseTypeEnum.FLUID.attributes: self.phaseNameDict[ PhaseTypeEnum.FLUID.type ].add( phaseName ) diff --git a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py index b495e747..c646f5f0 100644 --- a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py +++ b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py @@ -274,7 +274,7 @@ def doExtractAndMerge( self: Self, input: vtkMultiBlockDataSet, output: vtkMulti output (vtkMultiBlockDataSet): output volume mesh Returns: - bool: True if extraction and merge successfully eneded, False otherwise + bool: True if extraction and merge successfully ended, False otherwise """ # extract blocks blockExtractor: GeosBlockExtractor = GeosBlockExtractor( input ) @@ -303,14 +303,11 @@ def mergeBlocksFilter( self: Self, Returns: vtkMultiBlockDataSet: Multiblock mesh composed of internal merged blocks. """ - mergeBlockFilter: GeosBlockMerge = GeosBlockMerge() + mergeBlockFilter: GeosBlockMerge = GeosBlockMerge( input ) mergeBlockFilter.SetLogger( self.m_logger ) - mergeBlockFilter.SetInputDataObject( input ) if convertSurfaces: mergeBlockFilter.ConvertSurfaceMeshOn() - else: - mergeBlockFilter.ConvertSurfaceMeshOff() - mergeBlockFilter.Update() - mergedBlocks: vtkMultiBlockDataSet = mergeBlockFilter.GetOutputDataObject( 0 ) + mergeBlockFilter.applyFilter() + mergedBlocks: vtkMultiBlockDataSet = mergeBlockFilter.getOutput() assert mergedBlocks is not None, "Final merged MultiBlockDataSet is null." return mergedBlocks diff --git a/geos-utils/src/geos/utils/GeosOutputsConstants.py b/geos-utils/src/geos/utils/GeosOutputsConstants.py index a2fc72cb..52d50a3f 100644 --- a/geos-utils/src/geos/utils/GeosOutputsConstants.py +++ b/geos-utils/src/geos/utils/GeosOutputsConstants.py @@ -102,9 +102,9 @@ class GeosMeshSuffixEnum( Enum ): BIOT_COEFFICIENT_SUFFIX = "_biotCoefficient" # fluid attributes suffix - PHASE_DENSITY_SUFFIX = "_phaseDensity" + PHASE_DENSITY_SUFFIX = "_density" PHASE_MASS_DENSITY_SUFFIX = "_phaseMassDensity" - PHASE_VISCOSITY_SUFFIX = "_phaseViscosity" + PHASE_VISCOSITY_SUFFIX = "_viscosity" PHASE_FRACTION_SUFFIX = "_phaseFraction" # surface attribute transfer suffix From eabc05805b78a3d38c99e3b09d2c6b5ed391f957 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 4 Nov 2025 15:28:10 +0100 Subject: [PATCH 29/34] harmonized the code (logger, varriable name ...) --- .../post_processing/GeosBlockMerge.py | 126 +++++++++++------- 1 file changed, 76 insertions(+), 50 deletions(-) diff --git a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py index bdc34e3f..af4c2f10 100644 --- a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py +++ b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py @@ -2,7 +2,7 @@ # # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. # SPDX-FileContributor: Martin Lemay, Romain Baville # ruff: noqa: E402 # disable Module level import not at top of file - +import logging from geos.utils.GeosOutputsConstants import ( PHASE_SEP, PhaseTypeEnum, @@ -17,6 +17,8 @@ vtkMultiBlockDataSet, vtkPolyData, vtkUnstructuredGrid, + vtkCellTypes, + vtkTriangle, ) from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter @@ -60,92 +62,111 @@ output :vtkMultiBlockDataSet = mergeBlockFilter.GetOutputDataObject(0) """ +loggerTitle: str = "GEOS Block Merge" + class GeosBlockMerge(): - def __init__( self: Self, inputMesh: vtkMultiBlockDataSet ) -> None: + def __init__( + self: Self, + inputMesh: vtkMultiBlockDataSet, + convertFaultToSurface: bool = False, + speHandler: bool = False, + ) -> None: """VTK Filter that perform GEOS rank merge. - The filter returns a multiblock mesh composed of elementary blocks. + Args: + inputMesh (vtkMultiBlockDataSet): The mesh with the blocks to merge. + convertFaultToSurface (bool, optional): True if the merged block need to be convert to vtp, False otherwise. + Defaults to False. + speHandler (bool, optional): True to use a specific handler, False to use the internal handler. + Defaults to False. """ - self.m_inputMesh: vtkMultiBlockDataSet = inputMesh - self.m_outputMesh: vtkMultiBlockDataSet = vtkMultiBlockDataSet() + self.inputMesh: vtkMultiBlockDataSet = inputMesh + self.convertFaultToSurface: bool = convertFaultToSurface - self.m_convertFaultToSurface: bool = False + self.outputMesh: vtkMultiBlockDataSet = vtkMultiBlockDataSet() self.phaseNameDict: dict[ str, set[ str ] ] = { PhaseTypeEnum.ROCK.type: set(), PhaseTypeEnum.FLUID.type: set(), } - # set logger - self.m_logger: Logger = getLogger( "Geos Block Merge Filter" ) + # Logger. + self.logger: Logger + if not speHandler: + self.logger = getLogger( loggerTitle, True ) + else: + self.logger = logging.getLogger( loggerTitle ) + self.logger.setLevel( logging.INFO ) + + def setLoggerHandler( self: Self, handler: logging.Handler ) -> None: + """Set a specific handler for the filter logger. - def SetLogger( self: Self, logger: Logger ) -> None: - """Set the 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: - logger (Logger): logger + handler (logging.Handler): The handler to add. """ - self.m_logger = logger - - def ConvertSurfaceMeshOn( self: Self ) -> None: - """Activate surface conversion from vtkUnstructuredGrid to vtkPolyData.""" - self.m_convertFaultToSurface = True - - def ConvertSurfaceMeshOff( self: Self ) -> None: - """Deactivate surface conversion from vtkUnstructuredGrid to vtkPolyData.""" - self.m_convertFaultToSurface = False + 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 getOutput ( self: Self ) -> vtkMultiBlockDataSet: """Get the mesh with the composite blocks merged.""" - return self.m_outputMesh + return self.outputMesh def applyFilter( self: Self ) -> None: """Merge all elementary node that belong to a same parent node.""" + self.logger.info( f"Apply filter { self.logger.name }." ) + try: # Display phase names self.commutePhaseNames() for phase, phaseNames in self.phaseNameDict.items(): if len( phaseNames ) > 0: - self.m_logger.info( f"Identified { phase } phase(s) are: { phaseNames }." ) + self.logger.info( f"Identified { phase } phase(s) are: { phaseNames }." ) else: - self.m_logger.info( f"No { phase } phase has been identified." ) + self.logger.info( f"No { phase } phase has been identified." ) - # Merge all the composite blocks - compositeBlockIndexesToMerge: dict[ str, int ] = getElementaryCompositeBlockIndexes( self.m_inputMesh ) + # Parse all the composite blocks + compositeBlockIndexesToMerge: dict[ str, int ] = getElementaryCompositeBlockIndexes( self.inputMesh ) nbBlocks: int = len( compositeBlockIndexesToMerge ) - self.m_outputMesh.SetNumberOfBlocks( nbBlocks ) + self.outputMesh.SetNumberOfBlocks( nbBlocks ) for newIndex, ( blockName, blockIndex ) in enumerate( compositeBlockIndexesToMerge.items() ): - # Set the name of the merged block - self.m_outputMesh.GetMetaData( newIndex ).Set( vtkCompositeDataSet.NAME(), blockName ) + # Set the name of the composite block + self.outputMesh.GetMetaData( newIndex ).Set( vtkCompositeDataSet.NAME(), blockName ) # Merge blocks - blockToMerge: vtkMultiBlockDataSet = extractBlock( self.m_inputMesh, blockIndex ) - volumeMesh: vtkUnstructuredGrid = mergeBlocks( blockToMerge, keepPartialAttributes=True, logger=self.m_logger ) + blockToMerge: vtkMultiBlockDataSet = extractBlock( self.inputMesh, blockIndex ) + volumeMesh: vtkUnstructuredGrid = mergeBlocks( blockToMerge, keepPartialAttributes=True, logger=self.logger ) # Create index attribute keeping the index in initial mesh - if not createConstantAttribute( volumeMesh, [ blockIndex ], PostProcessingOutputsEnum.BLOCK_INDEX.attributeName, onPoints=False, logger=self.m_logger ): - self.m_logger.warning( "BlockIndex attribute was not created." ) + if not createConstantAttribute( volumeMesh, [ blockIndex ], PostProcessingOutputsEnum.BLOCK_INDEX.attributeName, onPoints=False, logger=self.logger ): + self.logger.warning( "BlockIndex attribute was not created." ) # Rename attributes self.renameAttributes( volumeMesh ) # Convert the volume mesh to a surface mesh - if self.m_convertFaultToSurface: + if self.convertFaultToSurface: surfaceMesh: vtkPolyData = self.convertBlockToSurface( volumeMesh ) assert surfaceMesh is not None, "Surface extraction from block failed." - surfaceMesh.ShallowCopy( computeNormals( surfaceMesh, logger=self.m_logger ) ) + surfaceMesh.ShallowCopy( computeNormals( surfaceMesh, logger=self.logger ) ) assert surfaceMesh is not None, "Normal calculation failed." - surfaceMesh.ShallowCopy( computeTangents( surfaceMesh, logger=self.m_logger ) ) + surfaceMesh.ShallowCopy( computeTangents( surfaceMesh, logger=self.logger ) ) assert surfaceMesh is not None, "Tangent calculation failed." - # Add the merged block to the output - self.m_outputMesh.SetBlock( newIndex, surfaceMesh ) + # Add the merged block to the output mesh + self.outputMesh.SetBlock( newIndex, surfaceMesh ) else: - self.m_outputMesh.SetBlock( newIndex, volumeMesh ) + self.outputMesh.SetBlock( newIndex, volumeMesh ) + + self.logger.info( "The filter succeeded." ) except ( ValueError, TypeError, RuntimeError ) as e: - self.m_logger.critical( "Geos block merge failed due to:" ) - self.m_logger.critical( e, exc_info=True ) + self.logger.critical( f"The filter failed.\n{ e }" ) return @@ -153,31 +174,33 @@ def renameAttributes( self: Self, mesh: vtkUnstructuredGrid, ) -> None: - """Rename attributes to harmonize throughout the mesh. + """Rename attributes to harmonize GEOS output, see more geos.utils.OutputsConstants.py. Args: - mesh (vtkMultiBlockDataSet): input mesh - phaseClassification (dict[str, PhaseTypeEnum]): phase classification - detected from attributes + mesh (vtkUnstructuredGrid): The mesh with the attribute to rename. """ + # All the attributes to rename are on cells for attributeName in getAttributeSet( mesh, False ): for suffix, newName in getRockSuffixRenaming().items(): if suffix in attributeName: + # Fluid and Rock density attribute have the same suffix, only the rock density need to be renamed if suffix == "_density": for phaseName in self.phaseNameDict[ PhaseTypeEnum.ROCK.type ]: if phaseName in attributeName: renameAttribute( mesh, attributeName, newName, False ) - elif suffix != "_density": + else: renameAttribute( mesh, attributeName, newName, False ) def commutePhaseNames( self: Self ) -> None: """Get the names of the phases in the mesh from Cell attributes.""" - for name in getAttributeSet( self.m_inputMesh, False ): + # All the phase attributes are on cells + for name in getAttributeSet( self.inputMesh, False ): if PHASE_SEP in name: - index = name.rindex( PHASE_SEP ) - phaseName: str = name[ :index ] - suffixName: str = name[ index: ] - if suffixName == "_density": + phaseName: str + suffixName: str + phaseName, suffixName = name.split( PHASE_SEP ) + # Fluid and Rock density attribute have the same suffix, common fluid name are used to separated the two phases + if suffixName == "density": if any( phaseName in fluidPrefix.value for fluidPrefix in list( FluidPrefixEnum ) ): self.phaseNameDict[ PhaseTypeEnum.FLUID.type ].add( phaseName ) else: @@ -203,6 +226,9 @@ def convertBlockToSurface( self: Self, block: vtkUnstructuredGrid ) -> vtkPolyDa Returns: vtkPolyData: extracted surface """ + cellTypes: list[ vtkCellTypes ] = block.GetDistinctCellTypesArray() + assert [ vtkTriangle ] == cellTypes, "The surface mesh must be a triangulated surface only." + extractSurfaceFilter: vtkDataSetSurfaceFilter = vtkDataSetSurfaceFilter() extractSurfaceFilter.SetInputData( block ) # fast mode should be used for rendering only From ab7bfcbd05d9300eb6740ce4ef626d47064bf017 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 4 Nov 2025 15:52:30 +0100 Subject: [PATCH 30/34] Clean the doc --- .../post_processing/GeosBlockMerge.py | 102 +++++++++--------- 1 file changed, 51 insertions(+), 51 deletions(-) diff --git a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py index af4c2f10..f3203b8d 100644 --- a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py +++ b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py @@ -3,63 +3,54 @@ # SPDX-FileContributor: Martin Lemay, Romain Baville # ruff: noqa: E402 # disable Module level import not at top of file import logging -from geos.utils.GeosOutputsConstants import ( - PHASE_SEP, - PhaseTypeEnum, - FluidPrefixEnum, - PostProcessingOutputsEnum, - getRockSuffixRenaming, -) -from geos.utils.Logger import Logger, getLogger from typing_extensions import Self -from vtkmodules.vtkCommonDataModel import ( - vtkCompositeDataSet, - vtkMultiBlockDataSet, - vtkPolyData, - vtkUnstructuredGrid, - vtkCellTypes, - vtkTriangle, -) + +from vtkmodules.vtkCommonDataModel import ( vtkCompositeDataSet, vtkMultiBlockDataSet, vtkPolyData, vtkUnstructuredGrid, + vtkCellTypes, vtkTriangle ) from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter -from geos.mesh.utils.multiblockHelpers import getElementaryCompositeBlockIndexes, extractBlock +from geos.utils.Logger import ( Logger, getLogger ) +from geos.utils.GeosOutputsConstants import ( PHASE_SEP, PhaseTypeEnum, FluidPrefixEnum, PostProcessingOutputsEnum, + getRockSuffixRenaming ) from geos.mesh.utils.arrayHelpers import getAttributeSet -from geos.mesh.utils.arrayModifiers import createConstantAttribute, renameAttribute +from geos.mesh.utils.arrayModifiers import ( createConstantAttribute, renameAttribute ) from geos.mesh.utils.multiblockModifiers import mergeBlocks -from geos.mesh.utils.genericHelpers import ( - computeNormals, - computeTangents, -) +from geos.mesh.utils.multiblockHelpers import ( getElementaryCompositeBlockIndexes, extractBlock ) +from geos.mesh.utils.genericHelpers import ( computeNormals, computeTangents ) __doc__ = """ -GeosBlockMerge module is a vtk filter that allows to merge Geos ranks, rename -stress and porosity attributes, and identify fluids and rock phases. +GeosBlockMerge is a vtk filter that allows to merge for a GEOS domain the ranks per region, identify "Fluids" and "Rock" phases and rename "Rock" attributes. Filter input and output types are vtkMultiBlockDataSet. +.. Important:: + This filter deals with the domain mesh of GEOS. This domain needs to be extracted before. + See geos-processing/src/geos/processing/post_processing/GeosBlockExtractor.py to see the type of input requires by this filter. + To use the filter: .. code-block:: python - from filters.GeosBlockMerge import GeosBlockMerge - - # filter inputs - logger :Logger - input :vtkMultiBlockDataSet - - # instantiate the filter - mergeBlockFilter :GeosBlockMerge = GeosBlockMerge() - # set the logger - mergeBlockFilter.SetLogger(logger) - # set input data object - mergeBlockFilter.SetInputDataObject(input) - # ConvertSurfaceMeshOff or ConvertSurfaceMeshOn to (de)activate the conversion - # of vtkUnstructuredGrid to surface withNormals and Tangents calculation. - mergeBlockFilter.ConvertSurfaceMeshOff() - # do calculations - mergeBlockFilter.Update() - # get output object - output :vtkMultiBlockDataSet = mergeBlockFilter.GetOutputDataObject(0) + from geos.processing.post_processing.GeosBlockMerge import GeosBlockMerge + + # Filter inputs. + inputMesh: vtkMultiBlockDataSet + # Optional inputs. + convertFaultToSurface: bool # Defaults to False + speHandler: bool # Defaults to False + + # Instantiate the filter + mergeBlockFilter: GeosBlockMerge = GeosBlockMerge( inputMesh, convertFaultToSurface, speHandler ) + + # Set the handler of yours (only if speHandler is True). + yourHandler: logging.Handler + mergeBlockFilter.setLoggerHandler( yourHandler ) + + # Do calculations + mergeBlockFilter.applyFilter() + + # Get the multiBlockDataSet with one dataSet per region + outputMesh: vtkMultiBlockDataSet = mergeBlockFilter.getOutput() """ loggerTitle: str = "GEOS Block Merge" @@ -73,7 +64,12 @@ def __init__( convertFaultToSurface: bool = False, speHandler: bool = False, ) -> None: - """VTK Filter that perform GEOS rank merge. + """VTK Filter that merge ranks of GEOS output mesh. + + for all the composite blocks of the input mesh: + - Ranks are merged + - "Rock" attributes are renamed + - Volume mesh are convert to surface if needed Args: inputMesh (vtkMultiBlockDataSet): The mesh with the blocks to merge. @@ -115,12 +111,12 @@ def setLoggerHandler( self: Self, handler: logging.Handler ) -> None: "The logger already has an handler, to use yours set the argument 'speHandler' to True during the filter initialization." ) - def getOutput ( self: Self ) -> vtkMultiBlockDataSet: + def getOutput( self: Self ) -> vtkMultiBlockDataSet: """Get the mesh with the composite blocks merged.""" return self.outputMesh def applyFilter( self: Self ) -> None: - """Merge all elementary node that belong to a same parent node.""" + """Apply the filter on the mesh.""" self.logger.info( f"Apply filter { self.logger.name }." ) try: @@ -142,10 +138,15 @@ def applyFilter( self: Self ) -> None: # Merge blocks blockToMerge: vtkMultiBlockDataSet = extractBlock( self.inputMesh, blockIndex ) - volumeMesh: vtkUnstructuredGrid = mergeBlocks( blockToMerge, keepPartialAttributes=True, logger=self.logger ) + volumeMesh: vtkUnstructuredGrid = mergeBlocks( blockToMerge, + keepPartialAttributes=True, + logger=self.logger ) # Create index attribute keeping the index in initial mesh - if not createConstantAttribute( volumeMesh, [ blockIndex ], PostProcessingOutputsEnum.BLOCK_INDEX.attributeName, onPoints=False, logger=self.logger ): + if not createConstantAttribute( volumeMesh, [ blockIndex ], + PostProcessingOutputsEnum.BLOCK_INDEX.attributeName, + onPoints=False, + logger=self.logger ): self.logger.warning( "BlockIndex attribute was not created." ) # Rename attributes @@ -159,13 +160,13 @@ def applyFilter( self: Self ) -> None: assert surfaceMesh is not None, "Normal calculation failed." surfaceMesh.ShallowCopy( computeTangents( surfaceMesh, logger=self.logger ) ) assert surfaceMesh is not None, "Tangent calculation failed." - # Add the merged block to the output mesh + # Add the merged block to the output mesh self.outputMesh.SetBlock( newIndex, surfaceMesh ) else: self.outputMesh.SetBlock( newIndex, volumeMesh ) self.logger.info( "The filter succeeded." ) - except ( ValueError, TypeError, RuntimeError ) as e: + except ( ValueError, TypeError, RuntimeError, AssertionError ) as e: self.logger.critical( f"The filter failed.\n{ e }" ) return @@ -219,7 +220,6 @@ def convertBlockToSurface( self: Self, block: vtkUnstructuredGrid ) -> vtkPolyDa .. TODO:: need to convert quadrangular to triangulated surfaces first - Args: block (vtkUnstructuredGrid): block from which to extract the surface From 789cd31e0a137a647b37b551bb7a497af6d3340e Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 4 Nov 2025 15:58:24 +0100 Subject: [PATCH 31/34] Update GeosBlockMerge calling --- .../pv/plugins/PVExtractMergeBlocksVolume.py | 10 ++++++---- .../PVExtractMergeBlocksVolumeSurface.py | 17 ++++++++--------- .../PVExtractMergeBlocksVolumeSurfaceWell.py | 17 ++++++++--------- .../plugins/PVExtractMergeBlocksVolumeWell.py | 17 ++++++++--------- 4 files changed, 30 insertions(+), 31 deletions(-) diff --git a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py index c646f5f0..382f8bc8 100644 --- a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py +++ b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolume.py @@ -32,6 +32,9 @@ 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, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py __doc__ = """ PVExtractMergeBlocksVolume is a Paraview plugin that allows to merge ranks @@ -303,10 +306,9 @@ def mergeBlocksFilter( self: Self, Returns: vtkMultiBlockDataSet: Multiblock mesh composed of internal merged blocks. """ - mergeBlockFilter: GeosBlockMerge = GeosBlockMerge( input ) - mergeBlockFilter.SetLogger( self.m_logger ) - if convertSurfaces: - mergeBlockFilter.ConvertSurfaceMeshOn() + mergeBlockFilter: GeosBlockMerge = GeosBlockMerge( input, convertSurfaces, True ) + if not mergeBlockFilter.logger.hasHandlers(): + mergeBlockFilter.setLoggerHandler( VTKHandler() ) mergeBlockFilter.applyFilter() mergedBlocks: vtkMultiBlockDataSet = mergeBlockFilter.getOutput() assert mergedBlocks is not None, "Final merged MultiBlockDataSet is null." diff --git a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurface.py b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurface.py index 33ac4f29..a670df2b 100644 --- a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurface.py +++ b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurface.py @@ -33,6 +33,9 @@ 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, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py __doc__ = """ PVExtractMergeBlocksVolumeSurface is a Paraview plugin that allows to merge @@ -325,14 +328,10 @@ def mergeBlocksFilter( self: Self, Returns: vtkMultiBlockDataSet: Multiblock mesh composed of internal merged blocks. """ - mergeBlockFilter = GeosBlockMerge() - mergeBlockFilter.SetLogger( self.m_logger ) - mergeBlockFilter.SetInputDataObject( input ) - if convertSurfaces: - mergeBlockFilter.ConvertSurfaceMeshOn() - else: - mergeBlockFilter.ConvertSurfaceMeshOff() - mergeBlockFilter.Update() - mergedBlocks: vtkMultiBlockDataSet = mergeBlockFilter.GetOutputDataObject( 0 ) + mergeBlockFilter: GeosBlockMerge = GeosBlockMerge( input, convertSurfaces, True ) + if not mergeBlockFilter.logger.hasHandlers(): + mergeBlockFilter.setLoggerHandler( VTKHandler() ) + mergeBlockFilter.applyFilter() + mergedBlocks: vtkMultiBlockDataSet = mergeBlockFilter.getOutput() assert mergedBlocks is not None, "Final merged MultiBlockDataSet is null." return mergedBlocks diff --git a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurfaceWell.py b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurfaceWell.py index ec388e48..d3eff7af 100644 --- a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurfaceWell.py +++ b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeSurfaceWell.py @@ -33,6 +33,9 @@ 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, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py __doc__ = """ PVExtractMergeBlocksVolumeSurfaceWell is a Paraview plugin that allows to merge @@ -353,14 +356,10 @@ def mergeBlocksFilter( self: Self, Returns: vtkMultiBlockDataSet: Multiblock mesh composed of internal merged blocks. """ - mergeBlockFilter = GeosBlockMerge() - mergeBlockFilter.SetLogger( self.m_logger ) - mergeBlockFilter.SetInputDataObject( input ) - if convertSurfaces: - mergeBlockFilter.ConvertSurfaceMeshOn() - else: - mergeBlockFilter.ConvertSurfaceMeshOff() - mergeBlockFilter.Update() - mergedBlocks: vtkMultiBlockDataSet = mergeBlockFilter.GetOutputDataObject( 0 ) + mergeBlockFilter: GeosBlockMerge = GeosBlockMerge( input, convertSurfaces, True ) + if not mergeBlockFilter.logger.hasHandlers(): + mergeBlockFilter.setLoggerHandler( VTKHandler() ) + mergeBlockFilter.applyFilter() + mergedBlocks: vtkMultiBlockDataSet = mergeBlockFilter.getOutput() assert mergedBlocks is not None, "Final merged MultiBlockDataSet is null." return mergedBlocks diff --git a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeWell.py b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeWell.py index 9ae83c7f..d1da166b 100644 --- a/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeWell.py +++ b/geos-pv/src/geos/pv/plugins/PVExtractMergeBlocksVolumeWell.py @@ -10,6 +10,9 @@ 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, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py from typing_extensions import Self from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet @@ -334,14 +337,10 @@ def mergeBlocksFilter( self: Self, Returns: vtkMultiBlockDataSet: Multiblock mesh composed of internal merged blocks. """ - mergeBlockFilter = GeosBlockMerge() - mergeBlockFilter.SetLogger( self.m_logger ) - mergeBlockFilter.SetInputDataObject( input ) - if convertSurfaces: - mergeBlockFilter.ConvertSurfaceMeshOn() - else: - mergeBlockFilter.ConvertSurfaceMeshOff() - mergeBlockFilter.Update() - mergedBlocks: vtkMultiBlockDataSet = mergeBlockFilter.GetOutputDataObject( 0 ) + mergeBlockFilter: GeosBlockMerge = GeosBlockMerge( input, convertSurfaces, True ) + if not mergeBlockFilter.logger.hasHandlers(): + mergeBlockFilter.setLoggerHandler( VTKHandler() ) + mergeBlockFilter.applyFilter() + mergedBlocks: vtkMultiBlockDataSet = mergeBlockFilter.getOutput() assert mergedBlocks is not None, "Final merged MultiBlockDataSet is null." return mergedBlocks From 8f08be87d27ee6d5a2fd09bbd67bed7761d58114 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 4 Nov 2025 16:14:30 +0100 Subject: [PATCH 32/34] fix computePhaseNames --- .../geos/processing/post_processing/GeosBlockMerge.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py index f3203b8d..cf285e29 100644 --- a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py +++ b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py @@ -121,7 +121,7 @@ def applyFilter( self: Self ) -> None: try: # Display phase names - self.commutePhaseNames() + self.computePhaseNames() for phase, phaseNames in self.phaseNameDict.items(): if len( phaseNames ) > 0: self.logger.info( f"Identified { phase } phase(s) are: { phaseNames }." ) @@ -192,7 +192,7 @@ def renameAttributes( else: renameAttribute( mesh, attributeName, newName, False ) - def commutePhaseNames( self: Self ) -> None: + def computePhaseNames( self: Self ) -> None: """Get the names of the phases in the mesh from Cell attributes.""" # All the phase attributes are on cells for name in getAttributeSet( self.inputMesh, False ): @@ -201,14 +201,14 @@ def commutePhaseNames( self: Self ) -> None: suffixName: str phaseName, suffixName = name.split( PHASE_SEP ) # Fluid and Rock density attribute have the same suffix, common fluid name are used to separated the two phases - if suffixName == "density": + if f"{ PHASE_SEP }{ suffixName }" == "_density": if any( phaseName in fluidPrefix.value for fluidPrefix in list( FluidPrefixEnum ) ): self.phaseNameDict[ PhaseTypeEnum.FLUID.type ].add( phaseName ) else: self.phaseNameDict[ PhaseTypeEnum.ROCK.type ].add( phaseName ) - elif suffixName in PhaseTypeEnum.ROCK.attributes: + elif f"{ PHASE_SEP }{ suffixName }" in PhaseTypeEnum.ROCK.attributes: self.phaseNameDict[ PhaseTypeEnum.ROCK.type ].add( phaseName ) - elif suffixName in PhaseTypeEnum.FLUID.attributes: + elif f"{ PHASE_SEP }{ suffixName }" in PhaseTypeEnum.FLUID.attributes: self.phaseNameDict[ PhaseTypeEnum.FLUID.type ].add( phaseName ) return From 13ea5b5a0256d8015bc710a4e49040f5a0e675e2 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 4 Nov 2025 18:40:59 +0100 Subject: [PATCH 33/34] fix convertBlockToSurface --- .../src/geos/processing/post_processing/GeosBlockMerge.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py index cf285e29..53588c6b 100644 --- a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py +++ b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py @@ -226,9 +226,6 @@ def convertBlockToSurface( self: Self, block: vtkUnstructuredGrid ) -> vtkPolyDa Returns: vtkPolyData: extracted surface """ - cellTypes: list[ vtkCellTypes ] = block.GetDistinctCellTypesArray() - assert [ vtkTriangle ] == cellTypes, "The surface mesh must be a triangulated surface only." - extractSurfaceFilter: vtkDataSetSurfaceFilter = vtkDataSetSurfaceFilter() extractSurfaceFilter.SetInputData( block ) # fast mode should be used for rendering only From 77556393f4c3c00e63276e7528f73a417cb1d55e Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 4 Nov 2025 18:43:31 +0100 Subject: [PATCH 34/34] fix ci --- .../src/geos/processing/post_processing/GeosBlockMerge.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py index 53588c6b..b1188ddf 100644 --- a/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py +++ b/geos-processing/src/geos/processing/post_processing/GeosBlockMerge.py @@ -5,8 +5,8 @@ import logging from typing_extensions import Self -from vtkmodules.vtkCommonDataModel import ( vtkCompositeDataSet, vtkMultiBlockDataSet, vtkPolyData, vtkUnstructuredGrid, - vtkCellTypes, vtkTriangle ) +from vtkmodules.vtkCommonDataModel import ( vtkCompositeDataSet, vtkMultiBlockDataSet, vtkPolyData, + vtkUnstructuredGrid ) from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter from geos.utils.Logger import ( Logger, getLogger )