diff --git a/docs/geos_mesh_docs/processing.rst b/docs/geos_mesh_docs/processing.rst index 003cc973d..2fccb6819 100644 --- a/docs/geos_mesh_docs/processing.rst +++ b/docs/geos_mesh_docs/processing.rst @@ -3,6 +3,14 @@ Processing filters The `processing` module of `geos-mesh` package contains filters to process meshes. +geos.mesh.processing.AttributeMapping filter +--------------------------------------------- + +.. automodule:: geos.mesh.processing.AttributeMapping + :members: + :undoc-members: + :show-inheritance: + geos.mesh.processing.CreateConstantAttributePerRegion filter ------------------------------------------------------------- diff --git a/docs/geos_posp_docs/PVplugins.rst b/docs/geos_posp_docs/PVplugins.rst index 293e0a65d..dbb329ee5 100644 --- a/docs/geos_posp_docs/PVplugins.rst +++ b/docs/geos_posp_docs/PVplugins.rst @@ -10,11 +10,6 @@ The plugins include: * Processing plugins to compute additional geomechanical properties; * Visualization plugins to plot Mohr's circles and cross plots using Paraview Python View. -PVAttributeMapping plugin ------------------------------------ - -.. automodule:: PVplugins.PVAttributeMapping - PVExtractMergeBlocksVolume plugin ------------------------------------------- @@ -85,10 +80,3 @@ PVplugins.PVSurfaceGeomechanics module .. automodule:: PVplugins.PVSurfaceGeomechanics - -PVTransferAttributesVolumeSurface plugin --------------------------------------------------- - -.. automodule:: PVplugins.PVTransferAttributesVolumeSurface - - diff --git a/docs/geos_posp_docs/filters.rst b/docs/geos_posp_docs/filters.rst index da4b0e559..30ec1c8d8 100644 --- a/docs/geos_posp_docs/filters.rst +++ b/docs/geos_posp_docs/filters.rst @@ -3,22 +3,6 @@ vtk Filters This package defines vtk filters that allows to process Geos outputs. -geos_posp.filters.AttributeMappingFromCellCoords module --------------------------------------------------------- - -.. automodule:: geos_posp.filters.AttributeMappingFromCellCoords - :members: - :undoc-members: - :show-inheritance: - -geos_posp.filters.AttributeMappingFromCellId module ------------------------------------------------------------ - -.. automodule:: geos_posp.filters.AttributeMappingFromCellId - :members: - :undoc-members: - :show-inheritance: - geos_posp.filters.GeomechanicsCalculator module --------------------------------------------------- @@ -50,19 +34,3 @@ geos_posp.filters.SurfaceGeomechanics module :members: :undoc-members: :show-inheritance: - -geos_posp.filters.TransferAttributesVolumeSurface module ------------------------------------------------------------- - -.. automodule:: geos_posp.filters.TransferAttributesVolumeSurface - :members: - :undoc-members: - :show-inheritance: - -geos_posp.filters.VolumeSurfaceMeshMapper module ----------------------------------------------------- - -.. automodule:: geos_posp.filters.VolumeSurfaceMeshMapper - :members: - :undoc-members: - :show-inheritance: \ No newline at end of file diff --git a/docs/geos_pv_docs/processing.rst b/docs/geos_pv_docs/processing.rst index 7655ede61..dde23c0f5 100644 --- a/docs/geos_pv_docs/processing.rst +++ b/docs/geos_pv_docs/processing.rst @@ -1,6 +1,10 @@ Post-/Pre-processing ========================= +PVAttributeMapping +-------------------- + +.. automodule:: geos.pv.plugins.PVAttributeMapping PVCreateConstantAttributePerRegion ----------------------------------- diff --git a/geos-mesh/src/geos/mesh/processing/AttributeMapping.py b/geos-mesh/src/geos/mesh/processing/AttributeMapping.py new file mode 100644 index 000000000..465f986b5 --- /dev/null +++ b/geos-mesh/src/geos/mesh/processing/AttributeMapping.py @@ -0,0 +1,199 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Raphaël Vinour, Martin Lemay, Romain Baville +# ruff: noqa: E402 # disable Module level import not at top of file +import numpy as np +import numpy.typing as npt +import logging + +from geos.utils.Logger import ( Logger, getLogger ) +from typing_extensions import Self, Union, Set, List, Dict + +from vtkmodules.vtkCommonDataModel import ( + vtkDataSet, + vtkMultiBlockDataSet, +) + +from geos.mesh.utils.arrayModifiers import transferAttributeWithElementMap +from geos.mesh.utils.arrayHelpers import ( computeElementMapping, getAttributeSet, isAttributeGlobal ) + +__doc__ = """ +AttributeMapping is a vtk filter that transfers global attributes from a source mesh to a final mesh with same point/cell coordinates. The final mesh is updated directly, without creation of a copy. + +Input meshes can be vtkDataSet or vtkMultiBlockDataSet. + +.. Warning:: + For one application of the filter, the attributes to transfer should all be located on the same piece (all on points or all on cells). + +.. Note:: + For cell, the coordinates of the points in the cell are compared. + If one of the two meshes is a surface and the other a volume, all the points of the surface must be points of the volume. + +To use a logger handler of yours, set the variable 'speHandler' to True and add it using the member function setLoggerHandler. + +To use the filter: + +.. code-block:: python + + from geos.mesh.processing.AttributeMapping import AttributeMapping + + # Filter inputs. + meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet ] + meshTo: Union[ vtkDataSet, vtkMultiBlockDataSet ] + attributeNames: Set[ str ] + # Optional inputs. + onPoints: bool # defaults to False + speHandler: bool # defaults to False + + # Instantiate the filter + filter: AttributeMapping = AttributeMapping( meshFrom, + meshTo, + attributeNames, + onPoints, + speHandler, + ) + + # Set the handler of yours (only if speHandler is True). + yourHandler: logging.Handler + filter.setLoggerHandler( yourHandler ) + + # Do calculations. + filter.applyFilter() +""" + +loggerTitle: str = "Attribute Mapping" + + +class AttributeMapping: + + def __init__( + self: Self, + meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet ], + meshTo: Union[ vtkDataSet, vtkMultiBlockDataSet ], + attributeNames: Set[ str ], + onPoints: bool = False, + speHandler: bool = False, + ) -> None: + """Transfer global attributes from a source mesh to a final mesh, mapping the piece of the attributes to transfer. + + Args: + meshFrom (Union[ vtkDataSet, vtkMultiBlockDataSet ]): The source mesh with attributes to transfer. + meshTo (Union[ vtkDataSet, vtkMultiBlockDataSet ]): The final mesh where to transfer attributes. + attributeNames (Set[str]): Names of the attributes to transfer. + onPoints (bool): True if attributes are on points, False if they are on cells. + Defaults to False. + speHandler (bool, optional): True to use a specific handler, False to use the internal handler. + Defaults to False. + """ + self.meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet ] = meshFrom + self.meshTo: Union[ vtkDataSet, vtkMultiBlockDataSet ] = meshTo + self.attributeNames: Set[ str ] = attributeNames + self.onPoints: bool = onPoints + #TODO/refact (@RomainBaville) make it an enum + self.piece: str = "points" if self.onPoints else "cells" + + # cell map + self.ElementMap: Dict[ int, npt.NDArray[ np.int64 ] ] = {} + + # 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. + + 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 getElementMap( self: Self ) -> Dict[ int, npt.NDArray[ np.int64 ] ]: + """Getter of the element mapping dictionary. + + If attribute to transfer are on points it will be a pointMap, else it will be a cellMap. + + Returns: + self.elementMap (Dict[int, npt.NDArray[np.int64]]): The element mapping dictionary. + """ + return self.ElementMap + + def applyFilter( self: Self ) -> bool: + """Transfer global attributes from a source mesh to a final mesh, mapping the piece of the attributes to transfer. + + Returns: + boolean (bool): True if calculation successfully ended, False otherwise. + """ + self.logger.info( f"Apply filter { self.logger.name }." ) + + if len( self.attributeNames ) == 0: + self.logger.warning( f"Please enter at least one { self.piece } attribute to transfer." ) + self.logger.warning( f"The filter { self.logger.name } has not been used." ) + return False + + attributesInMeshFrom: Set[ str ] = getAttributeSet( self.meshFrom, self.onPoints ) + wrongAttributeNames: Set[ str ] = self.attributeNames.difference( attributesInMeshFrom ) + if len( wrongAttributeNames ) > 0: + self.logger.error( + f"The { self.piece } attributes { wrongAttributeNames } are not present in the source mesh." ) + self.logger.error( f"The filter { self.logger.name } failed." ) + return False + + attributesInMeshTo: Set[ str ] = getAttributeSet( self.meshTo, self.onPoints ) + attributesAlreadyInMeshTo: Set[ str ] = self.attributeNames.intersection( attributesInMeshTo ) + if len( attributesAlreadyInMeshTo ) > 0: + self.logger.error( + f"The { self.piece } attributes { attributesAlreadyInMeshTo } are already present in the final mesh." ) + self.logger.error( f"The filter { self.logger.name } failed." ) + return False + + if isinstance( self.meshFrom, vtkMultiBlockDataSet ): + partialAttributes: List[ str ] = [] + for attributeName in self.attributeNames: + if not isAttributeGlobal( self.meshFrom, attributeName, self.onPoints ): + partialAttributes.append( attributeName ) + + if len( partialAttributes ) > 0: + self.logger.error( + f"All { self.piece } attributes to transfer must be global, { partialAttributes } are partials." ) + self.logger.error( f"The filter { self.logger.name } failed." ) + + self.ElementMap = computeElementMapping( self.meshFrom, self.meshTo, self.onPoints ) + sharedElement: bool = False + for key in self.ElementMap: + if np.any( self.ElementMap[ key ] > -1 ): + sharedElement = True + + if not sharedElement: + self.logger.warning( f"The two meshes do not have any shared { self.piece }." ) + self.logger.warning( f"The filter { self.logger.name } has not been used." ) + return False + + for attributeName in self.attributeNames: + if not transferAttributeWithElementMap( self.meshFrom, self.meshTo, self.ElementMap, attributeName, + self.onPoints, self.logger ): + self.logger.error( f"The attribute { attributeName } has not been mapped." ) + self.logger.error( f"The filter { self.logger.name } failed." ) + return False + + # Log the output message. + self._logOutputMessage() + + return True + + def _logOutputMessage( self: Self ) -> None: + """Create and log result messages of the filter.""" + self.logger.info( f"The filter { self.logger.name } succeeded." ) + self.logger.info( + f"The { self.piece } attributes { self.attributeNames } have been transferred from the source mesh to the final mesh with the { self.piece } mapping." + ) diff --git a/geos-mesh/src/geos/mesh/utils/arrayHelpers.py b/geos-mesh/src/geos/mesh/utils/arrayHelpers.py index 802a552bf..31f6ecdb1 100644 --- a/geos-mesh/src/geos/mesh/utils/arrayHelpers.py +++ b/geos-mesh/src/geos/mesh/utils/arrayHelpers.py @@ -19,12 +19,256 @@ ArrayHelpers module contains several utilities methods to get information on arrays in VTK datasets. These methods include: + - mesh element localization mapping by indexes - array getters, with conversion into numpy array or pandas dataframe - boolean functions to check whether an array is present in the dataset - bounds getter for vtu and multiblock datasets """ +def computeElementMapping( + meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet ], + meshTo: Union[ vtkDataSet, vtkMultiBlockDataSet ], + points: bool, +) -> dict[ int, npt.NDArray ]: + """Compute the map of points/cells between the source mesh and the final mesh. + + If the source mesh is a vtkDataSet, its flat index (flatIdDataSetFrom) is set to 0. + If the final mesh is a vtkDataSet, its flat index (flatIdDataSetTo) is set to 0. + + The elementMap is a dictionary where: + - Keys are the flat index of all the datasets of the final mesh. + - Items are arrays of size (nb elements in datasets of the final mesh, 2). + + For each element (idElementTo) of each dataset (flatIdDataSetTo) of the final mesh, + if the coordinates of an element (idElementFrom) of one dataset (flatIdDataSetFrom) of the source mesh + are the same as the coordinates of the element of the final mesh, + elementMap[flatIdDataSetTo][idElementTo] = [flatIdDataSetFrom, idElementFrom] + else, elementMap[flatIdDataSetTo][idElementTo] = [-1, -1]. + + For cells, the coordinates of the points in the cell are compared. + If one of the two meshes is a surface and the other a volume, all the points of the surface must be points of the volume. + + Args: + meshFrom (Union[vtkDataSet, vtkMultiBlockDataSet]): The source mesh with the element to map. + meshTo (Union[vtkDataSet, vtkMultiBlockDataSet]): The final mesh with the reference element coordinates. + points (bool): True if elements to map are points, False if they are cells. + + Returns: + elementMap (dict[int, npt.NDArray[np.int64]]): The map of points/cells between the two meshes. + """ + elementMap: dict[ int, npt.NDArray ] = {} + if isinstance( meshTo, vtkDataSet ): + UpdateElementMappingToDataSet( meshFrom, meshTo, elementMap, points ) + elif isinstance( meshTo, vtkMultiBlockDataSet ): + UpdateElementMappingToMultiBlockDataSet( meshFrom, meshTo, elementMap, points ) + + return elementMap + + +def UpdateElementMappingToMultiBlockDataSet( + meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet ], + multiBlockDataSetTo: vtkMultiBlockDataSet, + elementMap: dict[ int, npt.NDArray ], + points: bool, +) -> None: + """Update the map of points/cells between the source mesh and the final mesh. + + If the source mesh is a vtkDataSet, its flat index (flatIdDataSetFrom) is set to 0. + + Add the mapping for of the final mesh: + - Keys are the flat index of all the datasets of the final mesh. + - Items are arrays of size (nb elements in datasets of the final mesh, 2). + + For each element (idElementTo) of each dataset (flatIdDataSetTo) of final mesh, + if the coordinates of an element (idElementFrom) of one dataset (flatIdDataSetFrom) of the source mesh + are the same as the coordinates of the element of the final mesh, + elementMap[flatIdDataSetTo][idElementTo] = [flatIdDataSetFrom, idElementFrom] + else, elementMap[flatIdDataSetTo][idElementTo] = [-1, -1]. + + For cells, the coordinates of the points in the cell are compared. + If one of the two meshes is a surface and the other a volume, all the points of the surface must be points of the volume. + + Args: + meshFrom (Union[vtkDataSet, vtkMultiBlockDataSet]): The source mesh with the element to map. + multiBlockDataSetTo (vtkMultiBlockDataSet): The final mesh with the reference element coordinates. + elementMap (dict[int, npt.NDArray[np.int64]]): The map of points/cells to update. + points (bool): True if elements to map are points, False if they are cells. + """ + listFlatIdDataSetTo: list[ int ] = getBlockElementIndexesFlatten( multiBlockDataSetTo ) + for flatIdDataSetTo in listFlatIdDataSetTo: + dataSetTo: vtkDataSet = vtkDataSet.SafeDownCast( multiBlockDataSetTo.GetDataSet( flatIdDataSetTo ) ) + UpdateElementMappingToDataSet( meshFrom, dataSetTo, elementMap, points, flatIdDataSetTo=flatIdDataSetTo ) + + +def UpdateElementMappingToDataSet( + meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet ], + dataSetTo: vtkDataSet, + elementMap: dict[ int, npt.NDArray ], + points: bool, + flatIdDataSetTo: int = 0, +) -> None: + """Update the map of points/cells between the source mesh and the final mesh. + + If the source mesh is a vtkDataSet, its flat index (flatIdDataSetFrom) is set to 0. + + Add the mapping for the final mesh: + - The key is the flat index of the final mesh. + - The item is an array of size (nb elements in the final mesh, 2). + + For each element (idElementTo) of the final mesh, + if the coordinates of an element (idElementFrom) of one dataset (flatIdDataSetFrom) of the source mesh + are the same as the coordinates of the element of the final mesh, + elementMap[flatIdDataSetTo][idElementTo] = [flatIdDataSetFrom, idElementFrom] + else, elementMap[flatIdDataSetTo][idElementTo] = [-1, -1]. + + For cells, the coordinates of the points in the cell are compared. + If one of the two meshes is a surface and the other a volume, all the points of the surface must be points of the volume. + + Args: + meshFrom (Union[vtkDataSet, vtkMultiBlockDataSet]): The source mesh with the element to map. + dataSetTo (vtkDataSet): The final mesh with the reference element coordinates. + elementMap (dict[int, npt.NDArray[np.int64]]): The map of points/cells to update. + points (bool): True if elements to map are points, False if they are cells. + flatIdDataSetTo (int, Optional): The flat index of the final mesh considered as a dataset of a vtkMultiblockDataSet. + Defaults to 0 for final meshes who are not datasets of vtkMultiBlockDataSet. + """ + nbElementsTo: int = dataSetTo.GetNumberOfPoints() if points else dataSetTo.GetNumberOfCells() + elementMap[ flatIdDataSetTo ] = np.full( ( nbElementsTo, 2 ), -1, np.int64 ) + if isinstance( meshFrom, vtkDataSet ): + UpdateDictElementMappingFromDataSetToDataSet( meshFrom, + dataSetTo, + elementMap, + points, + flatIdDataSetTo=flatIdDataSetTo ) + elif isinstance( meshFrom, vtkMultiBlockDataSet ): + UpdateElementMappingFromMultiBlockDataSetToDataSet( meshFrom, + dataSetTo, + elementMap, + points, + flatIdDataSetTo=flatIdDataSetTo ) + + +def UpdateElementMappingFromMultiBlockDataSetToDataSet( + multiBlockDataSetFrom: vtkMultiBlockDataSet, + dataSetTo: vtkDataSet, + elementMap: dict[ int, npt.NDArray ], + points: bool, + flatIdDataSetTo: int = 0, +) -> None: + """Update the map of points/cells between the source mesh and the final mesh. + + For each element (idElementTo) of the final mesh not yet mapped (elementMap[flatIdDataSetTo][idElementTo] = [-1, -1]), + if the coordinates of an element (idElementFrom) of one dataset (flatIdDataSetFrom) of the source mesh + are the same as the coordinates of the element of the final mesh, + the map of points/cells is update: elementMap[flatIdDataSetTo][idElementTo] = [flatIdDataSetFrom, idElementFrom]. + + For cells, the coordinates of the points in the cell are compared. + If one of the two meshes is a surface and the other a volume, all the points of the surface must be points of the volume. + + Args: + multiBlockDataSetFrom (vtkMultiBlockDataSet): The source mesh with the element to map. + dataSetTo (vtkDataSet): The final mesh with the reference element coordinates. + elementMap (dict[int, npt.NDArray[np.int64]]): The map of points/cells to update with; + The key is the flat index of the final mesh. + The item is an array of size (nb elements in the final mesh, 2). + points (bool): True if elements to map are points, False if they are cells. + flatIdDataSetTo (int, Optional): The flat index of the final mesh considered as a dataset of a vtkMultiblockDataSet. + Defaults to 0 for final meshes who are not datasets of vtkMultiBlockDataSet. + """ + listFlatIDataSetFrom: list[ int ] = getBlockElementIndexesFlatten( multiBlockDataSetFrom ) + for flatIdDataSetFrom in listFlatIDataSetFrom: + dataSetFrom: vtkDataSet = vtkDataSet.SafeDownCast( multiBlockDataSetFrom.GetDataSet( flatIdDataSetFrom ) ) + UpdateDictElementMappingFromDataSetToDataSet( dataSetFrom, + dataSetTo, + elementMap, + points, + flatIdDataSetFrom=flatIdDataSetFrom, + flatIdDataSetTo=flatIdDataSetTo ) + + +def UpdateDictElementMappingFromDataSetToDataSet( + dataSetFrom: vtkDataSet, + dataSetTo: vtkDataSet, + elementMap: dict[ int, npt.NDArray[ np.int64 ] ], + points: bool, + flatIdDataSetFrom: int = 0, + flatIdDataSetTo: int = 0, +) -> None: + """Update the map of points/cells between the source mesh and the final mesh. + + For each element (idElementTo) of the final mesh not yet mapped (elementMap[flatIdDataSetTo][idElementTo] = [-1, -1]), + if the coordinates of an element (idElementFrom) of the source mesh + are the same as the coordinates of the element of the final mesh, + the map of points/cells is update: elementMap[flatIdDataSetTo][idElementTo] = [flatIdDataSetFrom, idElementFrom]. + + For cells, the coordinates of the points in the cell are compared. + If one of the two meshes is a surface and the other a volume, all the points of the surface must be points of the volume. + + Args: + dataSetFrom (vtkDataSet): The source mesh with the element to map. + dataSetTo (vtkDataSet): The final mesh with the reference element coordinates. + elementMap (dict[int, npt.NDArray[np.int64]]): The map of points/cells to update with; + The key is the flat index of the final mesh. + The item is an array of size (nb elements in the final mesh, 2). + points (bool): True if elements to map are points, False if they are cells. + flatIdDataSetFrom (int, Optional): The flat index of the source mesh considered as a dataset of a vtkMultiblockDataSet. + Defaults to 0 for source meshes who are not datasets of vtkMultiBlockDataSet. + flatIdDataSetTo (int, Optional): The flat index of the final mesh considered as a dataset of a vtkMultiblockDataSet. + Defaults to 0 for final meshes who are not datasets of vtkMultiBlockDataSet. + """ + idElementsFromFund: list[ int ] = [] + nbElementsTo: int = len( elementMap[ flatIdDataSetTo ] ) + nbElementsFrom: int = dataSetFrom.GetNumberOfPoints() if points else dataSetFrom.GetNumberOfCells() + for idElementTo in range( nbElementsTo ): + # Test if the element of the final mesh is already mapped. + if -1 in elementMap[ flatIdDataSetTo ][ idElementTo ]: + coordElementTo: list[ tuple[ float, ...] ] = [] + if points: + coordElementTo.append( dataSetTo.GetPoint( idElementTo ) ) + else: + # Get the coordinates of each points of the cell. + nbPointsTo: int = dataSetTo.GetCell( idElementTo ).GetNumberOfPoints() + cellPointsTo: vtkPoints = dataSetTo.GetCell( idElementTo ).GetPoints() + for idPointTo in range( nbPointsTo ): + coordElementTo.append( cellPointsTo.GetPoint( idPointTo ) ) + + idElementFrom: int = 0 + ElementFromFund: bool = False + while idElementFrom < nbElementsFrom and not ElementFromFund: + # Test if the element of the source mesh is already mapped. + if idElementFrom not in idElementsFromFund: + coordElementFrom: list[ tuple[ float, ...] ] = [] + if points: + coordElementFrom.append( dataSetFrom.GetPoint( idElementFrom ) ) + else: + # Get the coordinates of each points of the cell. + nbPointsFrom: int = dataSetFrom.GetCell( idElementFrom ).GetNumberOfPoints() + cellPointsFrom: vtkPoints = dataSetFrom.GetCell( idElementFrom ).GetPoints() + for idPointFrom in range( nbPointsFrom ): + coordElementFrom.append( cellPointsFrom.GetPoint( idPointFrom ) ) + + pointShared: bool = True + if dataSetTo.GetClassName() == dataSetFrom.GetClassName(): + if coordElementTo != coordElementFrom: + pointShared = False + elif isinstance( dataSetTo, vtkPolyData ): + for coordPointsTo in coordElementTo: + if coordPointsTo not in coordElementFrom: + pointShared = False + elif isinstance( dataSetFrom, vtkPolyData ): + for coordPointsFrom in coordElementFrom: + if coordPointsFrom not in coordElementTo: + pointShared = False + + if pointShared: + elementMap[ flatIdDataSetTo ][ idElementTo ] = [ flatIdDataSetFrom, idElementFrom ] + ElementFromFund = True + idElementsFromFund.append( idElementFrom ) + + idElementFrom += 1 + + def has_array( mesh: vtkUnstructuredGrid, array_names: list[ str ] ) -> bool: """Checks if input mesh contains at least one of input data arrays. @@ -652,24 +896,28 @@ def getComponentNamesMultiBlock( return () -def getAttributeValuesAsDF( surface: vtkPolyData, attributeNames: tuple[ str, ...] ) -> pd.DataFrame: +def getAttributeValuesAsDF( surface: vtkPolyData, + attributeNames: tuple[ str, ...], + onPoints: bool = False ) -> pd.DataFrame: """Get attribute values from input surface. Args: surface (vtkPolyData): Mesh where to get attribute values. attributeNames (tuple[str,...]): Tuple of attribute names to get the values. + onPoints (bool, optional): True if attributes are on points, False if they are on cells. + Defaults to False. Returns: pd.DataFrame: DataFrame containing property names as columns. """ - nbRows: int = surface.GetNumberOfCells() + nbRows: int = surface.GetNumberOfPoints() if onPoints else surface.GetNumberOfCells() data: pd.DataFrame = pd.DataFrame( np.full( ( nbRows, len( attributeNames ) ), np.nan ), columns=attributeNames ) for attributeName in attributeNames: - if not isAttributeInObject( surface, attributeName, False ): + if not isAttributeInObject( surface, attributeName, onPoints ): logging.warning( f"Attribute {attributeName} is not in the mesh." ) continue - array: npt.NDArray[ np.float64 ] = getArrayInObject( surface, attributeName, False ) + array: npt.NDArray[ np.float64 ] = getArrayInObject( surface, attributeName, onPoints ) if len( array.shape ) > 1: for i in range( array.shape[ 1 ] ): diff --git a/geos-mesh/src/geos/mesh/utils/arrayModifiers.py b/geos-mesh/src/geos/mesh/utils/arrayModifiers.py index 2724acda7..165f2de12 100644 --- a/geos-mesh/src/geos/mesh/utils/arrayModifiers.py +++ b/geos-mesh/src/geos/mesh/utils/arrayModifiers.py @@ -40,6 +40,7 @@ isAttributeGlobal, getVtkArrayTypeInObject, getVtkArrayTypeInMultiBlock, + getVtkDataTypeInObject, getNumberOfComponentsMultiBlock, ) from geos.mesh.utils.multiblockHelpers import getBlockElementIndexesFlatten @@ -50,6 +51,8 @@ These methods include: - filling partial VTK arrays with values (useful for block merge) - creation of new VTK array, empty or with a given data array + - copy VTK array from a source mesh to a final mesh + - transfer VTK array from a source mesh to a final mesh with a element map - transfer from VTK point data to VTK cell data """ @@ -680,6 +683,177 @@ def copyAttributeDataSet( return createAttribute( dataSetTo, npArray, attributeNameTo, componentNames, onPoints, vtkArrayType, logger ) +def transferAttributeToDataSetWithElementMap( + meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet ], + dataSetTo: vtkDataSet, + elementMap: dict[ int, npt.NDArray[ np.int64 ] ], + attributeName: str, + onPoints: bool, + flatIdDataSetTo: int = 0, + logger: Union[ Logger, Any ] = None, +) -> bool: + """Transfer attributes from the source mesh to the final mesh using a map of points/cells. + + If the source mesh is a vtkDataSet, its flat index (flatIdDataSetFrom) is set to 0. + + The map of points/cells used to transfer the attribute is a dictionary where: + - The key is the flat index of the final mesh. + - The item is an array of size (nb elements in the final mesh, 2). + + If an element (idElementTo) of the final mesh is mapped with no element of the source mesh: + - elementMap[flatIdDataSetTo][idElementTo] = [-1, -1]. + - The value of the attribute for this element depends of the type of the value of the attribute (0 for unit, -1 for int, nan for float). + + If an element (idElementTo) of the final mesh is mapped with an element (idElementFrom) of one of the dataset (flatIdDataSetFrom) of the source mesh: + - elementMap[flatIdDataSetTo][idElementTo] = [flatIdDataSetFrom, idElementFrom]. + - The value of the attribute for this element is the value of the element (idElementFrom) of the dataset (flatIdDataSetFrom) of the source mesh. + + Args: + meshFrom (Union[vtkDataSet, vtkMultiBlockDataSet]): The source mesh with the attribute to transfer. + dataSetTo (vtkDataSet): The final mesh where to transfer the attribute. + elementMap (dict[int, npt.NDArray[np.int64]]): The map of points/cells. + attributeName (str): The name of the attribute to transfer. + onPoints (bool): True if the attribute is on points, False if it is on cells. + flatIdDataSetTo (int, Optional): The flat index of the final mesh considered as a dataset of a vtkMultiblockDataSet. + Defaults to 0 for final meshes who are not datasets of vtkMultiBlockDataSet. + logger (Union[Logger, None], optional): A logger to manage the output messages. + Defaults to None, an internal logger is used. + + Returns: + bool: True if transfer successfully ended. + """ + # Check if an external logger is given. + if logger is None: + logger = getLogger( "transferAttributeToDataSetWithElementMap", True ) + + if flatIdDataSetTo not in elementMap: + logger.error( f"The map is incomplete, there is no data for the final mesh (flat index { flatIdDataSetTo })." ) + return False + + nbElementsTo: int = dataSetTo.GetNumberOfPoints() if onPoints else dataSetTo.GetNumberOfCells() + if len( elementMap[ flatIdDataSetTo ] ) != nbElementsTo: + logger.error( + f"The map is wrong, there is { nbElementsTo } elements in the final mesh (flat index { flatIdDataSetTo })\ + but { len( elementMap[ flatIdDataSetTo ] ) } elements in the map." ) + return False + + componentNames: tuple[ str, ...] = getComponentNames( meshFrom, attributeName, onPoints ) + nbComponents: int = len( componentNames ) + + vtkDataType: int = getVtkDataTypeInObject( meshFrom, attributeName, onPoints ) + defaultValue: Any + if vtkDataType in ( VTK_FLOAT, VTK_DOUBLE ): + defaultValue = np.nan + elif vtkDataType in ( VTK_CHAR, VTK_SIGNED_CHAR, VTK_SHORT, VTK_LONG, VTK_INT, VTK_LONG_LONG, VTK_ID_TYPE ): + defaultValue = -1 + elif vtkDataType in ( VTK_BIT, VTK_UNSIGNED_CHAR, VTK_UNSIGNED_SHORT, VTK_UNSIGNED_LONG, VTK_UNSIGNED_INT, + VTK_UNSIGNED_LONG_LONG ): + defaultValue = 0 + + typeMapping: dict[ int, type ] = vnp.get_vtk_to_numpy_typemap() + valueType: type = typeMapping[ vtkDataType ] + + arrayTo: npt.NDArray[ Any ] + if nbComponents > 1: + defaultValue = [ defaultValue ] * nbComponents + arrayTo = np.full( ( nbElementsTo, nbComponents ), defaultValue, dtype=valueType ) + else: + arrayTo = np.array( [ defaultValue for _ in range( nbElementsTo ) ], dtype=valueType ) + + for idElementTo in range( nbElementsTo ): + valueToTransfer: Any = defaultValue + idElementFrom: int = elementMap[ flatIdDataSetTo ][ idElementTo ][ 1 ] + if idElementFrom != -1: + dataFrom: Union[ vtkPointData, vtkCellData ] + if isinstance( meshFrom, vtkDataSet ): + dataFrom = meshFrom.GetPointData() if onPoints else meshFrom.GetCellData() + elif isinstance( meshFrom, vtkMultiBlockDataSet ): + flatIdDataSetFrom: int = elementMap[ flatIdDataSetTo ][ idElementTo ][ 0 ] + dataSetFrom: vtkDataSet = vtkDataSet.SafeDownCast( meshFrom.GetDataSet( flatIdDataSetFrom ) ) + dataFrom = dataSetFrom.GetPointData() if onPoints else dataSetFrom.GetCellData() + + arrayFrom: npt.NDArray[ Any ] = vnp.vtk_to_numpy( dataFrom.GetArray( attributeName ) ) + valueToTransfer = arrayFrom[ idElementFrom ] + + arrayTo[ idElementTo ] = valueToTransfer + + return createAttribute( dataSetTo, + arrayTo, + attributeName, + componentNames, + onPoints=onPoints, + vtkDataType=vtkDataType, + logger=logger ) + + +def transferAttributeWithElementMap( + meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet ], + meshTo: Union[ vtkDataSet, vtkMultiBlockDataSet ], + elementMap: dict[ int, npt.NDArray[ np.int64 ] ], + attributeName: str, + onPoints: bool, + logger: Union[ Logger, Any ] = None, +) -> bool: + """Transfer attributes from the source mesh to the final mesh using a map of points/cells. + + If the source mesh is a vtkDataSet, its flat index (flatIdDataSetFrom) is set to 0. + If the final mesh is a vtkDataSet, its flat index (flatIdDataSetTo) is set to 0. + + The map of points/cells used to transfer the attribute is a dictionary where: + - Keys are the flat index of all the datasets of the final mesh. + - Items are arrays of size (nb elements in datasets of the final mesh, 2). + + If an element (idElementTo) of one dataset (flatIdDataSetTo) of the final mesh is mapped with no element of the source mesh: + - elementMap[flatIdDataSetTo][idElementTo] = [-1, -1]. + - The value of the attribute for this element depends of the type of the value of the attribute (0 for unit, -1 for int, nan for float). + + If an element (idElementTo) of one dataset (flatIdDataSetTo) of the final mesh is mapped with an element (idElementFrom) of one of the dataset (flatIdDataSetFrom) of the source mesh: + - elementMap[flatIdDataSetTo][idElementTo] = [flatIdDataSetFrom, idElementFrom]. + - The value of the attribute for this element is the value of the element (idElementFrom) of the dataset (flatIdDataSetFrom) of the source mesh. + + Args: + meshFrom (Union[vtkDataSet, vtkMultiBlockDataSet]): The source mesh with the attribute to transfer. + meshTo (Union[vtkDataSet, vtkMultiBlockDataSet]): The final mesh where to transfer the attribute. + elementMap (dict[int, npt.NDArray[np.int64]]): The map of points/cells. + attributeName (str): The name of the attribute to transfer. + onPoints (bool): True if the attribute is on points, False if it is on cells. + Defaults to 0 for final meshes who are not datasets of vtkMultiBlockDataSet. + logger (Union[Logger, None], optional): A logger to manage the output messages. + Defaults to None, an internal logger is used. + + Returns: + bool: True if transfer successfully ended. + """ + # Check if an external logger is given. + if logger is None: + logger = getLogger( "transferAttributeWithElementMap", True ) + + if isinstance( meshTo, vtkDataSet ): + return transferAttributeToDataSetWithElementMap( meshFrom, + meshTo, + elementMap, + attributeName, + onPoints, + logger=logger ) + elif isinstance( meshTo, vtkMultiBlockDataSet ): + listFlatIdDataSetTo: list[ int ] = getBlockElementIndexesFlatten( meshTo ) + for flatIdDataSetTo in listFlatIdDataSetTo: + dataSetTo: vtkDataSet = vtkDataSet.SafeDownCast( meshTo.GetDataSet( flatIdDataSetTo ) ) + if not transferAttributeToDataSetWithElementMap( meshFrom, + dataSetTo, + elementMap, + attributeName, + onPoints, + flatIdDataSetTo=flatIdDataSetTo, + logger=logger ): + logger.error( + f"The attribute transfer has failed for the dataset with the flat index { flatIdDataSetTo } of the final mesh." + ) + logger.warning( "The final mesh may has been modify for the other datasets." ) + return False + return True + + def renameAttribute( object: Union[ vtkMultiBlockDataSet, vtkDataSet ], attributeName: str, diff --git a/geos-mesh/src/geos/mesh/utils/multiblockModifiers.py b/geos-mesh/src/geos/mesh/utils/multiblockModifiers.py index 5f00afb82..f579af291 100644 --- a/geos-mesh/src/geos/mesh/utils/multiblockModifiers.py +++ b/geos-mesh/src/geos/mesh/utils/multiblockModifiers.py @@ -3,7 +3,7 @@ # SPDX-FileContributor: Martin Lemay from typing import Union from vtkmodules.vtkCommonDataModel import ( vtkCompositeDataSet, vtkDataObjectTreeIterator, vtkMultiBlockDataSet, - vtkUnstructuredGrid ) + vtkDataSet ) from vtkmodules.vtkFiltersCore import vtkAppendDataSets from geos.mesh.utils.arrayModifiers import fillAllPartialAttributes @@ -14,7 +14,7 @@ def mergeBlocks( input: Union[ vtkMultiBlockDataSet, vtkCompositeDataSet ], keepPartialAttributes: bool = False, -) -> vtkUnstructuredGrid: +) -> vtkDataSet: """Merge all blocks of a multi block mesh. Args: @@ -38,7 +38,7 @@ def mergeBlocks( iter.VisitOnlyLeavesOn() iter.GoToFirstItem() while iter.GetCurrentDataObject() is not None: - block: vtkUnstructuredGrid = vtkUnstructuredGrid.SafeDownCast( iter.GetCurrentDataObject() ) + block: vtkDataSet = vtkDataSet.SafeDownCast( iter.GetCurrentDataObject() ) af.AddInputData( block ) iter.GoToNextItem() af.Update() diff --git a/geos-mesh/tests/conftest.py b/geos-mesh/tests/conftest.py index ea0810593..5d4b74d51 100644 --- a/geos-mesh/tests/conftest.py +++ b/geos-mesh/tests/conftest.py @@ -5,12 +5,12 @@ # ruff: noqa: E402 # disable Module level import not at top of file import os import pytest -from typing import Union, Any +from typing import Union, Any, Tuple, Dict import numpy as np import numpy.typing as npt from vtkmodules.vtkCommonDataModel import vtkDataSet, vtkMultiBlockDataSet, vtkPolyData -from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridReader, vtkXMLMultiBlockDataReader +from vtkmodules.vtkIOXML import vtkXMLGenericDataObjectReader @pytest.fixture @@ -158,22 +158,23 @@ def _get_dataset( datasetType: str ) -> Union[ vtkMultiBlockDataSet, vtkPolyData Returns: (vtkMultiBlockDataSet, vtkPolyData, vtkDataSet): The vtk object. """ - reader: Union[ vtkXMLMultiBlockDataReader, vtkXMLUnstructuredGridReader ] + reader: vtkXMLGenericDataObjectReader = vtkXMLGenericDataObjectReader() if datasetType == "multiblock": - reader: vtkXMLMultiBlockDataReader = vtkXMLMultiBlockDataReader() vtkFilename = "data/displacedFault.vtm" elif datasetType == "emptymultiblock": - reader: vtkXMLMultiBlockDataReader = vtkXMLMultiBlockDataReader() vtkFilename = "data/displacedFaultempty.vtm" + elif datasetType == "fracture": + vtkFilename = "data/fracture_res5_id.vtu" + elif datasetType == "emptyFracture": + vtkFilename = "data/fracture_res5_id_empty.vtu" elif datasetType == "dataset": - reader: vtkXMLUnstructuredGridReader = vtkXMLUnstructuredGridReader() vtkFilename = "data/domain_res5_id.vtu" elif datasetType == "emptydataset": - reader: vtkXMLUnstructuredGridReader = vtkXMLUnstructuredGridReader() vtkFilename = "data/domain_res5_id_empty.vtu" elif datasetType == "polydata": - reader: vtkXMLUnstructuredGridReader = vtkXMLUnstructuredGridReader() - vtkFilename = "data/triangulatedSurface.vtu" + vtkFilename = "data/fracture_res5_id.vtp" + elif datasetType == "emptypolydata": + vtkFilename = "data/fracture_res5_id_empty.vtp" datapath: str = os.path.join( os.path.dirname( os.path.realpath( __file__ ) ), vtkFilename ) reader.SetFileName( datapath ) reader.Update() @@ -181,3 +182,84 @@ def _get_dataset( datasetType: str ) -> Union[ vtkMultiBlockDataSet, vtkPolyData return reader.GetOutput() return _get_dataset + + +@pytest.fixture +def getElementMap() -> Any: + """Get the element indexes mapping dictionary using the function _get_elementMap() between two meshes. + + Returns: + elementMap (Dict[int, npt.NDArray[np.int64]]): The cell mapping dictionary. + """ + + def _get_elementMap( meshFromName: str, meshToName: str, points: bool ) -> Dict[ int, npt.NDArray[ np.int64 ] ]: + """Get the element indexes mapping dictionary between two meshes. + + Args: + meshFromName (str): The name of the meshFrom. + meshToName (str): The name of the meshTo. + points (bool): True if elements to map is points, False if it is cells. + + Returns: + elementMap (Dict[int, npt.NDArray[np.int64]]): The element mapping dictionary. + """ + elementMap: Dict[ int, npt.NDArray[ np.int64 ] ] = {} + nbElements: Tuple[ int, int ] = ( 4092, 212 ) if points else ( 1740, 156 ) + if meshFromName == "multiblock": + if meshToName == "emptymultiblock": + elementMap[ 1 ] = np.array( [ [ 1, element ] for element in range( nbElements[ 0 ] ) ] ) + elementMap[ 3 ] = np.array( [ [ 3, element ] for element in range( nbElements[ 1 ] ) ] ) + elif meshToName == "emptyFracture": + elementMap[ 0 ] = np.array( [ [ 3, element ] for element in range( nbElements[ 1 ] ) ] ) + elif meshFromName == "dataset": + if meshToName == "emptydataset": + elementMap[ 0 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 0 ] ) ] ) + elif meshToName == "emptyFracture": + elementMap[ 0 ] = np.full( ( nbElements[ 1 ], 2 ), -1, np.int64 ) + elif meshToName == "emptypolydata": + elementMap[ 0 ] = np.array( [ [ 0, 0 ], [ 0, 1 ], [ 0, 2 ], [ 0, 3 ], [ 0, 4 ], [ 0, 5 ], [ 0, 6 ], + [ 0, 7 ], [ 0, 8 ], [ 0, 9 ], [ 0, 10 ], [ 0, 11 ], [ 0, 12 ], [ 0, 13 ], + [ 0, 14 ], [ 0, 15 ], [ 0, 16 ], [ 0, 17 ], [ 0, 18 ], [ 0, + 19 ], [ 0, 20 ], + [ 0, 21 ], [ 0, 22 ], [ 0, 23 ], [ 0, 48 ], [ 0, 50 ], [ 0, + 51 ], [ 0, 54 ], + [ 0, 56 ], [ 0, 57 ], [ 0, 58 ], [ 0, 59 ], [ 0, 60 ], [ 0, + 61 ], [ 0, 62 ], + [ 0, 63 ], [ 0, 64 ], [ 0, 65 ], [ 0, 66 ], [ 0, 67 ], [ 0, + 68 ], [ 0, 69 ], + [ 0, 70 ], [ 0, 71 ], [ 0, 72 ], [ 0, 73 ], [ 0, 74 ], + [ 0, 75 ], [ 0, 76 ], [ 0, 77 ], [ 0, 78 ], [ 0, 79 ], [ 0, 580 ], + [ 0, 581 ], [ 0, 582 ], [ 0, 583 ], [ 0, 584 ], [ 0, 585 ], [ 0, 586 ], + [ 0, 587 ], [ 0, 588 ], [ 0, 589 ], [ 0, 590 ], [ 0, 591 ], [ 0, 592 ], + [ 0, 593 ], [ 0, 594 ], [ 0, 595 ], [ 0, 596 ], [ 0, 597 ], [ 0, 598 ], + [ 0, 599 ], [ 0, 600 ], [ 0, 601 ], [ 0, 602 ], [ 0, 603 ], [ 0, 628 ], + [ 0, 630 ], [ 0, 631 ], [ 0, 634 ], [ 0, 636 ], [ 0, 637 ], [ 0, 638 ], + [ 0, 639 ], [ 0, 640 ], [ 0, 641 ], [ 0, 642 ], [ 0, 643 ], [ 0, 644 ], + [ 0, 645 ], [ 0, 646 ], [ 0, 647 ], [ 0, 648 ], [ 0, 649 ], [ 0, 650 ], + [ 0, 651 ], [ 0, 652 ], [ 0, 653 ], [ 0, 654 ], [ 0, 655 ], [ 0, 656 ], + [ 0, 657 ], [ 0, 658 ], [ 0, 659 ], [ 0, 1160 ], [ 0, 1161 ], [ 0, 1162 ], + [ 0, 1163 ], [ 0, 1164 ], [ 0, 1165 ], [ 0, 1166 ], [ 0, 1167 ], + [ 0, 1168 ], [ 0, 1169 ], [ 0, 1170 ], [ 0, 1171 ], [ 0, 1172 ], + [ 0, 1173 ], [ 0, 1174 ], [ 0, 1175 ], [ 0, 1176 ], [ 0, 1177 ], + [ 0, 1178 ], [ 0, 1179 ], [ 0, 1180 ], [ 0, 1181 ], [ 0, 1182 ], + [ 0, 1183 ], [ 0, 1208 ], [ 0, 1210 ], [ 0, 1211 ], [ 0, 1214 ], + [ 0, 1216 ], [ 0, 1217 ], [ 0, 1218 ], [ 0, 1219 ], [ 0, 1220 ], + [ 0, 1221 ], [ 0, 1222 ], [ 0, 1223 ], [ 0, 1224 ], [ 0, 1225 ], + [ 0, 1226 ], [ 0, 1227 ], [ 0, 1228 ], [ 0, 1229 ], [ 0, 1230 ], + [ 0, 1231 ], [ 0, 1232 ], [ 0, 1233 ], [ 0, 1234 ], [ 0, 1235 ], + [ 0, 1236 ], [ 0, 1237 ], [ 0, 1238 ], [ 0, 1239 ] ] ) + elif meshToName == "emptymultiblock": + elementMap[ 1 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 0 ] ) ] ) + elementMap[ 3 ] = np.full( ( nbElements[ 1 ], 2 ), -1, np.int64 ) + elif meshFromName == "fracture": + if meshToName == "emptyFracture": + elementMap[ 0 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 1 ] ) ] ) + elif meshToName == "emptymultiblock": + elementMap[ 1 ] = np.full( ( nbElements[ 0 ], 2 ), -1, np.int64 ) + elementMap[ 3 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 1 ] ) ] ) + elif meshFromName == "polydata" and meshToName == "emptypolydata": + elementMap[ 0 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 1 ] ) ] ) + + return elementMap + + return _get_elementMap diff --git a/geos-mesh/tests/data/displacedFaultSurface.vtm b/geos-mesh/tests/data/displacedFaultSurface.vtm new file mode 100644 index 000000000..8998b1d96 --- /dev/null +++ b/geos-mesh/tests/data/displacedFaultSurface.vtm @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/geos-mesh/tests/data/displacedFaultemptySurface.vtm b/geos-mesh/tests/data/displacedFaultemptySurface.vtm new file mode 100644 index 000000000..5946185af --- /dev/null +++ b/geos-mesh/tests/data/displacedFaultemptySurface.vtm @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/geos-mesh/tests/data/fracture_res5_id.vtp b/geos-mesh/tests/data/fracture_res5_id.vtp new file mode 100644 index 000000000..0c6b93c19 --- /dev/null +++ b/geos-mesh/tests/data/fracture_res5_id.vtp @@ -0,0 +1,64 @@ + + + + + + + + + 3905.8931117 + + + 5326.4624283 + + + + + + + + + + + + + + + 0 + + + 2304.8861143 + + + + + 0 + + + 2304.8861143 + + + + + + + + + + + + + + + + + + + + + + + + _AQAAAAAAAAAAgAAAAAAAAEANAAAAAAAABgMAAAAAAAA=eF4tyMVSFgAAhVFduHDh+Ah2d2GD2I0NdnchJnZiYmCLgYmtKHZid3djgo0dG2f8z918c27B7Jn+LyVzoIX4BBfmk1yET3FRPs3F+AwX57N8Tkv4S+p5fym+wKX5IpfhS1yWL3M5vsJBfJWvaXl/Bb3ur8g3uBLf5Mp8i6vwba7KdziY7/I9DfFX0/v+UH7A1fkh1+BHXJMfcy1+wrX5KT/TOv66muqvx8+5Pr/gBvySG/IrbsSvuTG/4TQN8zfRdH9TfsvN+B035/fcgj9wS/7IrfgTf9Zwf4Rm+FvzF27DX7ktf+N2/J1nZgm0vX8Wd+BY7siddLa/M8/hudrFP4+7chx34/ncnRdwD17IPbmXLvL35sW8RPv4l3JfXsb9OJ7783IewCt4IEfqSv8gXsUJGuVfzYN5DQ/htTyU1/EwXs/DeYRu8EdzIm/Ukf5NPIo382jewmN4K4/lbTyOx+t2/wTewTt1oj+JJ/Eunsy7eQoncwzv4ak8Tff6p/M+3q8z/Ad4Jh/kWXyIY/kwz+YjPIfn6lH/PD7Gcdwya6DzuRUv4HCO0IX+1ryI2/BiXqJt/Uu5HS/j9hzPHXg5d+ROusLfmVdyF17FCdrVv5q78Rruzmu5B6/jntxL1/t78wbuw4m8Ufv6N3E/3sz9eQsP4K08kCN1m38Qb+co3sE7dbA/iYfwLh7Ku3kYJ/NwHqF7/NG8l0fyPt6vo/wHeDQf5DF8iMfyYR7H4/WIfwIf5Yl8jI/rJH8KT+YTPIVPcgyf4qk8TU/7p/MZPqs5sgWaU8/5c/F5vqC5/Xn0ov+S5vVf5nx8hfPzVS7ABfWavxBf5xta2F9Eb/pvaVH/bS7Gd7g43+USXFLv+UvxfX6gpf1l9KH/kZb1P+Zy/ISD+CmX5wr6zF+RU7kSP+fK/IJfahX/K67KrzmY33AIV9M0fyinc3V+yzX4Hb/Xmv4PXIs/cm3+xHW4rn721+MMrs9fuAF/5W/a0P+dG/EPbsw/OYyb6C9/U/7NzfgPN+e//A+qS/z/AQAAAAAAAAAAgAAAAAAAAKAGAAAAAAAAbgEAAAAAAAA=eF4txdciEAAAAEBRUmlpK9q0aEpbey/tTUMb0d57T6VBO+2NSGjvoflDHrp7uYCA/6o40EGu6moOdnWHuIZrupZDXdt1XNf1XN9hbuCGbuTGbuKmbuZwN3cLRzjSLd3Krd3Gbd3O7R3laHdwR3dyZ3dxjGPd1d3c3T3c070c596Odx/3dT/39wAP9CAneLCHeKiHebhHeKRHebTHeKzHebwneKInebITPcVTPc3TPcMzPcuzPcdzPc/zvcBJTvZCL/JiL3GKl3qZl3uFV3qVVzvVaU73Gmc402u9zuu9wRu9yZu9xVu9zdu9wzu9y7u9x3u9z/t9wAd9yId9xEd9zMd9wid9ylk+7TPO9lmf83lfcI5zfdGXfNlXfNXXfN03nOebvuXbvuO7vuf7fuCHfuTHfuKnzneBC/3MRS72c5f4hUtd5nK/9Cu/9hu/9Tu/9wd/9Cd/9hd/9Td/9w9X+Kd/+bf/+K//uRLqf1dfAQAAAAAAAAAAgAAAAAAAAOAEAAAAAAAAEgEAAAAAAAA=eF4txddCCAAAAMAiozSkoaGh0NAeqGhrSNEg7SFkJKFhlIhoaCBUP9tDdy8XEHAo0Ed81EE+5uM+4ZMOdohPOdRhDneETzvSZxzlaMc41mcd53gnONHnnORkpzjV553mdF/wRV9yhjOd5Wxfdo5zned8F7jQRS52iUt9xVd9zWUud4Wv+4YrXeVq17jWda73TTe40U1u9i23+LZb3eY7vut2d7jTXb7n++72A/e4133u94AHPeRhj3jUDz3mR37sJx73Uz/zc7/whF960q885dd+47ee9oxnPed3fu8P/uh5L/iTF/3ZX7zkr/7mZX/3D6941Wte909veNNb3vYv//Yf7/iv//m/d73nfR8ARZMvOw==AQAAAAAAAAAAgAAAAAAAAPAJAAAAAAAAagIAAAAAAAA=eF5t1aFqK1EUheGJqokoEdVVl3mKoSS3b5AHqC4E+gjja6qihxEjbkxUdZM2LpS62FA/UFlG3Z1N9v5XoVEfw5xz1hwW2UVhv2a8K+JXXuD1CD98v6XrHj994uaA13v8ssEfz/i4wl8NLpb48hFf1/juXvLMJUMl55Zy1kT2H16TN++4bsULXFyJl9vi199yiq9myXohbvHNOy4GfJz8Tb+UuKlwPcd39/i6xpePuFjirwYfV/jjWc7d4PVeMhzw06fk6fHDt6wd3abLC9yM1bvVvz+z8zvYupe27qWte2nrXtq6x9o9tu6lrXtp617aupe27qWte2nrXtq6R565ZKjk3FLOmsj+w2vausc+rXiBrXuSbYuLQp5P5f2Z7CNusXVP8mDrnuTH1j35Xmzdk/vB1j3Jhq17cv/Yuse5G2zdI8MBW/fI02PrHmtHt2nrHvuM1daxLt7B3r0u9nxLe/e6yIC9e7F2j717XXwj9u51cSfYu9fFHWLvXhd3LnnmkqGSc0s5ayL7n7p3tncv9mnFC+zdy2xbfOpePp/K+zPZR9xi717mwd69zI+9e/m92LuX94O9e5kNe/fy/rF3L87dYO9eZDhg717k6bF3L9aeune2dy/2Gat3+b/vM3eItdhn7hBnYZ+5Q2TDPnOH+BbsM3eIb8c+c2P2LLHP3CHuFvvMjTxzyVDJuaWcNZH9TzP3bJ+5sU8rXmCfuZlt+2NG8nwq789kH3GLfeZmHuwzN/Njn7n5vdhnbt4P9pmb2bDP3DhrhX3mxrkb7DM3Mhywz9zI02OfubH2NHPP9pkb+4zd/wFdP06cAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAQAAAAAAAAAAgAAAAAAAAIATAAAAAAAAZwMAAAAAAAA=eF5d1cMWGGYABeHUSm3bbmqmtu3Utp3atpvadpsitfVCXXS+Rf+7mXPmzP5OmvDfNo6nxVPjRnHS4E+JGw4df3LcYOj4k+L6Q8efGNcbOv6EuO7Q8cfHdYaOPy6uPXT8lLjW0PHHxjWHjj8mrjF0/NFx9aHjj4qrDR1/ZFx16Pgj4ipDxx8eVx46/rC40tDxh8YVh44/JK4wdPzBcfmh4w+Kyw0df2Bcduj4A+IyQ8fvH5ceOn6/uNTQ8fvGJYeO3ycuMXT83nHxoeP3iosNHb9nXHTo+D3iIkPH7x4XHjp+t7jQ0PG7xgWHjt8lLjB0/M5x/qHjd4rzDR2/Y5x36Pgd4sSh47eP8wwdPznOPXT8dnGuoeO3jXMOHb9NnGPo+K3j7EPHbxVnGzp+yzjr0PFbxFmGjt88zjx0/GZxpqHjN53w/+n4TaI/8TtPxSejP9HxT0R/ouMfj/5Exz8W/YmOfzT6Ex3/SPQnOv7h6E90/ENxytDxD0Z/ouMfiP5Ex98f/YmOvy/6Ex1/b/QnOv6e6E90/N3Rn+j4u6I/0fF3Rn+i4++I/kTH3x79iY6/LfoTHX9r9Cc6/pboT3T8zdGf6Piboj/R8TdGf6Ljb4j+RMdfH/2Jjr8u+hMdf230Jzp+avQnOv6a6E90/NXRn+j4q6I/0fFXRn+i46+I/kTHXx79iY6/LPoTHX9p9Cc6/pLoT3T8xXHy0PEXRX+i4y+M/kTHXxD9iY4/P/oTHX9e9Cc6/tzoT3T8OdGf6Pizoz/R8WdFf6Ljz4z+RMefEf2Mjj89+hO/80/8O/oTHf9X9Cc6/s/oT3T8H9Gf6Pjfoz/R8b9Ff6Ljf43+RMf/Ev2Jjv85+hMd/1P0Jzr+x+hPdPwP0Z/o+O+jP9Hx30V/ouO/jf5Ex38T/YmOnxH9iY7/OvoTHf9V9Cc6/svoT3T8F9Gf6Pjp0Z/o+M+jP9Hxn0V/ouM/jf5Ex38S/YmO/zj6Ex3/UfQnOv7DOHXo+A+iP9Hx70d/ouPfi/5Ex78b/YmOfyf6Ex3/dvQnOv6t6E90/JvRn+j4N6I/0fGvR3+i41+L/kTHvxr9iY5/JfoTHf9y9Cc6/qXoT3T8i9Gf6PgXoj/R8c9Hf6Ljn4v+RMc/G/2Jjp8W/YmOfyb6GR3/dPwXoR4BOA==AQAAAAAAAAAAgAAAAAAAAOAEAAAAAAAADgEAAAAAAAA=eF4txRFwAgAAAMC2C4IgCIIgCIIgCIIgCIIgCIIgCIIgCIIgCIIgCIIgCIIgCAaDwSAIgiAIgiAIBkEQDPqXDwbeQg474qhjjjvhpFNOO+Osc8674KJLLrviqmuuu+GmW26746577nvgoUcee+KpZ5574aVXXnvjrb/87R//eue9Dz765LMvvvrmu//88NMvBz7eBR1y2BFHHXPcCSedctoZZ51z3gUXXXLZFVddc90NN91y2x133XPfAw898tgTTz3z3AsvvfLaG2/95W//+Nc7733w0SefffHVN9/954effjnw+S7okMOOOOqY40446ZTTzjjrnPMuuOiSy6646prrbrjpltvu+B9fwUXT + + diff --git a/geos-mesh/tests/data/fracture_res5_id_empty.vtp b/geos-mesh/tests/data/fracture_res5_id_empty.vtp new file mode 100644 index 000000000..e3b8bddd8 --- /dev/null +++ b/geos-mesh/tests/data/fracture_res5_id_empty.vtp @@ -0,0 +1,54 @@ + + + + + + + + + + + + + + + + + 0 + + + 2304.8861143 + + + + + 0 + + + 2304.8861143 + + + + + + + + + + + + + + + + + + + + + + + + _AQAAAAAAAAAAgAAAAAAAAKAGAAAAAAAAbgEAAAAAAAA=eF4txdciEAAAAEBRUmlpK9q0aEpbey/tTUMb0d57T6VBO+2NSGjvoflDHrp7uYCA/6o40EGu6moOdnWHuIZrupZDXdt1XNf1XN9hbuCGbuTGbuKmbuZwN3cLRzjSLd3Krd3Gbd3O7R3laHdwR3dyZ3dxjGPd1d3c3T3c070c596Odx/3dT/39wAP9CAneLCHeKiHebhHeKRHebTHeKzHebwneKInebITPcVTPc3TPcMzPcuzPcdzPc/zvcBJTvZCL/JiL3GKl3qZl3uFV3qVVzvVaU73Gmc402u9zuu9wRu9yZu9xVu9zdu9wzu9y7u9x3u9z/t9wAd9yId9xEd9zMd9wid9ylk+7TPO9lmf83lfcI5zfdGXfNlXfNXXfN03nOebvuXbvuO7vuf7fuCHfuTHfuKnzneBC/3MRS72c5f4hUtd5nK/9Cu/9hu/9Tu/9wd/9Cd/9hd/9Td/9w9X+Kd/+bf/+K//uRLqf1dfAQAAAAAAAAAAgAAAAAAAAOAEAAAAAAAAEgEAAAAAAAA=eF4txddCCAAAAMAiozSkoaGh0NAeqGhrSNEg7SFkJKFhlIhoaCBUP9tDdy8XEHAo0Ed81EE+5uM+4ZMOdohPOdRhDneETzvSZxzlaMc41mcd53gnONHnnORkpzjV553mdF/wRV9yhjOd5Wxfdo5zned8F7jQRS52iUt9xVd9zWUud4Wv+4YrXeVq17jWda73TTe40U1u9i23+LZb3eY7vut2d7jTXb7n++72A/e4133u94AHPeRhj3jUDz3mR37sJx73Uz/zc7/whF960q885dd+47ee9oxnPed3fu8P/uh5L/iTF/3ZX7zkr/7mZX/3D6941Wte909veNNb3vYv//Yf7/iv//m/d73nfR8ARZMvOw==AQAAAAAAAAAAgAAAAAAAAPAJAAAAAAAAagIAAAAAAAA=eF5t1aFqK1EUheGJqokoEdVVl3mKoSS3b5AHqC4E+gjja6qihxEjbkxUdZM2LpS62FA/UFlG3Z1N9v5XoVEfw5xz1hwW2UVhv2a8K+JXXuD1CD98v6XrHj994uaA13v8ssEfz/i4wl8NLpb48hFf1/juXvLMJUMl55Zy1kT2H16TN++4bsULXFyJl9vi199yiq9myXohbvHNOy4GfJz8Tb+UuKlwPcd39/i6xpePuFjirwYfV/jjWc7d4PVeMhzw06fk6fHDt6wd3abLC9yM1bvVvz+z8zvYupe27qWte2nrXtq6x9o9tu6lrXtp617aupe27qWte2nrXtq6R565ZKjk3FLOmsj+w2vausc+rXiBrXuSbYuLQp5P5f2Z7CNusXVP8mDrnuTH1j35Xmzdk/vB1j3Jhq17cv/Yuse5G2zdI8MBW/fI02PrHmtHt2nrHvuM1daxLt7B3r0u9nxLe/e6yIC9e7F2j717XXwj9u51cSfYu9fFHWLvXhd3LnnmkqGSc0s5ayL7n7p3tncv9mnFC+zdy2xbfOpePp/K+zPZR9xi717mwd69zI+9e/m92LuX94O9e5kNe/fy/rF3L87dYO9eZDhg717k6bF3L9aeune2dy/2Gat3+b/vM3eItdhn7hBnYZ+5Q2TDPnOH+BbsM3eIb8c+c2P2LLHP3CHuFvvMjTxzyVDJuaWcNZH9TzP3bJ+5sU8rXmCfuZlt+2NG8nwq789kH3GLfeZmHuwzN/Njn7n5vdhnbt4P9pmb2bDP3DhrhX3mxrkb7DM3Mhywz9zI02OfubH2NHPP9pkb+4zd/wFdP06cAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAQAAAAAAAAAAgAAAAAAAAIATAAAAAAAAZwMAAAAAAAA=eF5d1cMWGGYABeHUSm3bbmqmtu3Utp3atpvadpsitfVCXXS+Rf+7mXPmzP5OmvDfNo6nxVPjRnHS4E+JGw4df3LcYOj4k+L6Q8efGNcbOv6EuO7Q8cfHdYaOPy6uPXT8lLjW0PHHxjWHjj8mrjF0/NFx9aHjj4qrDR1/ZFx16Pgj4ipDxx8eVx46/rC40tDxh8YVh44/JK4wdPzBcfmh4w+Kyw0df2Bcduj4A+IyQ8fvH5ceOn6/uNTQ8fvGJYeO3ycuMXT83nHxoeP3iosNHb9nXHTo+D3iIkPH7x4XHjp+t7jQ0PG7xgWHjt8lLjB0/M5x/qHjd4rzDR2/Y5x36Pgd4sSh47eP8wwdPznOPXT8dnGuoeO3jXMOHb9NnGPo+K3j7EPHbxVnGzp+yzjr0PFbxFmGjt88zjx0/GZxpqHjN53w/+n4TaI/8TtPxSejP9HxT0R/ouMfj/5Exz8W/YmOfzT6Ex3/SPQnOv7h6E90/ENxytDxD0Z/ouMfiP5Ex98f/YmOvy/6Ex1/b/QnOv6e6E90/N3Rn+j4u6I/0fF3Rn+i4++I/kTH3x79iY6/LfoTHX9r9Cc6/pboT3T8zdGf6Piboj/R8TdGf6Ljb4j+RMdfH/2Jjr8u+hMdf230Jzp+avQnOv6a6E90/NXRn+j4q6I/0fFXRn+i46+I/kTHXx79iY6/LPoTHX9p9Cc6/pLoT3T8xXHy0PEXRX+i4y+M/kTHXxD9iY4/P/oTHX9e9Cc6/tzoT3T8OdGf6Pizoz/R8WdFf6Ljz4z+RMefEf2Mjj89+hO/80/8O/oTHf9X9Cc6/s/oT3T8H9Gf6Pjfoz/R8b9Ff6Ljf43+RMf/Ev2Jjv85+hMd/1P0Jzr+x+hPdPwP0Z/o+O+jP9Hx30V/ouO/jf5Ex38T/YmOnxH9iY7/OvoTHf9V9Cc6/svoT3T8F9Gf6Pjp0Z/o+M+jP9Hxn0V/ouM/jf5Ex38S/YmO/zj6Ex3/UfQnOv7DOHXo+A+iP9Hx70d/ouPfi/5Ex78b/YmOfyf6Ex3/dvQnOv6t6E90/JvRn+j4N6I/0fGvR3+i41+L/kTHvxr9iY5/JfoTHf9y9Cc6/qXoT3T8i9Gf6PgXoj/R8c9Hf6Ljn4v+RMc/G/2Jjp8W/YmOfyb6GR3/dPwXoR4BOA==AQAAAAAAAAAAgAAAAAAAAOAEAAAAAAAADgEAAAAAAAA=eF4txRFwAgAAAMC2C4IgCIIgCIIgCIIgCIIgCIIgCIIgCIIgCIIgCIIgCIIgCAaDwSAIgiAIgiAIBkEQDPqXDwbeQg474qhjjjvhpFNOO+Osc8674KJLLrviqmuuu+GmW26746577nvgoUcee+KpZ5574aVXXnvjrb/87R//eue9Dz765LMvvvrmu//88NMvBz7eBR1y2BFHHXPcCSedctoZZ51z3gUXXXLZFVddc90NN91y2x133XPfAw898tgTTz3z3AsvvfLaG2/95W//+Nc7733w0SefffHVN9/954effjnw+S7okMOOOOqY40446ZTTzjjrnPMuuOiSy6646prrbrjpltvu+B9fwUXT + + diff --git a/geos-mesh/tests/test_AttributeMapping.py b/geos-mesh/tests/test_AttributeMapping.py new file mode 100644 index 000000000..019492f18 --- /dev/null +++ b/geos-mesh/tests/test_AttributeMapping.py @@ -0,0 +1,35 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Romain Baville +# ruff: noqa: E402 # disable Module level import not at top of file +import pytest + +from typing import Union, Any +from geos.mesh.utils.arrayModifiers import fillAllPartialAttributes +from geos.mesh.processing.AttributeMapping import AttributeMapping + +from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet, vtkDataSet + + +@pytest.mark.parametrize( "meshFromName, meshToName, attributeNames, onPoints", [ + ( "fracture", "emptyFracture", { "collocated_nodes" }, True ), + ( "multiblock", "emptyFracture", { "FAULT" }, False ), + ( "multiblock", "emptymultiblock", { "FAULT" }, False ), + ( "dataset", "emptymultiblock", { "FAULT" }, False ), + ( "dataset", "emptydataset", { "FAULT" }, False ), +] ) +def test_AttributeMapping( + dataSetTest: Any, + meshFromName: str, + meshToName: str, + attributeNames: set[ str ], + onPoints: bool, +) -> None: + """Test mapping an attribute between two meshes.""" + meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet ] = dataSetTest( meshFromName ) + meshTo: Union[ vtkDataSet, vtkMultiBlockDataSet ] = dataSetTest( meshToName ) + if isinstance( meshFrom, vtkMultiBlockDataSet ): + fillAllPartialAttributes( meshFrom ) + + filter = AttributeMapping( meshFrom, meshTo, attributeNames, onPoints ) + assert filter.applyFilter() diff --git a/geos-mesh/tests/test_arrayHelpers.py b/geos-mesh/tests/test_arrayHelpers.py index cbfaf50d4..e1d58ab38 100644 --- a/geos-mesh/tests/test_arrayHelpers.py +++ b/geos-mesh/tests/test_arrayHelpers.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 # mypy: disable-error-code="operator, attr-defined" @@ -19,6 +19,38 @@ from geos.mesh.utils.arrayModifiers import createConstantAttribute +@pytest.mark.parametrize( "meshFromName, meshToName, points", [ + ( "multiblock", "emptymultiblock", False ), + ( "multiblock", "emptyFracture", False ), + ( "dataset", "emptyFracture", False ), + ( "dataset", "emptypolydata", False ), + ( "fracture", "emptyFracture", True ), + ( "fracture", "emptyFracture", False ), + ( "fracture", "emptymultiblock", False ), + ( "polydata", "emptypolydata", False ), +] ) +def test_computeElementMapping( + dataSetTest: vtkDataSet, + getElementMap: dict[ int, npt.NDArray[ np.int64 ] ], + meshFromName: str, + meshToName: str, + points: bool, +) -> None: + """Test getting the map between two meshes element.""" + meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet ] = dataSetTest( meshFromName ) + meshTo: Union[ vtkDataSet, vtkMultiBlockDataSet ] = dataSetTest( meshToName ) + elementMapComputed: dict[ int, npt.NDArray[ np.int64 ] ] = arrayHelpers.computeElementMapping( + meshFrom, meshTo, points ) + elementMapTest: dict[ int, npt.NDArray[ np.int64 ] ] = getElementMap( meshFromName, meshToName, points ) + + keysComputed: list[ int ] = list( elementMapComputed.keys() ) + keysTest: list[ int ] = list( elementMapTest.keys() ) + assert keysComputed == keysTest + + for key in keysTest: + assert np.all( elementMapComputed[ key ] == elementMapTest[ key ] ) + + @pytest.mark.parametrize( "onpoints, expected", [ ( True, { 'GLOBAL_IDS_POINTS': 1, 'collocated_nodes': 2, @@ -274,18 +306,14 @@ def test_getComponentNamesMultiBlock( assert obtained == expected -@pytest.mark.parametrize( "attributeNames, expected_columns", [ - ( ( "CellAttribute1", ), ( "CellAttribute1_0", "CellAttribute1_1", "CellAttribute1_2" ) ), - ( ( - "CellAttribute1", - "CellAttribute2", - ), ( "CellAttribute2", "CellAttribute1_0", "CellAttribute1_1", "CellAttribute1_2" ) ), +@pytest.mark.parametrize( "attributeNames, onPoints, expected_columns", [ + ( ( "collocated_nodes", ), True, ( "collocated_nodes_0", "collocated_nodes_1" ) ), ] ) -def test_getAttributeValuesAsDF( dataSetTest: vtkPolyData, attributeNames: Tuple[ str, ...], +def test_getAttributeValuesAsDF( dataSetTest: vtkPolyData, attributeNames: Tuple[ str, ...], onPoints: bool, expected_columns: Tuple[ str, ...] ) -> None: """Test getting an attribute from a polydata as a dataframe.""" - polydataset: vtkPolyData = dataSetTest( "polydata" ) - data: pd.DataFrame = arrayHelpers.getAttributeValuesAsDF( polydataset, attributeNames ) + polydataset: vtkPolyData = vtkPolyData.SafeDownCast( dataSetTest( "polydata" ) ) + data: pd.DataFrame = arrayHelpers.getAttributeValuesAsDF( polydataset, attributeNames, onPoints ) obtained_columns = data.columns.values.tolist() assert obtained_columns == list( expected_columns ) diff --git a/geos-mesh/tests/test_arrayModifiers.py b/geos-mesh/tests/test_arrayModifiers.py index f4446de6f..e03a5b57f 100644 --- a/geos-mesh/tests/test_arrayModifiers.py +++ b/geos-mesh/tests/test_arrayModifiers.py @@ -481,6 +481,59 @@ def test_copyAttributeDataSet( assert vtkDataTypeCopied == vtkDataTypeTest +@pytest.mark.parametrize( "meshFromName, meshToName, attributeName, onPoints, defaultValueTest", [ + ( "fracture", "emptyFracture", "collocated_nodes", True, [ -1, -1 ] ), + ( "multiblock", "emptyFracture", "FAULT", False, -1 ), + ( "multiblock", "emptymultiblock", "FAULT", False, -1 ), + ( "dataset", "emptymultiblock", "FAULT", False, -1 ), + ( "dataset", "emptydataset", "FAULT", False, -1 ), +] ) +def test_transferAttributeWithElementMap( + dataSetTest: Any, + getElementMap: dict[ int, npt.NDArray[ np.int64 ] ], + meshFromName: str, + meshToName: str, + attributeName: str, + onPoints: bool, + defaultValueTest: Any, +) -> None: + """Test to transfer attributes from the source mesh to the final mesh using a map of points/cells.""" + meshFrom: Union[ vtkMultiBlockDataSet, vtkDataSet ] = dataSetTest( meshFromName ) + if isinstance( meshFrom, vtkMultiBlockDataSet ): + arrayModifiers.fillAllPartialAttributes( meshFrom ) + + meshTo: Union[ vtkMultiBlockDataSet, vtkDataSet ] = dataSetTest( meshToName ) + elementMap: dict[ int, npt.NDArray[ np.int64 ] ] = getElementMap( meshFromName, meshToName, onPoints ) + + assert arrayModifiers.transferAttributeWithElementMap( meshFrom, meshTo, elementMap, attributeName, onPoints ) + + for flatIdDataSetTo in elementMap: + dataTo: Union[ vtkPointData, vtkCellData ] + if isinstance( meshTo, vtkDataSet ): + dataTo = meshTo.GetPointData() if onPoints else meshTo.GetCellData() + elif isinstance( meshTo, vtkMultiBlockDataSet ): + dataSetTo: vtkDataSet = vtkDataSet.SafeDownCast( meshTo.GetDataSet( flatIdDataSetTo ) ) + dataTo = dataSetTo.GetPointData() if onPoints else dataSetTo.GetCellData() + + arrayTo: npt.NDArray[ Any ] = vnp.vtk_to_numpy( dataTo.GetArray( attributeName ) ) + for idElementTo in range( len( arrayTo ) ): + idElementFrom: int = elementMap[ flatIdDataSetTo ][ idElementTo ][ 1 ] + if idElementFrom == -1: + assert arrayTo[ idElementTo ] == defaultValueTest + + else: + dataFrom: Union[ vtkPointData, vtkCellData ] + if isinstance( meshFrom, vtkDataSet ): + dataFrom = meshFrom.GetPointData() if onPoints else meshFrom.GetCellData() + elif isinstance( meshFrom, vtkMultiBlockDataSet ): + flatIdDataSetFrom: int = elementMap[ flatIdDataSetTo ][ idElementTo ][ 0 ] + dataSetFrom: vtkDataSet = vtkDataSet.SafeDownCast( meshFrom.GetDataSet( flatIdDataSetFrom ) ) + dataFrom = dataSetFrom.GetPointData() if onPoints else dataSetFrom.GetCellData() + + arrayFrom: npt.NDArray[ Any ] = vnp.vtk_to_numpy( dataFrom.GetArray( attributeName ) ) + assert np.all( arrayTo[ idElementTo ] == arrayFrom[ idElementFrom ] ) + + @pytest.mark.parametrize( "attributeName, onPoints", [ ( "CellAttribute", False ), ( "PointAttribute", True ), diff --git a/geos-posp/src/PVplugins/PVAttributeMapping.py b/geos-posp/src/PVplugins/PVAttributeMapping.py deleted file mode 100644 index 39b17b51f..000000000 --- a/geos-posp/src/PVplugins/PVAttributeMapping.py +++ /dev/null @@ -1,279 +0,0 @@ -# SPDX-License-Identifier: Apache-2.0 -# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -# SPDX-FileContributor: Raphaël Vinour, Martin Lemay -# ruff: noqa: E402 # disable Module level import not at top of file -import os -import sys -from typing import Union - -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_posp.filters.AttributeMappingFromCellCoords import ( - AttributeMappingFromCellCoords, ) -from geos.mesh.utils.arrayModifiers import fillPartialAttributes -from geos.mesh.utils.multiblockModifiers import mergeBlocks -from geos.mesh.utils.arrayHelpers import ( - getAttributeSet, ) -from geos_posp.visu.PVUtils.checkboxFunction import ( # type: ignore[attr-defined] - createModifiedCallback, ) -from geos_posp.visu.PVUtils.paraviewTreatments import getArrayChoices -from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] - VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, -) -from vtkmodules.vtkCommonCore import ( - vtkDataArraySelection, - vtkInformation, - vtkInformationVector, -) -from vtkmodules.vtkCommonDataModel import ( - vtkCompositeDataSet, - vtkDataObjectTreeIterator, - vtkDataSet, - vtkMultiBlockDataSet, - vtkUnstructuredGrid, -) - -__doc__ = """ -Map the attributes from a source mesh to a client mesh. - -Input and output are vtkUnstructuredGrid. - -To use it: - -* Load the module in Paraview: Tools>Manage Plugins...>Load new>PVAttributeMapping. -* Select the server mesh. -* Search and Select Attribute Mapping Filter. -* Select the client mesh. -* Select the attributes to transfer and Apply. - -""" - - -@smproxy.filter( name="PVAttributeMapping", label="Attribute Mapping" ) -@smhint.xml( '' ) -@smproperty.input( name="Client", port_index=1, label="Client mesh" ) -@smdomain.datatype( - dataTypes=[ "vtkUnstructuredGrid", "vtkMultiBlockDataSet" ], - composite_data_supported=True, -) -@smproperty.input( name="Server", port_index=0, label="Server mesh" ) -@smdomain.datatype( - dataTypes=[ "vtkUnstructuredGrid" ], - composite_data_supported=False, -) -class PVAttributeMapping( VTKPythonAlgorithmBase ): - - def __init__( self: Self ) -> None: - """Map the properties of a server mesh to a client mesh.""" - super().__init__( nInputPorts=2, nOutputPorts=1, outputType="vtkUnstructuredGrid" ) - - # boolean to check if first use of the filter for attribute list initialization - self.m_firstUse = True - - # list of attribute names to transfer - self.m_attributes: vtkDataArraySelection = vtkDataArraySelection() - self.m_attributes.AddObserver( 0, createModifiedCallback( self ) ) - - # logger - self.m_logger: Logger = getLogger( "Attribute Mapping" ) - - def SetLogger( self: Self, logger: Logger ) -> None: - """Set filter logger. - - Args: - logger (Logger): logger - """ - self.m_logger = logger - self.Modified() - - @smproperty.dataarrayselection( name="AttributesToTransfer" ) - def a02GetAttributeToTransfer( self: Self ) -> vtkDataArraySelection: - """Get selected attribute names to transfer. - - Returns: - vtkDataArraySelection: selected attribute names. - """ - return self.m_attributes - - 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. - """ - # only at initialization step, no change later - if self.m_firstUse: - # get cell ids - inData = self.GetInputData( inInfoVec, 0, 0 ) - assert isinstance( inData, ( vtkDataSet, vtkMultiBlockDataSet ) ), "Input object type is not supported." - - # update vtkDAS - attributeNames: set[ str ] = getAttributeSet( inData, False ) - for attributeName in attributeNames: - if not self.m_attributes.ArrayExists( attributeName ): - self.m_attributes.AddArray( attributeName ) - - self.m_firstUse = False - return 1 - - def RequestDataObject( - self: Self, - request: vtkInformation, - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestDataObject. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - inData = self.GetInputData( inInfoVec, 1, 0 ) - outData = self.GetOutputData( outInfoVec, 0 ) - assert inData is not None - if outData is None or ( not outData.IsA( inData.GetClassName() ) ): - outData = inData.NewInstance() - outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) - return super().RequestDataObject( request, inInfoVec, outInfoVec ) # type: ignore[no-any-return] - - 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: - serverMesh: Union[ vtkUnstructuredGrid, vtkMultiBlockDataSet, - vtkCompositeDataSet ] = self.GetInputData( inInfoVec, 0, 0 ) - clientMesh: Union[ vtkUnstructuredGrid, vtkMultiBlockDataSet, - vtkCompositeDataSet ] = self.GetInputData( inInfoVec, 1, 0 ) - outData: Union[ vtkUnstructuredGrid, vtkMultiBlockDataSet, - vtkCompositeDataSet ] = self.GetOutputData( outInfoVec, 0 ) - - assert serverMesh is not None, "Input server mesh is null." - assert clientMesh is not None, "Input client mesh is null." - assert outData is not None, "Output pipeline is null." - - outData.ShallowCopy( clientMesh ) - attributeNames: set[ str ] = set( getArrayChoices( self.a02GetAttributeToTransfer() ) ) - for attributeName in attributeNames: - fillPartialAttributes( serverMesh, attributeName, False ) - - mergedServerMesh: vtkUnstructuredGrid - if isinstance( serverMesh, vtkUnstructuredGrid ): - mergedServerMesh = serverMesh - elif isinstance( serverMesh, ( vtkMultiBlockDataSet, vtkCompositeDataSet ) ): - mergedServerMesh = mergeBlocks( serverMesh ) - else: - raise ValueError( "Server mesh data type is not supported. " + - "Use either vtkUnstructuredGrid or vtkMultiBlockDataSet" ) - - if isinstance( outData, vtkUnstructuredGrid ): - self.doTransferAttributes( mergedServerMesh, outData, attributeNames ) - elif isinstance( outData, ( vtkMultiBlockDataSet, vtkCompositeDataSet ) ): - self.transferAttributesMultiBlockDataSet( mergedServerMesh, outData, attributeNames ) - else: - raise ValueError( "Client mesh data type is not supported. " + - "Use either vtkUnstructuredGrid or vtkMultiBlockDataSet" ) - - outData.Modified() - - mess: str = "Attributes were successfully transferred ." - self.m_logger.info( mess ) - except AssertionError as e: - mess1: str = "Attribute transfer failed due to:" - self.m_logger.error( mess1 ) - self.m_logger.error( e, exc_info=True ) - return 0 - except Exception as e: - mess0: str = "Attribute transfer failed due to:" - self.m_logger.critical( mess0 ) - self.m_logger.critical( e, exc_info=True ) - return 0 - return 1 - - def doTransferAttributes( - self: Self, - serverMesh: vtkUnstructuredGrid, - clientMesh: vtkUnstructuredGrid, - attributeNames: set[ str ], - ) -> bool: - """Transfer attributes between two vtkUnstructuredGrids. - - Args: - serverMesh (vtkUnstructuredGrid): server mesh - clientMesh (vtkUnstructuredGrid): client mesh - attributeNames (set[str]): set of attribut names to transfer - - Returns: - bool: True if attributes were successfully transferred. - - """ - filter: AttributeMappingFromCellCoords = AttributeMappingFromCellCoords() - filter.AddInputDataObject( 0, serverMesh ) - filter.AddInputDataObject( 1, clientMesh ) - filter.SetTransferAttributeNames( attributeNames ) - filter.Update() - clientMesh.ShallowCopy( filter.GetOutputDataObject( 0 ) ) - return True - - def transferAttributesMultiBlockDataSet( - self: Self, - serverBlock: vtkUnstructuredGrid, - clientMesh: Union[ vtkMultiBlockDataSet, vtkCompositeDataSet ], - attributeNames: set[ str ], - ) -> bool: - """Transfer attributes from a vtkUnstructuredGrid to a multiblock. - - Args: - serverBlock (vtkUnstructuredGrid): server mesh - clientMesh (vtkMultiBlockDataSet | vtkCompositeDataSet): client mesh - attributeNames (set[str]): set of attribut names to transfer - - Returns: - bool: True if attributes were successfully transferred. - - """ - iter: vtkDataObjectTreeIterator = vtkDataObjectTreeIterator() - iter.SetDataSet( clientMesh ) - iter.VisitOnlyLeavesOn() - iter.GoToFirstItem() - while iter.GetCurrentDataObject() is not None: - clientBlock: vtkUnstructuredGrid = vtkUnstructuredGrid.SafeDownCast( iter.GetCurrentDataObject() ) - self.doTransferAttributes( serverBlock, clientBlock, attributeNames ) - clientBlock.Modified() - iter.GoToNextItem() - return True diff --git a/geos-posp/src/PVplugins/PVTransferAttributesVolumeSurface.py b/geos-posp/src/PVplugins/PVTransferAttributesVolumeSurface.py deleted file mode 100644 index bbc5c1ad3..000000000 --- a/geos-posp/src/PVplugins/PVTransferAttributesVolumeSurface.py +++ /dev/null @@ -1,273 +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 - -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_posp.filters.TransferAttributesVolumeSurface import ( - TransferAttributesVolumeSurface, ) -from geos.mesh.utils.multiblockHelpers import ( - getBlockElementIndexesFlatten, - getBlockFromFlatIndex, -) -from geos.mesh.utils.arrayHelpers import getAttributeSet -from geos.mesh.utils.multiblockModifiers import mergeBlocks -from geos_posp.visu.PVUtils.checkboxFunction import ( # type: ignore[attr-defined] - createModifiedCallback, ) -from geos_posp.visu.PVUtils.paraviewTreatments import getArrayChoices -from paraview.simple import ( # type: ignore[import-not-found] - FindSource, GetActiveSource, GetSources, servermanager, -) -from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] - VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, -) -from vtkmodules.vtkCommonCore import ( - vtkDataArray, - vtkInformation, - vtkInformationVector, -) -from vtkmodules.vtkCommonCore import ( - vtkDataArraySelection as vtkDAS, ) -from vtkmodules.vtkCommonDataModel import ( - vtkDataObject, - vtkMultiBlockDataSet, - vtkPolyData, - vtkUnstructuredGrid, -) - -__doc__ = """ -PVTransferAttributesVolumeSurface is a Paraview plugin that allows to map face ids from -surface mesh with cell ids from volume mesh. - -Input and output types are vtkMultiBlockDataSet. - -To use it: - -* Load the module in Paraview: Tools>Manage Plugins...>Load new>PVTransferAttributesVolumeSurface. -* Select first a volume mesh, then a surface mesh, both of type vtkMultiBlockDataSet -* Search and Apply 'Geos Volume Surface Mesh Mapper' Filter. - -""" - - -@smproxy.filter( - name="PVTransferAttributesVolumeSurface", - label="Geos Transfer Attributes From Volume to Surface", -) -@smhint.xml( '' ) -@smproperty.input( name="SurfaceMesh", port_index=1 ) -@smdomain.datatype( - dataTypes=[ "vtkMultiBlockDataSet", "vtkUnstructuredGrid", "vtkPolyData" ], - composite_data_supported=True, -) -@smproperty.input( name="VolumeMesh", port_index=0 ) -@smdomain.datatype( - dataTypes=[ "vtkMultiBlockDataSet", "vtkUnstructuredGrid" ], - composite_data_supported=True, -) -class PVTransferAttributesVolumeSurface( 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=2, nOutputPorts=1, outputType="vtkMultiBlockDataSet" ) - - #: input volume mesh - self.m_volumeMesh: vtkMultiBlockDataSet - #: output surface mesh - self.m_outputSurfaceMesh: vtkMultiBlockDataSet - # volume mesh source - self.m_sourceNameVolume: str = [ k for ( k, v ) in GetSources().items() if v == GetActiveSource() ][ 0 ][ 0 ] - - # list of attributes - self.m_attributes: vtkDAS = self.createAttributeVTKDas() - # logger - self.m_logger: Logger = getLogger( "Volume Surface Mesh Mapper Filter" ) - - def SetLogger( self: Self, logger: Logger ) -> None: - """Set filter logger. - - Args: - logger (Logger): logger - """ - self.m_logger = logger - - def createAttributeVTKDas( self: Self ) -> vtkDAS: - """Create the vtkDataArraySelection with cell attribute names. - - Returns: - vtkDAS: vtkDataArraySelection with attribute names - """ - source = FindSource( self.m_sourceNameVolume ) - dataset = servermanager.Fetch( source ) - attributes: set[ str ] = getAttributeSet( dataset, False ) - attributesDAS: vtkDAS = vtkDAS() - for attribute in attributes: - attributesDAS.AddArray( attribute, False ) - attributesDAS.AddObserver( "ModifiedEvent", createModifiedCallback( self ) ) # type: ignore[arg-type] - return attributesDAS - - @smproperty.dataarrayselection( name="AttributesToTransfer" ) - def a02GetAttributeToTransfer( self: Self ) -> vtkDAS: - """Get selected attribute names to transfer. - - Returns: - vtkDataArraySelection: selected attribute names. - """ - return self.m_attributes - - 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" ) - else: - 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 RequestDataObject( - self: Self, - request: vtkInformation, - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestDataObject. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - inData1 = self.GetInputData( inInfoVec, 0, 0 ) - inData2 = self.GetInputData( inInfoVec, 1, 0 ) - outData = self.GetOutputData( outInfoVec, 0 ) - assert inData1 is not None - assert inData2 is not None - if outData is None or ( not outData.IsA( inData2.GetClassName() ) ): - outData = inData2.NewInstance() - outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) - return super().RequestDataObject( request, inInfoVec, outInfoVec ) # type: ignore[no-any-return] - - 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: - self.m_logger.info( f"Apply filter {__name__}" ) - self.SetProgressText( "Computation in progress..." ) - self.SetProgressShiftScale( 50, 100 ) - - self.m_volumeMesh = vtkMultiBlockDataSet.GetData( inInfoVec[ 0 ] ) - surfaceMesh: vtkMultiBlockDataSet = vtkMultiBlockDataSet.GetData( inInfoVec[ 1 ] ) - self.m_outputSurfaceMesh = self.GetOutputData( outInfoVec, 0 ) - - assert self.m_volumeMesh is not None, "Input Volume mesh is null." - assert surfaceMesh is not None, "Input Surface mesh is null." - assert self.m_outputSurfaceMesh is not None, "Output pipeline is null." - - self.m_outputSurfaceMesh.ShallowCopy( surfaceMesh ) - self.transferAttributes() - - self.m_logger.info( "Attributes were successfully transfered." ) - except AssertionError as e: - mess: str = "Attribute transfer failed due to:" - self.m_logger.error( mess ) - self.m_logger.error( e, exc_info=True ) - return 0 - except Exception as e: - mess0: str = "Attribute transfer failed due to:" - self.m_logger.critical( mess0 ) - self.m_logger.critical( e, exc_info=True ) - return 0 - return 1 - - def transferAttributes( self: Self ) -> bool: - """Do transfer attributes from volume to surface. - - Returns: - bool: True if transfer is successfull, False otherwise. - """ - attributeNames: set[ str ] = set( getArrayChoices( self.a02GetAttributeToTransfer() ) ) - volumeMeshMerged: vtkUnstructuredGrid = mergeBlocks( self.m_volumeMesh ) - surfaceBlockIndexes: list[ int ] = getBlockElementIndexesFlatten( self.m_outputSurfaceMesh ) - for blockIndex in surfaceBlockIndexes: - surfaceBlock0: vtkDataObject = getBlockFromFlatIndex( self.m_outputSurfaceMesh, blockIndex ) - assert surfaceBlock0 is not None, "surfaceBlock0 is undefined." - surfaceBlock: vtkPolyData = vtkPolyData.SafeDownCast( surfaceBlock0 ) - assert surfaceBlock is not None, "surfaceBlock is undefined." - - assert isinstance( surfaceBlock, vtkPolyData ), "Wrong object type." - - # do transfer of attributes - filter: TransferAttributesVolumeSurface = TransferAttributesVolumeSurface() - filter.AddInputDataObject( 0, volumeMeshMerged ) - filter.AddInputDataObject( 1, surfaceBlock ) - filter.SetAttributeNamesToTransfer( attributeNames ) - filter.Update() - outputSurface: vtkUnstructuredGrid = 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() - - self.m_outputSurfaceMesh.Modified() - return True diff --git a/geos-posp/src/geos_posp/filters/AttributeMappingFromCellCoords.py b/geos-posp/src/geos_posp/filters/AttributeMappingFromCellCoords.py deleted file mode 100644 index 4d9500b11..000000000 --- a/geos-posp/src/geos_posp/filters/AttributeMappingFromCellCoords.py +++ /dev/null @@ -1,235 +0,0 @@ -# SPDX-License-Identifier: Apache-2.0 -# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -# SPDX-FileContributor: Raphaël Vinour, Martin Lemay -# ruff: noqa: E402 # disable Module level import not at top of file - -from collections.abc import MutableSequence - -import numpy as np -import numpy.typing as npt -from geos.utils.Logger import Logger, getLogger -from typing_extensions import Self -from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase -from vtkmodules.vtkCommonCore import ( - vtkDataArray, - vtkDoubleArray, - vtkInformation, - vtkInformationVector, -) -from vtkmodules.vtkCommonDataModel import ( - vtkCellData, - vtkCellLocator, - vtkUnstructuredGrid, -) -from geos.mesh.utils.arrayModifiers import createEmptyAttribute -from geos.mesh.utils.arrayHelpers import ( getVtkArrayInObject, computeCellCenterCoordinates ) - -__doc__ = """ -AttributeMappingFromCellCoords module is a vtk filter that map two identical mesh (or a mesh is -an extract from the other one) and create an attribute containing shared cell ids. - -Filter input and output types are vtkUnstructuredGrid. - -To use the filter: - -.. code-block:: python - - from filters.AttributeMappingFromCellCoords import AttributeMappingFromCellCoords - - # filter inputs - logger :Logger - input :vtkUnstructuredGrid - TransferAttributeName : str - - # instanciate the filter - filter :AttributeMappingFromCellCoords = AttributeMappingFromCellCoords() - # set the logger - filter.SetLogger(logger) - # set input data object - filter.SetInputDataObject(input) - # set Attribute to transfer - filter.SetTransferAttributeNames(AttributeName) - # set Attribute to compare - filter.SetIDAttributeName(AttributeName) - # do calculations - filter.Update() - # get output object - output :vtkPolyData = filter.GetOutputDataObject(0) - # get created attribute names - newAttributeNames :set[str] = filter.GetNewAttributeNames() -""" - - -class AttributeMappingFromCellCoords( VTKPythonAlgorithmBase ): - - def __init__( self: Self ) -> None: - """Map the properties of a source mesh to a receiver mesh.""" - super().__init__( nInputPorts=2, nOutputPorts=1, outputType="vtkUnstructuredGrid" ) - - # source mesh - self.m_clientMesh: vtkUnstructuredGrid - # input mesh - self.m_serverMesh: vtkUnstructuredGrid - # cell map - self.m_cellMap: npt.NDArray[ np.int64 ] = np.empty( 0 ).astype( int ) - - # Transfer Attribute name - self.m_transferedAttributeNames: set[ str ] = set() - # logger - self.m_logger: Logger = getLogger( "Attribute Mapping From Cell Coords Filter" ) - - def SetTransferAttributeNames( self: Self, transferredAttributeNames: set[ str ] ) -> None: - """Setter for transferredAttributeName. - - Args: - transferredAttributeNames (set[str]): set of names of the - attributes to transfer. - - """ - self.m_transferedAttributeNames.clear() - for name in transferredAttributeNames: - self.m_transferedAttributeNames.add( name ) - - def SetLogger( self: Self, logger: Logger ) -> None: - """Set filter logger. - - Args: - logger (Logger): logger - """ - self.m_logger = logger - self.Modified() - - def GetCellMap( self: Self ) -> npt.NDArray[ np.int64 ]: - """Getter of cell map.""" - return self.m_cellMap - - def RequestDataObject( - self: Self, - request: vtkInformation, - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestDataObject. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - inData = self.GetInputData( inInfoVec, 0, 0 ) - outData = self.GetOutputData( outInfoVec, 0 ) - assert inData is not None - if outData is None or ( not outData.IsA( inData.GetClassName() ) ): - outData = inData.NewInstance() - outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) - return super().RequestDataObject( request, inInfoVec, outInfoVec ) # type: ignore[no-any-return] - - 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: - self.m_serverMesh = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) - clientMesh = vtkUnstructuredGrid.GetData( inInfoVec[ 1 ] ) - self.m_clientMesh = self.GetOutputData( outInfoVec, 0 ) - - assert self.m_serverMesh is not None, "Server mesh is null." - assert clientMesh is not None, "Client mesh is null." - assert self.m_clientMesh is not None, "Output pipeline is null." - - self.m_clientMesh.ShallowCopy( clientMesh ) - - # create cell map - self.computeCellMapping() - - # transfer attributes if at least one corresponding cell - if np.any( self.m_cellMap > -1 ): - self.transferAttributes() - self.m_clientMesh.Modified() - else: - self.m_logger.warning( "Input and output meshes do not have any corresponding cells" ) - - except AssertionError as e: - mess1: str = "Mapping Transfer Coord failed due to:" - self.m_logger.error( mess1 ) - self.m_logger.error( e, exc_info=True ) - return 0 - except Exception as e: - mess0: str = "Mapping Transfer Coord failed due to:" - self.m_logger.critical( mess0 ) - self.m_logger.critical( e, exc_info=True ) - return 0 - mess2: str = "Mapping Transfer Coord were successfully computed." - self.m_logger.info( mess2 ) - - return 1 - - def computeCellMapping( self: Self ) -> bool: - """Create the cell map from client to server mesh cell indexes. - - For each cell index of the client mesh, stores the index of the cell - in the server mesh. - - Returns: - bool: True if the map was computed. - - """ - self.m_cellMap = np.full( self.m_clientMesh.GetNumberOfCells(), -1 ).astype( int ) - cellLocator: vtkCellLocator = vtkCellLocator() - cellLocator.SetDataSet( self.m_serverMesh ) - cellLocator.BuildLocator() - - cellCenters: vtkDataArray = computeCellCenterCoordinates( self.m_clientMesh ) - for i in range( self.m_clientMesh.GetNumberOfCells() ): - cellCoords: MutableSequence[ float ] = [ 0.0, 0.0, 0.0 ] - cellCenters.GetTuple( i, cellCoords ) - cellIndex: int = cellLocator.FindCell( cellCoords ) - self.m_cellMap[ i ] = cellIndex - return True - - def transferAttributes( self: Self ) -> bool: - """Transfer attributes from server to client meshes using cell mapping. - - Returns: - bool: True if transfer successfully ended. - - """ - for attributeName in self.m_transferedAttributeNames: - array: vtkDoubleArray = getVtkArrayInObject( self.m_serverMesh, attributeName, False ) - - dataType = array.GetDataType() - nbComponents: int = array.GetNumberOfComponents() - componentNames: list[ str ] = [] - if nbComponents > 1: - for i in range( nbComponents ): - componentNames.append( array.GetComponentName( i ) ) - newArray: vtkDataArray = createEmptyAttribute( self.m_clientMesh, attributeName, tuple( componentNames ), - dataType ) - nanValues: list[ float ] = [ np.nan for _ in range( nbComponents ) ] - for indexClient in range( self.m_clientMesh.GetNumberOfCells() ): - indexServer: int = self.m_cellMap[ indexClient ] - data: MutableSequence[ float ] = nanValues - if indexServer > -1: - array.GetTuple( indexServer, data ) - newArray.InsertNextTuple( data ) - - cellData: vtkCellData = self.m_clientMesh.GetCellData() - assert cellData is not None, "CellData is undefined." - cellData.AddArray( newArray ) - cellData.Modified() - return True diff --git a/geos-posp/src/geos_posp/filters/AttributeMappingFromCellId.py b/geos-posp/src/geos_posp/filters/AttributeMappingFromCellId.py deleted file mode 100644 index ceae63135..000000000 --- a/geos-posp/src/geos_posp/filters/AttributeMappingFromCellId.py +++ /dev/null @@ -1,188 +0,0 @@ -# SPDX-License-Identifier: Apache-2.0 -# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -# SPDX-FileContributor: Raphaël Vinour, Martin Lemay -# ruff: noqa: E402 # disable Module level import not at top of file -import numpy as np -import numpy.typing as npt -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 vtkUnstructuredGrid - -from geos.mesh.utils.arrayModifiers import createAttribute -from geos.mesh.utils.arrayHelpers import getArrayInObject - -__doc__ = """ -AttributeMappingFromCellId module is a vtk filter that transfer a attribute from a -server mesh to a client mesh. - -Filter input and output types are vtkPolyData. - -To use the filter: - -.. code-block:: python - - from filters.AttributeMappingFromCellId import AttributeMappingFromCellId - - # filter inputs - logger :Logger - input :vtkPolyData - TransferAttributeName : str - - # instanciate the filter - filter :AttributeMappingFromCellId = AttributeMappingFromCellId() - # set the logger - filter.SetLogger(logger) - # set input data object - filter.SetInputDataObject(input) - # set Attribute to transfer - filter.SetTransferAttributeName(AttributeName) - # set Attribute to compare - filter.SetIDAttributeName(AttributeName) - # do calculations - filter.Update() - # get output object - output :vtkPolyData = filter.GetOutputDataObject(0) - # get created attribute names - newAttributeNames :set[str] = filter.GetNewAttributeNames() -""" - - -class AttributeMappingFromCellId( VTKPythonAlgorithmBase ): - - def __init__( self: Self ) -> None: - """Map the properties of a source mesh to a receiver mesh.""" - super().__init__( nInputPorts=2, nOutputPorts=1, outputType="vtkUnstructuredGrid" ) - - # Transfer Attribute name - self.m_transferedAttributeName: str = "" - # ID Attribute name - self.m_idAttributeName: str = "" - # logger - self.m_logger: Logger = getLogger( "Attribute Mapping From Cell Id Filter" ) - - def SetLogger( self: Self, logger: Logger ) -> None: - """Set filter logger. - - Args: - logger (Logger): logger - """ - self.m_logger = logger - self.Modified() - - def SetTransferAttributeName( self: Self, name: str ) -> None: - """Set Transfer attribute name.""" - self.m_transferedAttributeName = name - self.Modified() - - def SetIDAttributeName( self: Self, name: str ) -> None: - """Set ID attribute name.""" - self.m_idAttributeName = name - self.Modified() - - def RequestDataObject( - self: Self, - request: vtkInformation, - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestDataObject. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - inData = self.GetInputData( inInfoVec, 0, 0 ) - outData = self.GetOutputData( outInfoVec, 0 ) - assert inData is not None - if outData is None or ( not outData.IsA( inData.GetClassName() ) ): - outData = inData.NewInstance() - outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) - return super().RequestDataObject( request, inInfoVec, outInfoVec ) # type: ignore[no-any-return] - - 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: - serverMesh: vtkUnstructuredGrid = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) - clientMesh: vtkUnstructuredGrid = vtkUnstructuredGrid.GetData( inInfoVec[ 1 ] ) - outData: vtkUnstructuredGrid = self.GetOutputData( outInfoVec, 0 ) - - assert serverMesh is not None, "Server mesh is null." - assert clientMesh is not None, "Client mesh is null." - assert outData is not None, "Output pipeline is null." - - outData.ShallowCopy( clientMesh ) - - cellIdMap: npt.NDArray[ np.int64 ] = self.getCellMap( serverMesh, clientMesh ) - - attributeArrayServer: npt.NDArray[ np.float64 ] = getArrayInObject( serverMesh, - self.m_transferedAttributeName, False ) - attributeArrayClient: npt.NDArray[ np.float64 ] = np.full_like( attributeArrayServer, np.nan ) - for i in range( serverMesh.GetNumberOfCells() ): - k: int = cellIdMap[ i ] - attributeArrayClient[ k ] = attributeArrayServer[ i ] - createAttribute( - clientMesh, - attributeArrayClient, - self.m_transferedAttributeName, - (), - False, - ) - outData.Modified() - - except AssertionError as e: - mess1: str = "Attribute mapping failed due to:" - self.m_logger.error( mess1 ) - self.m_logger.error( e, exc_info=True ) - return 0 - except Exception as e: - mess0: str = "Attribute mapping failed due to:" - self.m_logger.critical( mess0 ) - self.m_logger.critical( e, exc_info=True ) - return 0 - mess2: str = ( f"Attribute {self.m_transferedAttributeName} was successfully transferred." ) - self.m_logger.info( mess2 ) - - return 1 - - def getCellMap( self: Self, serverMesh: vtkUnstructuredGrid, - clientMesh: vtkUnstructuredGrid ) -> npt.NDArray[ np.int64 ]: - """Compute the mapping between both mesh from cell ids. - - Args: - serverMesh (vtkUnstructuredGrid): server mesh - clientMesh (vtkUnstructuredGrid): client mesh - - Returns: - npt.NDArray[np.int64]: the map where for each cell index of the - client mesh, the cell index of the server mesh. - - """ - cellIdArrayServer: npt.NDArray[ np.int64 ] = getArrayInObject( serverMesh, self.m_idAttributeName, - False ).astype( int ) - cellIdArrayClient: npt.NDArray[ np.int64 ] = getArrayInObject( clientMesh, self.m_idAttributeName, - False ).astype( int ) - cellMap: npt.NDArray[ np.int64 ] = np.zeros( serverMesh.GetNumberOfCells() ).astype( int ) - for i, cellId in enumerate( cellIdArrayServer ): - k: int = np.argwhere( cellIdArrayClient == cellId )[ 0 ] - cellMap[ i ] = k - return cellMap diff --git a/geos-posp/src/geos_posp/filters/TransferAttributesVolumeSurface.py b/geos-posp/src/geos_posp/filters/TransferAttributesVolumeSurface.py deleted file mode 100644 index 2640d6250..000000000 --- a/geos-posp/src/geos_posp/filters/TransferAttributesVolumeSurface.py +++ /dev/null @@ -1,306 +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 numpy as np -import numpy.typing as npt -import vtkmodules.util.numpy_support as vnp -from geos.utils.ConnectionSet import ConnectionSetCollection -from geos.utils.GeosOutputsConstants import GeosMeshSuffixEnum -from geos.utils.Logger import Logger, getLogger -from typing_extensions import Self -from vtk import VTK_DOUBLE # type: ignore[import-untyped] -from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase -from vtkmodules.vtkCommonCore import ( - vtkDataArray, - vtkInformation, - vtkInformationVector, -) -from vtkmodules.vtkCommonDataModel import vtkPolyData, vtkUnstructuredGrid - -from geos_posp.filters.VolumeSurfaceMeshMapper import VolumeSurfaceMeshMapper -from geos.mesh.utils.arrayHelpers import ( - getArrayInObject, - getComponentNames, - isAttributeInObject, -) - -__doc__ = """ -TransferAttributesVolumeSurface is a vtk filter that allows to transfer volume -mesh attributes to surface mesh. - -This filter transfer a cell attribute as a face attribute from each volume cell -adjacent to a surface face. Since a face can be adjacent to 2 cells (one at -each side), 2 attributes are created with the suffix '_plus' and '_minus' -depending on the normal vector direction of the face. - -Input and output surface types are vtkPolyData and input volume mesh is -vtkUnstructuredGrid. - -.. WARNING:: - This filter must be use very cautiously since its result may have no sense. - - -To use it: - -.. code-block:: python - - from filters.TransferAttributesVolumeSurface import TransferAttributesVolumeSurface - - filter :TransferAttributesVolumeSurface = TransferAttributesVolumeSurface() - # set input data objects - filter.AddInputDataObject(0, volumeMesh) - filter.AddInputDataObject(1, surfaceMesh) - # set attribute names to transfer - attributeNames :set[str] - filter.SetAttributeNamesToTransfer(attributeNames) - # do calculations - filter.Update() - # get filter output surface with new attributes - output :vtkPolyData = filter.GetOutputDataObject(0) -""" - - -class TransferAttributesVolumeSurface( VTKPythonAlgorithmBase ): - - def __init__( self: Self ) -> None: - """Vtk filter to transfer attributes from volume to surface mesh. - - Input volume is vtkUnstructuredGrid and input, output surface mesh - is vtkPolyData, and the list of names of the attributes to transfer. - """ - super().__init__( nInputPorts=2, nOutputPorts=1, outputType="vtkPolyData" ) - - #: input volume mesh attributes are from - self.m_volumeMesh: vtkUnstructuredGrid - #: input surface mesh where to transfer attributes - self.m_surfaceMesh: vtkPolyData - #: output surface mesh - self.m_outputSurfaceMesh: vtkPolyData - #: set of attribute names to transfer - self.m_attributeNames: set[ str ] = set() - #: create attribute names - self.m_newAttributeNames: set[ str ] = set() - # logger - self.m_logger: Logger = getLogger( "Attribute Transfer from Volume to Surface Filter" ) - - def GetAttributeNamesToTransfer( self: Self ) -> set[ str ]: - """Get the set of attribute names to transfer from volume to surface. - - Returns: - set[str]: set of attributes names. - """ - return self.m_attributeNames - - def SetAttributeNamesToTransfer( self: Self, names: set[ str ] ) -> None: - """Set the set of attribute names to transfer from volume to surface. - - Args: - names (set[str]): set of attributes names. - """ - self.m_attributeNames = names - - def GetNewAttributeNames( self: Self ) -> set[ str ]: - """Get the set of attribute names created in the output surface. - - Returns: - set[str]: Set of new attribute names in the surface. - """ - 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(), "vtkUnstructuredGrid" ) - else: - 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 RequestDataObject( - self: Self, - request: vtkInformation, - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestDataObject. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - inData1 = self.GetInputData( inInfoVec, 0, 0 ) - inData2 = self.GetInputData( inInfoVec, 1, 0 ) - outData = self.GetOutputData( outInfoVec, 0 ) - assert inData1 is not None - assert inData2 is not None - if outData is None or ( not outData.IsA( inData1.GetClassName() ) ): - outData = inData1.NewInstance() - outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) - return super().RequestDataObject( request, inInfoVec, outInfoVec ) # type: ignore[no-any-return] - - 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 volume mesh - self.m_volumeMesh = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) - # input surface - self.m_surfaceMesh = vtkPolyData.GetData( inInfoVec[ 1 ] ) - # output volume mesh - self.m_outputSurfaceMesh = self.GetOutputData( outInfoVec, 0 ) - self.m_outputSurfaceMesh.ShallowCopy( self.m_surfaceMesh ) - - # compute cell adjacency mapping - meshMap: ConnectionSetCollection = self.getMeshMapping() - # do transfer of attributes - self.doTransferAttributes( meshMap ) - self.m_logger.info( "Attribute transfer was successfully computed." ) - except AssertionError as e: - mess: str = "Attribute transfer failed due to:" - self.m_logger.error( mess ) - self.m_logger.error( e, exc_info=True ) - return 0 - except Exception as e: - mess0: str = "Attribute transfer failed due to:" - self.m_logger.critical( mess0 ) - self.m_logger.critical( e, exc_info=True ) - return 0 - return 1 - - def getMeshMapping( self: Self ) -> ConnectionSetCollection: - """Compute cell mapping between volume and surface mesh. - - Returns: - dict[int, dict[int, bool]]: dictionnary of face ids as keys and - volume cell ids and side as values. - """ - filter: VolumeSurfaceMeshMapper = VolumeSurfaceMeshMapper() - filter.AddInputDataObject( 0, self.m_volumeMesh ) - filter.AddInputDataObject( 1, self.m_surfaceMesh ) - filter.SetCreateAttribute( False ) - filter.Update() - return filter.GetSurfaceToVolumeConnectionSets() - - def doTransferAttributes( self: Self, meshMap: ConnectionSetCollection ) -> bool: - """Transfer all attributes from the set of attribute names. - - Except on boundaries, surfaces are bounded by cells along each side. - Two attributes are then created on the surface, one corresponding to - positive side cell values, one corresponding to negative side cell - values. - - Args: - meshMap (dict[int, dict[int, bool]]): map of surface face ids to - volume mesh cell ids and side. - - Returns: - bool: True if transfer successfully ended, False otherwise. - """ - for attributeName in self.m_attributeNames: - # negative side attribute - self.transferAttribute( attributeName, False, meshMap ) - # positive side attribute - self.transferAttribute( attributeName, True, meshMap ) - return True - - def transferAttribute( - self: Self, - attributeName: str, - surfaceSide: bool, - meshMap: ConnectionSetCollection, - ) -> bool: - """Transfer the attribute attributeName from volume to surface mesh. - - Created surface attribute will have the same name as input attribute - name with suffix "_Plus" for positive (True) side and "_Minus" for - negative (False) side. - - Args: - attributeName (str): Name of the attribute to transfer. - surfaceSide (bool): Side of the surface the attribute is from. - meshMap (dict[int, dict[int, bool]]): map of surface face ids to - volume mesh cell ids and side. - - Returns: - bool: True if transfer successfully ended, False otherwise. - """ - # get volume mesh attribute - if isAttributeInObject( self.m_volumeMesh, attributeName, False ): - attr: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_volumeMesh, attributeName, False ) - attrComponentNames: tuple[ str, ...] = getComponentNames( self.m_volumeMesh, attributeName, False ) - # creates attribute arrays on the surface - nbFaces: int = self.m_surfaceMesh.GetNumberOfCells() - nbComponents: int = ( len( attrComponentNames ) if len( attrComponentNames ) > 0 else 1 ) - suffix: str = ( GeosMeshSuffixEnum.SURFACE_PLUS_SUFFIX.value - if surfaceSide else GeosMeshSuffixEnum.SURFACE_MINUS_SUFFIX.value ) - surfaceAttributeName: str = attributeName + suffix - attributeValues: npt.NDArray[ np.float64 ] = np.full( ( nbFaces, nbComponents ), np.nan ) - - # for each face of the surface - for connectionSet in meshMap: - # for each cell of the volume mesh - for cellId, side in connectionSet.getConnectedCellIds().items(): - if side == surfaceSide: - attributeValues[ connectionSet.getCellIdRef() ] = attr[ cellId ] - - surfaceAttribute: vtkDataArray = vnp.numpy_to_vtk( attributeValues, deep=True, array_type=VTK_DOUBLE ) - surfaceAttribute.SetName( surfaceAttributeName ) - if surfaceAttribute.GetNumberOfComponents() > 1: - for i in range( surfaceAttribute.GetNumberOfComponents() ): - surfaceAttribute.SetComponentName( i, attrComponentNames[ i ] ) - - self.m_outputSurfaceMesh.GetCellData().AddArray( surfaceAttribute ) - self.m_outputSurfaceMesh.GetCellData().Modified() - self.m_outputSurfaceMesh.Modified() - - # add new attribute names to the set - self.m_newAttributeNames.add( surfaceAttributeName ) - else: - self.m_logger.warning( f"{attributeName} was skipped since it not" + " in the input volume mesh." ) - - return True diff --git a/geos-posp/src/geos_posp/filters/VolumeSurfaceMeshMapper.py b/geos-posp/src/geos_posp/filters/VolumeSurfaceMeshMapper.py deleted file mode 100644 index c4508bcce..000000000 --- a/geos-posp/src/geos_posp/filters/VolumeSurfaceMeshMapper.py +++ /dev/null @@ -1,423 +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 -from typing import Optional - -import numpy as np -import numpy.typing as npt -import vtkmodules.util.numpy_support as vnp -from geos.utils.ConnectionSet import ( - ConnectionSet, - ConnectionSetCollection, -) -from geos.utils.geometryFunctions import getCellSideAgainstPlane -from geos.utils.GeosOutputsConstants import PostProcessingOutputsEnum -from geos.utils.Logger import Logger, getLogger -from geos.utils.PhysicalConstants import ( - EPSILON, ) -from typing_extensions import Self -from vtk import VTK_INT # type: ignore[import-untyped] -from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase -from vtkmodules.vtkCommonCore import ( - vtkDataArray, - vtkIdList, - vtkInformation, - vtkInformationVector, - vtkPoints, -) -from vtkmodules.vtkCommonDataModel import vtkPolyData, vtkUnstructuredGrid - -__doc__ = """ -VolumeSurfaceMeshMapper is a vtk filter that collects the cell of a volume mesh -adjacents to the faces of a surface mesh. - -VolumeSurfaceMeshMapper inputs are a volume mesh of type vtkUnstructuredGrid -and a surface mesh of type vtkPolyData. - -To use the filter: - -.. code-block:: python - - from filters.VolumeSurfaceMeshMapper import VolumeSurfaceMeshMapper - - # filter inputs - logger :Logger - # input objects - volumeMesh :vtkUnstructuredGrid - surfaceMesh :vtkPolyData - - # instanciate the filter - volumeSurfaceMeshMapper :VolumeSurfaceMeshMapper = VolumeSurfaceMeshMapper() - # set parameters and logger - volumeSurfaceMeshMapper.SetCreateAttribute(True) - volumeSurfaceMeshMapper.SetAttributeName("Surface1_AdjacentCell") - volumeSurfaceMeshMapper.SetLogger(logger) - # set input objects - volumeSurfaceMeshMapper.AddInputDataObject(0, volumeMesh) - volumeSurfaceMeshMapper.AddInputDataObject(1, surfaceMesh) - # do calculations - volumeSurfaceMeshMapper.Update() - # get filter output mesh (same as input volume mesh with new attribute if - # volumeSurfaceMeshMapper.SetCreateAttribute(True)) - output :vtkUnstructuredGrid = volumeSurfaceMeshMapper.GetOutputDataObject(0) - # get filter output mapping - # surface faces to volume cell mapping - surfaceToVolumeMap :dict[int, dict[int, bool]] = volumeSurfaceMeshMapper.GetSurfaceToVolumeMap() - # volume cell to surface faces mapping - volumeToSurfaceMap :dict[int, tuple[set[int], bool]] = volumeSurfaceMeshMapper.GetVolumeToSurfaceMap() -""" - - -class VolumeSurfaceMeshMapper( VTKPythonAlgorithmBase ): - - def __init__( self: Self ) -> None: - """Vtk filter to compute cell adjacency between volume and surface meshes. - - Inputs are vtkUnstructuredGrid for volume mesh, vtkPolyData for - surface mesh, and attribute name. - - Ouputs are the volume mesh (vtkUnstructuredGrid) with a new attribute if - SetCreateAttribute was set to True, and a map of surface face indexess - as keys and adjacent volume cell indexes as values. - """ - super().__init__( nInputPorts=2, nOutputPorts=1, outputType="vtkPolyData" ) - - #: input volume mesh - self.m_volumeMesh: vtkUnstructuredGrid - #: input surface mesh - self.m_surfaceMesh: vtkPolyData - #: cell adjacency mapping that will be computed - # self.m_cellIdMap :dict[int, dict[int, bool]] = {} - self.m_connectionSets: ConnectionSetCollection = ConnectionSetCollection() - #: if set to True, will create an attribute in the volume mesh - self.m_createAttribute = False - #: name of the attribute to create in the volume mesh - self.m_attributeName: str = ( PostProcessingOutputsEnum.ADJACENT_CELL_SIDE.attributeName ) - #: logger - self.m_logger: Logger = getLogger( "Volume to Surface Mapper Filter" ) - - def SetLogger( self: Self, logger: Logger ) -> None: - """Set the logger. - - Args: - logger (Logger): logger - """ - self.m_logger = logger - - def GetAttributeName( self: Self ) -> str: - """Get the name of the attribute to create. - - Returns: - str: name of the attribute. - """ - return self.m_attributeName - - def SetAttributeName( self: Self, name: str ) -> None: - """Set the name of the attribute to create. - - Args: - name (str): name of the attribute. - """ - self.m_attributeName = name - - def GetCreateAttribute( self: Self ) -> bool: - """Get the value of the boolean to create the attribute. - - Returns: - bool: create attribute boolean value. - """ - return self.m_createAttribute - - def SetCreateAttribute( self: Self, value: bool ) -> None: - """Set the value of the boolean to create the attribute. - - Returns: - bool: True to create the attribute, False otherwise. - """ - self.m_createAttribute = value - - def GetSurfaceToVolumeConnectionSets( self: Self ) -> ConnectionSetCollection: - """Get the collection of surface to volume cell ids. - - Returns: - ConnectionSetCollection: collection of ConnectionSet. - """ - return self.m_connectionSets - - def GetVolumeToSurfaceConnectionSets( self: Self ) -> ConnectionSetCollection: - """Get the ConnectionSetCollection of volume to surface cell ids. - - Returns: - ConnectionSetCollection: reversed collection of connection set. - """ - return self.m_connectionSets.getReversedConnectionSetCollection() - - def RequestDataObject( - self: Self, - request: vtkInformation, - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestDataObject. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - inData = self.GetInputData( inInfoVec, 0, 0 ) - outData = self.GetOutputData( outInfoVec, 0 ) - assert inData is not None - if outData is None or ( not outData.IsA( inData.GetClassName() ) ): - outData = inData.NewInstance() - outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) - return super().RequestDataObject( request, inInfoVec, outInfoVec ) - - 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(), "vtkUnstructuredGrid" ) - else: - 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, - 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 volume mesh - self.m_volumeMesh = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) - # input surface - self.m_surfaceMesh = vtkPolyData.GetData( inInfoVec[ 1 ] ) - # output volume mesh - outData: vtkUnstructuredGrid = self.GetOutputData( outInfoVec, 0 ) - outData.ShallowCopy( self.m_volumeMesh ) - - # volume mesh properties - ptsVol: vtkPoints = self.m_volumeMesh.GetPoints() - assert ptsVol is not None, "Volume mesh points are undefined." - ptsCoordsVol: npt.NDArray[ np.float64 ] = vnp.vtk_to_numpy( ptsVol.GetData() ) - - # surface mesh properties - nbFaces: int = self.m_surfaceMesh.GetNumberOfCells() - ptsSurf: vtkPoints = self.m_surfaceMesh.GetPoints() - assert ptsSurf is not None, "Surface mesh points are undefined." - ptsCoordsSurf: npt.NDArray[ np.float64 ] = vnp.vtk_to_numpy( ptsSurf.GetData() ) - - # get cell ids for each face id - self.m_connectionSets.clear() - for faceId in range( nbFaces ): - # self.m_cellIdMap[faceId] = self.getAdjacentCells( - # faceId, ptsCoordsSurf, ptsCoordsVol) - cellIdsSide: dict[ int, bool ] = self.getAdjacentCells( faceId, ptsCoordsSurf, ptsCoordsVol ) - self.m_connectionSets.add( ConnectionSet( faceId, cellIdsSide ) ) - - if self.m_createAttribute: - self.createAttribute( outData ) - except AssertionError as e: - mess: str = "Surface to Volume mesh mapping failed due to:" - self.m_logger.error( mess ) - self.m_logger.error( e, exc_info=True ) - return 0 - except Exception as e: - mess0: str = "Surface to Volume mesh mapping failed due to:" - self.m_logger.critical( mess0 ) - self.m_logger.critical( e, exc_info=True ) - return 0 - return 1 - - def getAdjacentCells( - self: Self, - faceId: int, - ptsCoordsSurf: npt.NDArray[ np.float64 ], - ptsCoordsVol: npt.NDArray[ np.float64 ], - ) -> dict[ int, bool ]: - """Get the cells from volume mesh adjacent to the face cellIdSurf. - - Args: - faceId (int): id of the face of the surface mesh. - ptsCoordsSurf (npt.NDArray[np.float64]): coordinates of the points - of the surface mesh - ptsCoordsVol (npt.NDArray[np.float64]): coordinates of the points - of the volume mesh - - Returns: - dict[int, bool]: map of cell ids adjacent to the face as keys, and - side as values - """ - # Retrieve point ids of the face of the surface - facePtIds: vtkIdList = vtkIdList() - self.m_surfaceMesh.GetCellPoints( faceId, facePtIds ) - - # number of face vertices - nbPtsFace: int = facePtIds.GetNumberOfIds() - # coordinates of the vertices of the face - ptsCoordsFace: list[ npt.NDArray[ np.float64 ] ] = [ - ptsCoordsSurf[ facePtIds.GetId( p ) ] for p in range( nbPtsFace ) - ] - - # get the ids of all the cells that are adjacent to the face - cellIds: set[ int ] = self.getCellIds( ptsCoordsFace, ptsCoordsVol ) - - # get the side of each cell - cellIdsSide: dict[ int, bool ] = {} - for cellId in cellIds: - side: bool = self.getCellSide( cellId, faceId, ptsCoordsFace, ptsCoordsVol ) - cellIdsSide[ cellId ] = side - - return cellIdsSide - - def getCellIds( - self: Self, - ptsCoordsFace: list[ npt.NDArray[ np.float64 ] ], - ptsCoordsVol: npt.NDArray[ np.float64 ], - ) -> set[ int ]: - """Get the ids of all the cells that are adjacent to input face. - - A cell is adjacent to a face if it contains all the vertices of the face. - - .. WARNING:: - Face adjacency determination relies on a geometrical test of vertex - location, it is error prone because of computational precision. - - Args: - ptsCoordsFace (npt.NDArray[np.float64]): List of vertex coordinates - of the face. The shape of the array should be (N,1). - ptsCoordsVol (npt.NDArray[np.float64]): Coodinates of the vertices - of the volume mesh. The shape of the array should be (N,M), where M - is the number of vertices in the volume mesh. - - Returns: - set[int]: set of cell ids adjacent to the face. - """ - cellIds: list[ set[ int ] ] = [] - # get the ids of the cells that contains each vertex of the face. - for coords in ptsCoordsFace: - assert coords.shape[ 0 ] == ptsCoordsVol.shape[ 1 ] - diffCoords: npt.NDArray[ np.float64 ] = np.abs( ptsCoordsVol - coords ) - idPtsVol: tuple[ npt.NDArray[ np.int64 ], ...] = np.where( np.all( diffCoords < EPSILON, axis=-1 ) ) - - assert len( idPtsVol ) == 1 - idPt = idPtsVol[ 0 ] - - # In case of collocated points from the volume mesh - if len( idPt ) > 1: - self.m_logger.warning( "The volume mesh has collocated points: " + f"{idPt}. Clean the mesh first." ) - - # Retrieve all the cells attached to the point of the volume mesh - cellList: vtkIdList = vtkIdList() - self.m_volumeMesh.GetPointCells( idPt[ 0 ], cellList ) - cellIdsVertex: set[ int ] = { cellList.GetId( c ) for c in range( cellList.GetNumberOfIds() ) } - cellIds += [ cellIdsVertex ] - # keep the cells that contain all the vertices of the face - return set.intersection( *cellIds ) - - def getCellSide( - self: Self, - cellId: int, - faceId: int, - ptsCoordsFace: list[ npt.NDArray[ np.float64 ] ], - ptsCoordsVol: npt.NDArray[ np.float64 ], - ) -> bool: - """Get the side of the cell from volume mesh against the surface. - - Args: - cellId (int): cell id in the volume mesh. - faceId (int): cell id in the surface mesh. - ptsCoordsFace (list[npt.NDArray[np.float64]]): coordinates of the - vertices of the face. - ptsCoordsVol (npt.NDArray[np.float64]): coordinates of the vertices - of the cell. - - Returns: - bool: True if the cell is the same side as surface normal, False - otherwise. - """ - # Retrieve vertex coordinates of the cell of the volume mesh - cellPtIds: vtkIdList = vtkIdList() - self.m_volumeMesh.GetCellPoints( cellId, cellPtIds ) - cellPtsCoords: npt.NDArray[ np.float64 ] = np.array( - [ ptsCoordsVol[ cellPtIds.GetId( i ) ] for i in range( cellPtIds.GetNumberOfIds() ) ] ) - # get face normal vector # type: ignore[no-untyped-call] - normals: npt.NDArray[ np.float64 ] = vnp.vtk_to_numpy( self.m_surfaceMesh.GetCellData().GetNormals() ) - normalVec: npt.NDArray[ np.float64 ] = normals[ faceId ] - normalVec /= np.linalg.norm( normalVec ) # normalization - # get cell side - return getCellSideAgainstPlane( cellPtsCoords, ptsCoordsFace[ 0 ], normalVec ) - - def createAttribute( self: Self, mesh: vtkUnstructuredGrid ) -> bool: - """Create the cell adjacency attribute. - - Attribute yields -1 is a cell is not adjacent to input surface, 0 if - the cell is along negative side, 1 if the cell is along positive side. - - Args: - mesh (vtkUnstructuredGrid): mesh where to add the attribute - - Returns: - bool: True if the new attribute was successfully added, - False otherwise - """ - adjacentCellIdsMap: ConnectionSetCollection = ( self.GetVolumeToSurfaceConnectionSets() ) - array: npt.NDArray[ np.float64 ] = -1 * np.ones( mesh.GetNumberOfCells() ) - for i in range( mesh.GetNumberOfCells() ): - connectionSet: Optional[ ConnectionSet ] = adjacentCellIdsMap.get( i ) - # cell i is not in the collection - if connectionSet is None: - continue - # get connected cell sides - values: tuple[ bool, ...] = tuple( connectionSet.getConnectedCellIds().values() ) - # assume that all boolean values are the same, keep the first one - # convert boolean to int - array[ i ] = int( values[ 0 ] ) - newAttr: vtkDataArray = vnp.numpy_to_vtk( array, deep=True, array_type=VTK_INT ) - newAttr.SetName( self.m_attributeName ) - mesh.GetCellData().AddArray( newAttr ) - mesh.GetCellData().Modified() - mesh.Modified() - return True diff --git a/geos-pv/src/geos/pv/plugins/PVAttributeMapping.py b/geos-pv/src/geos/pv/plugins/PVAttributeMapping.py new file mode 100644 index 000000000..243ac7f9f --- /dev/null +++ b/geos-pv/src/geos/pv/plugins/PVAttributeMapping.py @@ -0,0 +1,199 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Raphaël Vinour, Martin Lemay, Romain Baville +# ruff: noqa: E402 # disable Module level import not at top of file +import sys +from pathlib import Path +from typing import Union +from typing_extensions import Self + +# 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.mesh.processing.AttributeMapping import AttributeMapping +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 vtkmodules.vtkCommonCore import ( + vtkInformation, + vtkInformationVector, +) +from vtkmodules.vtkCommonDataModel import ( + vtkCompositeDataSet, + vtkDataSet, + vtkMultiBlockDataSet, +) + +__doc__ = """ +AttributeMapping is a paraview plugin that transfers global attributes from a source mesh to a final mesh with same point/cell coordinates. + +Input and output meshes can be vtkDataSet or vtkMultiBlockDataSet. + +.. Warning:: + For one application of the plugin, the attributes to transfer should all be located on the same piece (all on points or all on cells). + +.. Note:: + For cell, the coordinates of the points in the cell are compared. + If one of the two meshes is a surface and the other a volume, all the points of the surface must be points of the volume. + +To use it: + +* Load the module in Paraview: Tools>Manage Plugins...>Load new>PVAttributeMapping. +* Select the mesh to transfer the global attributes (meshTo). +* Select Filters > 4- Geos Utils > Attribute Mapping. +* Select the source mesh with global attributes to transfer (meshFrom). +* Select the on witch element (onPoints/onCells) the attributes to transfer are. +* Select the global attributes to transfer from the source mesh to the final mesh. +* Apply. + +""" + + +@smproxy.filter( name="PVAttributeMapping", label="Attribute Mapping" ) +@smhint.xml( '' ) +@smproperty.input( name="meshFrom", port_index=1, label="Mesh From" ) +@smdomain.datatype( + dataTypes=[ "vtkDataSet", "vtkMultiBlockDataSet" ], + composite_data_supported=True, +) +@smproperty.input( name="meshTo", port_index=0, label="Mesh To" ) +@smdomain.datatype( + dataTypes=[ "vtkDataSet", "vtkMultiBlockDataSet" ], + composite_data_supported=True, +) +class PVAttributeMapping( VTKPythonAlgorithmBase ): + + def __init__( self: Self ) -> None: + """Map attributes of the source mesh (meshFrom) to the final mesh (meshTo).""" + super().__init__( nInputPorts=2, nOutputPorts=1, inputType="vtkObject", outputType="vtkObject" ) + + self.onPoints: bool = False + self.clearAttributeNames = True + self.attributeNames: list[ str ] = [] + + @smproperty.intvector( + name="AttributePiece", + default_values=1, + number_of_elements=1, + ) + @smdomain.xml( """ + + + + + + """ ) + def setAttributePiece( self: Self, piece: int ) -> None: + """Set attributes piece (points or cells). + + Args: + piece (int): 0 if on points, 1 if on cells. + """ + self.onPoints = bool( piece ) + self.Modified() + + @smproperty.stringvector( + name="SelectAttributeToTransfer", + label="Select Attribute To Transfer", + repeat_command=1, + number_of_elements_per_command="1", + element_types="2", + default_values="None", + ) + @smdomain.xml( """ + + + + + + + + Select attributes to transfer from the source mesh to the final mesh. + + + + + """ ) + def selectMultipleAttribute( self: Self, name: str ) -> None: + """Set the attribute to transfer from the source mesh to the final mesh. + + Args: + name (str): The name of the attribute to transfer. + """ + if self.clearAttributeNames: + self.attributeNames = [] + self.clearAttributeNames = False + + if name != "None": + self.attributeNames.append( name ) + self.Modified() + + def RequestDataObject( + self: Self, + request: vtkInformation, + inInfoVec: list[ vtkInformationVector ], + outInfoVec: vtkInformationVector, + ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestDataObject. + + Args: + request (vtkInformation): Request. + inInfoVec (list[vtkInformationVector]): Input objects. + outInfoVec (vtkInformationVector): Output objects. + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + inDataTo = self.GetInputData( inInfoVec, 0, 0 ) + outData = self.GetOutputData( outInfoVec, 0 ) + assert inDataTo is not None + if outData is None or ( not outData.IsA( inDataTo.GetClassName() ) ): + outData = inDataTo.NewInstance() + outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) + return super().RequestDataObject( request, inInfoVec, outInfoVec ) # type: ignore[no-any-return] + + 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. + """ + meshTo: Union[ vtkDataSet, vtkMultiBlockDataSet, vtkCompositeDataSet ] = self.GetInputData( inInfoVec, 0, 0 ) + meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet, vtkCompositeDataSet ] = self.GetInputData( inInfoVec, 1, 0 ) + outData: Union[ vtkDataSet, vtkMultiBlockDataSet, vtkCompositeDataSet ] = self.GetOutputData( outInfoVec, 0 ) + + assert meshTo is not None, "The final mesh (meshTo) where to transfer attributes is null." + assert meshFrom is not None, "The source mesh (meshFrom) with attributes to transfer is null." + assert outData is not None, "Output pipeline is null." + + outData.ShallowCopy( meshTo ) + + filter: AttributeMapping = AttributeMapping( meshFrom, outData, set( self.attributeNames ), self.onPoints, + True ) + if not filter.logger.hasHandlers(): + filter.setLoggerHandler( VTKHandler() ) + + filter.applyFilter() + self.clearAttributeNames = True + + return 1