refactor: Move GeosBlockExtractor from geos-posp to geos-proccessing

.github/workflows/python-package.yml

@@ -396,4 +396,3 @@ jobs:
           echo ""
           echo "✓ CI requirements satisfied - PR can be merged"
         fi
-

docs/geos_mesh_docs/processing.rst

@@ -27,6 +27,15 @@ geos.mesh.processing.FillPartialArrays filter
     :undoc-members:
     :show-inheritance:
 
+geos.mesh.processing.GeosBlockExtractor module
+-----------------------------------------------
+
+.. automodule:: geos.mesh.processing.GeosBlockExtractor
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
 geos.mesh.processing.meshQualityMetricHelpers module
 -----------------------------------------------------
 

docs/geos_posp_docs/filters.rst

@@ -3,14 +3,6 @@ vtk Filters
 
 This package defines vtk filters that allows to process Geos outputs.
 
-geos_posp.filters.GeosBlockExtractor module
------------------------------------------------
-
-.. automodule:: geos_posp.filters.GeosBlockExtractor
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
 geos_posp.filters.GeosBlockMerge module
 -------------------------------------------
 

geos-mesh/src/geos/mesh/processing/GeosBlockExtractor.py

@@ -0,0 +1,215 @@
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Martin Lemay, Romain Baville
+import logging
+from dataclasses import dataclass
+from typing_extensions import Self
+
+from geos.utils.Logger import ( Logger, getLogger )
+from geos.utils.GeosOutputsConstants import ( GeosDomainNameEnum )
+from geos.mesh.utils.arrayHelpers import ( getCellDimension )
+from geos.mesh.utils.multiblockHelpers import ( getBlockIndexFromName )
+
+from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet
+from vtkmodules.vtkFiltersExtraction import vtkExtractBlock
+
+__doc__ = """
+GeosBlockExtractor is a vtk filter that allows to extract the domain (volume, fault and well) from a GEOS output multiBlockDataset mesh.
+
+.. Important::
+    The input mesh must be an output of a GEOS simulation or contain at least three blocks labeled with the same domain names:
+        "CellElementRegion" for volume domain
+        "SurfaceElementRegion" for fault domain
+        "WellElementRegion" for well domain
+        See more https://geosx-geosx.readthedocs-hosted.com/en/latest/docs/sphinx/datastructure/ElementRegions.html?_sm_au_=iVVT5rrr5fN00R8sQ0WpHK6H8sjL6#xml-element-elementregions
+
+.. Note::
+    Volume domain is automatically extracted, by defaults Fault and Well domains are empty multiBlockDataSet.
+
+To use the filter:
+
+.. code-block:: python
+
+    from geos.mesh.processing.GeosBlockExtractor import GeosBlockExtractor
+
+    # Filter inputs.
+    geosMesh: vtkMultiBlockDataSet
+    # Optional inputs.
+    extractFault: bool # Defaults to False
+    extractWell: bool # Defaults to False
+    speHandler: bool # Defaults to False
+
+    # Instantiate the filter
+    filter: GeosBlockExtractor = GeosBlockExtractor( geosMesh, extractFault, extractWell, speHandler )
+
+    # Set the handler of yours (only if speHandler is True).
+    yourHandler: logging.Handler
+    filter.setLoggerHandler( yourHandler )
+
+    # Do calculations
+    filter.applyFilter()
+
+    # Get the multiBlockDataSet with blocks of the extracted domain.
+    geosDomainExtracted: vtkMultiBlockDataSet
+    geosDomainExtracted = filter.extractedGeosDomain.volume # For volume domain
+    geosDomainExtracted = filter.extractedGeosDomain.fault # For fault domain
+    geosDomainExtracted = filter.extractedGeosDomain.well # For well domain
+"""
+
+loggerTitle: str = "Geos Block Extractor Filter"
+
+
+class GeosExtractDomainBlock( vtkExtractBlock ):
+
+    def __init__( self: Self ) -> None:
+        """Extract blocks from a GEOS output multiBlockDataset mesh."""
+
+    def AddGeosDomainName( self: Self, geosDomainName: GeosDomainNameEnum ) -> None:
+        """Add the index of the GEOS domain to extract from its name.
+
+        Args:
+            geosDomainName (GeosDomainNameEnum): Name of the GEOS domain to extract.
+        """
+        domainBlockIndex: int = getBlockIndexFromName( self.GetInput(), geosDomainName.value )
+        return super().AddIndex( domainBlockIndex )
+
+
+class GeosBlockExtractor:
+
+    @dataclass
+    class ExtractedGeosDomain:
+        """The dataclass with the three GEOS domain mesh."""
+        _volume: vtkMultiBlockDataSet = vtkMultiBlockDataSet()
+        _fault: vtkMultiBlockDataSet = vtkMultiBlockDataSet()
+        _well: vtkMultiBlockDataSet = vtkMultiBlockDataSet()
+
+        @property
+        def volume( self: Self ) -> vtkMultiBlockDataSet:
+            """Get the mesh with the blocks of the GEOS CellElementRegion."""
+            return self._volume
+
+        @volume.setter
+        def volume( self: Self, multiBlockDataSet: vtkMultiBlockDataSet ) -> None:
+            cellDim: set[ int ] = getCellDimension( multiBlockDataSet )
+            if len( cellDim ) == 1 and 3 in cellDim:
+                self._volume.DeepCopy( multiBlockDataSet )
+            else:
+                raise TypeError( "The input mesh must be a volume mesh with cells dimension equal to 3." )
+
+        @property
+        def fault( self: Self ) -> vtkMultiBlockDataSet:
+            """Get the mesh with the blocks of the GEOS SurfaceElementRegion."""
+            return self._fault
+
+        @fault.setter
+        def fault( self: Self, multiBlockDataSet: vtkMultiBlockDataSet ) -> None:
+            cellDim: set[ int ] = getCellDimension( multiBlockDataSet )
+            if len( cellDim ) == 1 and 2 in cellDim:
+                self._fault.DeepCopy( multiBlockDataSet )
+            else:
+                raise TypeError( "The input mesh must be a surface mesh with cells dimension equal to 2." )
+
+        @property
+        def well( self: Self ) -> vtkMultiBlockDataSet:
+            """Get the mesh with the blocks of the GEOS WellElementRegion."""
+            return self._well
+
+        @well.setter
+        def well( self: Self, multiBlockDataSet: vtkMultiBlockDataSet ) -> None:
+            cellDim: set[ int ] = getCellDimension( multiBlockDataSet )
+            if len( cellDim ) == 1 and 1 in cellDim:
+                self._well.DeepCopy( multiBlockDataSet )
+            else:
+                raise TypeError( "The input mesh must be a segment mesh with cells dimension equal to 1." )
+
+        def setExtractedDomain( self: Self, geosDomainName: GeosDomainNameEnum,
+                                multiBlockDataSet: vtkMultiBlockDataSet ) -> None:
+            """Set the mesh to the correct domain.
+
+            Args:
+                geosDomainName (GeosDomainNameEnum): Name of the GEOS domain.
+                multiBlockDataSet (vtkMultiBlockDataSet): The mesh to set.
+            """
+            if geosDomainName.value == "CellElementRegion":
+                self.volume = multiBlockDataSet
+            elif geosDomainName.value == "SurfaceElementRegion":
+                self.fault = multiBlockDataSet
+            elif geosDomainName.value == "WellElementRegion":
+                self.well = multiBlockDataSet
+            else:
+                raise ValueError(
+                    f"The GEOS extractable domains are { GeosDomainNameEnum.VOLUME_DOMAIN_NAME.value }, { GeosDomainNameEnum.FAULT_DOMAIN_NAME.value } and { GeosDomainNameEnum.WELL_DOMAIN_NAME.value }."
+                )
+
+    extractedGeosDomain: ExtractedGeosDomain
+
+    def __init__(
+        self: Self,
+        geosMesh: vtkMultiBlockDataSet,
+        extractFault: bool = False,
+        extractWell: bool = False,
+        speHandler: bool = False,
+    ) -> None:
+        """Blocks from the ElementRegions from a GEOS output multiBlockDataset mesh.
+
+        Args:
+            geosMesh (vtkMultiBlockDataSet): The mesh from Geos.
+            extractFault (bool, Optional): True if SurfaceElementRegion needs to be extracted, False otherwise.
+                Defaults to False.
+            extractWell (bool, Optional): True if WellElementRegion needs to be extracted, False otherwise.
+                Defaults to False.
+            speHandler (bool, optional): True to use a specific handler, False to use the internal handler.
+                Defaults to False.
+        """
+        self.geosMesh: vtkMultiBlockDataSet = geosMesh
+        self.extractedGeosDomain = self.ExtractedGeosDomain()
+
+        self.domainToExtract: list[ GeosDomainNameEnum ] = [ GeosDomainNameEnum.VOLUME_DOMAIN_NAME ]
+        if extractFault:
+            self.domainToExtract.append( GeosDomainNameEnum.FAULT_DOMAIN_NAME )
+        if extractWell:
+            self.domainToExtract.append( GeosDomainNameEnum.WELL_DOMAIN_NAME )
+
+        # 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 applyFilter( self: Self ) -> None:
+        """Extract the volume, the fault or the well domain of the mesh from GEOS."""
+        self.logger.info( f"Apply filter { self.logger.name }." )
+
+        try:
+            extractGeosDomain: GeosExtractDomainBlock = GeosExtractDomainBlock()
+            extractGeosDomain.SetInputData( self.geosMesh )
+
+            for domain in self.domainToExtract:
+                extractGeosDomain.RemoveAllIndices()
+                extractGeosDomain.AddGeosDomainName( domain )
+                extractGeosDomain.Update()
+                self.extractedGeosDomain.setExtractedDomain( domain, extractGeosDomain.GetOutput() )
+
+            self.logger.info( "The filter succeeded." )
+
+        except ValueError as ve:
+            self.logger.error( f"The filter failed.\n{ ve }." )
+        except TypeError as te:
+            self.logger.error( f"The filter failed.\n{ te }." )

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

