feat: Add meshing tools and stats

geos-mesh/pyproject.toml

@@ -30,6 +30,7 @@ dependencies = [
     "tqdm >= 4.67",
     "numpy >= 2.2",
     "meshio >= 5.3",
+    "typing_extensions >= 4.12",
 ]
 
 [project.scripts]
@@ -47,7 +48,8 @@ build = [
   "build ~= 1.2"
 ]
 dev = [
-   "mypy", 
+   "mypy",
+   "ruff", 
    "yapf",
 ]
 test = [
@@ -56,9 +58,12 @@ test = [
 ]
 
 [tool.pytest.ini_options]
-addopts = [
-    "--import-mode=importlib",
-]
-pythonpath = [
-    "src",
-]
+addopts = "--import-mode=importlib"
+console_output_style = "count"
+pythonpath = ["src"]
+python_classes = "Test"
+python_files = "test*.py"
+python_functions = "test*"
+testpaths = ["tests"]
+norecursedirs = "bin"
+filterwarnings = []

geos-mesh/src/geos/mesh/doctor/checks/generate_fractures.py

@@ -13,7 +13,7 @@
 from vtkmodules.util.numpy_support import numpy_to_vtk, vtk_to_numpy
 from vtkmodules.util.vtkConstants import VTK_ID_TYPE
 from geos.mesh.doctor.checks.vtk_polyhedron import FaceStream
-from geos.mesh.vtk.helpers import has_invalid_field, to_vtk_id_list, vtk_iter
+from geos.mesh.vtk.helpers import has_invalid_arrays, to_vtk_id_list, vtk_iter
 from geos.mesh.vtk.io import VtkOutput, read_mesh, write_mesh
 """
 TypeAliases cannot be used with Python 3.9. A simple assignment like described there will be used:
@@ -557,7 +557,7 @@ def check( vtk_input_file: str, options: Options ) -> Result:
     try:
         mesh = read_mesh( vtk_input_file )
         # Mesh cannot contain global ids before splitting.
-        if has_invalid_field( mesh, [ "GLOBAL_IDS_POINTS", "GLOBAL_IDS_CELLS" ] ):
+        if has_invalid_arrays( mesh, [ "GLOBAL_IDS_POINTS", "GLOBAL_IDS_CELLS" ] ):
             err_msg: str = ( "The mesh cannot contain global ids for neither cells nor points. The correct procedure " +
                              " is to split the mesh and then generate global ids for new split meshes." )
             logging.error( err_msg )

geos-mesh/src/geos/mesh/model/CellTypeCounts.py

@@ -0,0 +1,107 @@
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Antoine Mazuyer, Martin Lemay
+import numpy as np
+import numpy.typing as npt
+from typing_extensions import Self
+from vtkmodules.vtkCommonDataModel import ( vtkCellTypes, VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_VERTEX, VTK_POLYHEDRON,
+                                            VTK_POLYGON, VTK_PYRAMID, VTK_HEXAHEDRON, VTK_WEDGE,
+                                            VTK_NUMBER_OF_CELL_TYPES )
+
+__doc__ = """
+CellTypeCounts stores the number of elements of each type.
+"""
+
+
+class CellTypeCounts():
+
+    def __init__( self: Self ) -> None:
+        """CellTypeCounts stores the number of cells of each type."""
+        self._counts: npt.NDArray[ np.int64 ] = np.zeros( VTK_NUMBER_OF_CELL_TYPES, dtype=float )
+
+    def __str__( self: Self ) -> str:
+        """Overload __str__ method.
+
+        Returns:
+            str: counts as string.
+        """
+        return self.print()
+
+    def __add__( self: Self, other: Self ) -> Self:
+        """Addition operator.
+
+        CellTypeCounts addition consists in suming counts.
+
+        Args:
+            other (Self): other CellTypeCounts object
+
+        Returns:
+            Self: new CellTypeCounts object
+        """
+        assert isinstance( other, CellTypeCounts ), "Other object must be a CellTypeCounts."
+        newCounts: CellTypeCounts = CellTypeCounts()
+        newCounts._counts = self._counts + other._counts
+        return newCounts
+
+    def addType( self: Self, cellType: int ) -> None:
+        """Increment the number of cell of input type.
+
+        Args:
+            cellType (int): cell type
+        """
+        self._counts[ cellType ] += 1
+        self._updateGeneralCounts( cellType, 1 )
+
+    def setTypeCount( self: Self, cellType: int, count: int ) -> None:
+        """Set the number of cells of input type.
+
+        Args:
+            cellType (int): cell type
+            count (int): number of cells
+        """
+        prevCount = self._counts[ cellType ]
+        self._counts[ cellType ] = count
+        self._updateGeneralCounts( cellType, count - prevCount )
+
+    def getTypeCount( self: Self, cellType: int ) -> int:
+        """Get the number of cells of input type.
+
+        Args:
+            cellType (int): cell type
+
+        Returns:
+            int: number of cells
+        """
+        return int( self._counts[ cellType ] )
+
+    def _updateGeneralCounts( self: Self, cellType: int, count: int ) -> None:
+        """Update generic type counters.
+
+        Args:
+            cellType (int): cell type
+            count (int): count increment
+        """
+        if ( cellType != VTK_POLYGON ) and ( vtkCellTypes.GetDimension( cellType ) == 2 ):
+            self._counts[ VTK_POLYGON ] += count
+        if ( cellType != VTK_POLYHEDRON ) and ( vtkCellTypes.GetDimension( cellType ) == 3 ):
+            self._counts[ VTK_POLYHEDRON ] += count
+
+    def print( self: Self ) -> str:
+        """Print counts string.
+
+        Returns:
+            str: counts string.
+        """
+        card: str = ""
+        card += "|                                   |              |\n"
+        card += "|               -                   |       -      |\n"
+        card += f"| **Total Number of Vertices**      | {int(self._counts[VTK_VERTEX]):12} |\n"
+        card += f"| **Total Number of Polygon**       | {int(self._counts[VTK_POLYGON]):12} |\n"
+        card += f"| **Total Number of Polyhedron**    | {int(self._counts[VTK_POLYHEDRON]):12} |\n"
+        card += f"| **Total Number of Cells**         | {int(self._counts[VTK_POLYHEDRON]+self._counts[VTK_POLYGON]):12} |\n"
+        card += "|               -                   |       -      |\n"
+        for cellType in ( VTK_TRIANGLE, VTK_QUAD ):
+            card += f"| **Total Number of {vtkCellTypes.GetClassNameFromTypeId(cellType):<13}** | {int(self._counts[cellType]):12} |\n"
+        for cellType in ( VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON ):
+            card += f"| **Total Number of {vtkCellTypes.GetClassNameFromTypeId(cellType):<13}** | {int(self._counts[cellType]):12} |\n"
+        return card

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

@@ -0,0 +1,184 @@
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Antoine Mazuyer, Martin Lemay
+import numpy as np
+from typing_extensions import Self
+from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase
+from vtkmodules.vtkCommonCore import (
+    vtkIntArray,
+    vtkInformation,
+    vtkInformationVector,
+    vtkPoints,
+    reference,
+    vtkIdList,
+)
+from vtkmodules.vtkCommonDataModel import (
+    vtkUnstructuredGrid,
+    vtkIncrementalOctreePointLocator,
+    vtkCellTypes,
+    vtkCell,
+)
+
+__doc__ = """
+MergeColocatedPoints module is a vtk filter that merges colocated points from input mesh.
+
+Filter input and output types are vtkUnstructuredGrid.
+
+.. Warning:: This operation uses geometrical tests that may not be accurate in case of very small cells.
+
+
+To use the filter:
+
+.. code-block:: python
+
+    from geos.mesh.processing.MergeColocatedPoints import MergeColocatedPoints
+
+    # filter inputs
+    input :vtkUnstructuredGrid
+
+    # instanciate the filter
+    filter :MergeColocatedPoints = MergeColocatedPoints()
+    # set input data object
+    filter.SetInputDataObject(input)
+    # do calculations
+    filter.Update()
+    # get output object
+    output :vtkUnstructuredGrid = filter.GetOutputDataObject(0)
+"""
+
+
+class MergeColocatedPoints( VTKPythonAlgorithmBase ):
+
+    def __init__( self: Self ) -> None:
+        """MergeColocatedPoints filter merges duplacted points of the input mesh."""
+        super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkUnstructuredGrid" )
+
+    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" )
+
+    def RequestDataObject(
+        self: Self,
+        request: vtkInformation,  # noqa: F841
+        inInfoVec: list[ vtkInformationVector ],  # noqa: F841
+        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 RequestData(
+        self: Self,
+        request: vtkInformation,  # noqa: F841
+        inInfoVec: list[ vtkInformationVector ],  # noqa: F841
+        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.
+        """
+        inData: vtkUnstructuredGrid = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] )
+        output: vtkUnstructuredGrid = self.GetOutputData( outInfoVec, 0 )
+        assert inData is not None, "Input mesh is undefined."
+        assert output is not None, "Output mesh is undefined."
+        vertexMap: list[ int ] = self.setMergePoints( inData, output )
+        self.setCells( inData, output, vertexMap )
+        return 1
+
+    def setMergePoints( self: Self, input: vtkUnstructuredGrid, output: vtkUnstructuredGrid ) -> list[ int ]:
+        """Merge duplicated points and set new points and attributes to output mesh.
+
+        Args:
+            input (vtkUnstructuredGrid): input mesh
+            output (vtkUnstructuredGrid): output mesh
+
+        Returns:
+            list[int]: list containing new point ids.
+        """
+        vertexMap: list[ int ] = []
+        newPoints: vtkPoints = vtkPoints()
+        # use point locator to check for colocated points
+        pointsLocator = vtkIncrementalOctreePointLocator()
+        pointsLocator.InitPointInsertion( newPoints, input.GetBounds() )
+        # create an array to count the number of colocated points
+        vertexCount: vtkIntArray = vtkIntArray()
+        vertexCount.SetName( "Count" )
+        ptId = reference( 0 )
+        countD: int = 0  # total number of colocated points
+        for v in range( input.GetNumberOfPoints() ):
+            inserted: bool = pointsLocator.InsertUniquePoint( input.GetPoints().GetPoint( v ), ptId )
+            if inserted:
+                vertexCount.InsertNextValue( 1 )
+            else:
+                vertexCount.SetValue( ptId, vertexCount.GetValue( ptId ) + 1 )
+                countD = countD + 1
+            vertexMap += [ ptId.get() ]
+
+        output.SetPoints( pointsLocator.GetLocatorPoints() )
+        # copy point attributes
+        output.GetPointData().DeepCopy( input.GetPointData() )
+        # add the array to points data
+        output.GetPointData().AddArray( vertexCount )
+        return vertexMap
+
+    def setCells( self: Self, input: vtkUnstructuredGrid, output: vtkUnstructuredGrid, vertexMap: list[ int ] ) -> bool:
+        """Set cell point ids and attributes to output mesh.
+
+        Args:
+            input (vtkUnstructuredGrid): input mesh
+            output (vtkUnstructuredGrid): output mesh
+            vertexMap (list[int)]): list containing new point ids
+
+        Returns:
+            bool: True if calculation successfully ended.
+        """
+        nbCells: int = input.GetNumberOfCells()
+        nbPoints: int = output.GetNumberOfPoints()
+        assert np.unique(
+            vertexMap ).size == nbPoints, "The size of the list of point ids must be equal to the number of points."
+        cellTypes: vtkCellTypes = vtkCellTypes()
+        input.GetCellTypes( cellTypes )
+        output.Allocate( nbCells )
+        # create mesh cells
+        for cellId in range( nbCells ):
+            cell: vtkCell = input.GetCell( cellId )
+            # create cells from point ids
+            cellsID: vtkIdList = vtkIdList()
+            for ptId in range( cell.GetNumberOfPoints() ):
+                ptIdOld: int = cell.GetPointId( ptId )
+                ptIdNew: int = vertexMap[ ptIdOld ]
+                cellsID.InsertNextId( ptIdNew )
+            output.InsertNextCell( cell.GetCellType(), cellsID )
+        # copy cell attributes
+        assert output.GetNumberOfCells() == nbCells, "Output and input mesh must have the same number of cells."
+        output.GetCellData().DeepCopy( input.GetCellData() )
+        return True