@@ -26,6 +26,58 @@
 """
 
 
+def getCellDimension( mesh: Union[ vtkMultiBlockDataSet, vtkDataSet ] ) -> set[ int ]:
+    """Get the set of the different cells dimension of a mesh.
+
+    Args:
+        mesh (Union[vtkMultiBlockDataSet, vtkDataSet]): The input mesh with the cells dimension to get.
+
+    Returns:
+        set[int]: The set of the different cells dimension in the input mesh.
+    """
+    if isinstance( mesh, vtkDataSet ):
+        return getCellDimensionDataSet( mesh )
+    elif isinstance( mesh, vtkMultiBlockDataSet ):
+        return getCellDimensionMultiBlockDataSet( mesh )
+    else:
+        raise TypeError( "The input mesh must be a vtkMultiBlockDataSet or a vtkDataSet." )
+
+
+def getCellDimensionMultiBlockDataSet( multiBlockDataSet: vtkMultiBlockDataSet ) -> set[ int ]:
+    """Get the set of the different cells dimension of a multiBlockDataSet.
+
+    Args:
+        multiBlockDataSet (vtkMultiBlockDataSet): The input mesh with the cells dimension to get.
+
+    Returns:
+        set[int]: The set of the different cells dimension in the input multiBlockDataSet.
+    """
+    cellDim: set[ int ] = set()
+    listFlatIdDataSet: list[ int ] = getBlockElementIndexesFlatten( multiBlockDataSet )
+    for flatIdDataSet in listFlatIdDataSet:
+        dataSet: vtkDataSet = vtkDataSet.SafeDownCast( multiBlockDataSet.GetDataSet( flatIdDataSet ) )
+        cellDim = cellDim.union( getCellDimensionDataSet( dataSet ) )
+    return cellDim
+
+
+def getCellDimensionDataSet( dataSet: vtkDataSet ) -> set[ int ]:
+    """Get the set of the different cells dimension of a dataSet.
+
+    Args:
+        dataSet (vtkDataSet): The input mesh with the cells dimension to get.
+
+    Returns:
+        set[int]: The set of the different cells dimension in the input dataSet.
+    """
+    cellDim: set[ int ] = set()
+    cellIter = dataSet.NewCellIterator()
+    cellIter.InitTraversal()
+    while not cellIter.IsDoneWithTraversal():
+        cellDim.add( cellIter.GetCellDimension() )
+        cellIter.GoToNextCell()
+    return cellDim
+
+
 def computeElementMapping(
     meshFrom: Union[ vtkDataSet, vtkMultiBlockDataSet ],
     meshTo: Union[ vtkDataSet, vtkMultiBlockDataSet ],

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

@@ -190,29 +190,6 @@ def getBlockFromFlatIndex( multiBlockDataSet: Union[ vtkMultiBlockDataSet, vtkCo
     return None
 
 
-def getBlockFromName( multiBlockDataSet: Union[ vtkMultiBlockDataSet, vtkCompositeDataSet ],
-                      blockName: str ) -> Union[ None, vtkDataObject ]:
-    """Get the block named blockName from the vtkMultiBlockDataSet.
-
-    Args:
-        multiBlockDataSet (vtkMultiBlockDataSet | vtkCompositeDataSet): MultiBlockDataSet with the block to get.
-        blockName (str): The name of the block to get.
-
-    Returns:
-        Union[None, vtkDataObject]: The block name blockName if it exists, None otherwise
-    """
-    # initialize data object tree iterator
-    iterator: vtkDataObjectTreeIterator = vtkDataObjectTreeIterator()
-    iterator.SetDataSet( multiBlockDataSet )
-    iterator.VisitOnlyLeavesOff()
-    iterator.GoToFirstItem()
-    while iterator.GetCurrentDataObject() is not None:
-        if iterator.GetCurrentMetaData().Get( vtkMultiBlockDataSet.NAME() ) == blockName:
-            return iterator.GetCurrentDataObject()
-        iterator.GoToNextItem()
-    return None
-
-
 def extractBlock( multiBlockDataSet: vtkMultiBlockDataSet, blockIndex: int ) -> vtkMultiBlockDataSet:
     """Extract the block with index blockIndex from multiBlockDataSet.
 

geos-mesh/tests/conftest.py

@@ -10,7 +10,7 @@
 import numpy.typing as npt
 
 from vtkmodules.vtkCommonDataModel import vtkDataSet, vtkMultiBlockDataSet, vtkPolyData
-from vtkmodules.vtkIOXML import vtkXMLGenericDataObjectReader, vtkXMLMultiBlockDataReader
+from vtkmodules.vtkIOXML import vtkXMLGenericDataObjectReader
 
 
 @pytest.fixture
@@ -178,6 +178,10 @@ def _get_dataset( datasetType: str ) -> Union[ vtkMultiBlockDataSet, vtkPolyData
             vtkFilename = "data/fracture_res5_id.vtp"
         elif datasetType == "emptypolydata":
             vtkFilename = "data/fracture_res5_id_empty.vtp"
+        elif datasetType == "meshGeosExtractBlockTmp":
+            vtkFilename = "data/meshGeosExtractBlockTmp.vtm"
+        elif datasetType == "well":
+            vtkFilename = "data/well.vtu"
         datapath: str = os.path.join( os.path.dirname( os.path.realpath( __file__ ) ), vtkFilename )
         reader.SetFileName( datapath )
         reader.Update()

geos-mesh/tests/data/meshGeosExtractBlockTmp.vtm

@@ -0,0 +1,16 @@
+<?xml version="1.0"?>
+<VTKFile type="vtkMultiBlockDataSet" version="1.0">
+    <vtkMultiBlockDataSet>
+        <Block name="CellElementRegion">
+            <DataSet name="domain" file="domain_res5_id.vtu"/>
+            <DataSet name="emptyDomain" file="domain_res5_id_empty.vtu"/>
+        </Block>
+        <Block name="SurfaceElementRegion">
+            <DataSet name="fracture" file="fracture_res5_id.vtu"/>
+            <DataSet name="emptyFracture" file="fracture_res5_id_empty.vtu"/>
+        </Block>
+        <Block name="WellElementRegion">
+            <DataSet name="well" file="well.vtu"/>
+        </Block>
+    </vtkMultiBlockDataSet>
+</VTKFile>