diff --git a/docs/geos_mesh_docs/model.rst b/docs/geos_mesh_docs/model.rst new file mode 100644 index 00000000..ead57332 --- /dev/null +++ b/docs/geos_mesh_docs/model.rst @@ -0,0 +1,13 @@ +Model +^^^^^^^ + +The `model` module of `geos-mesh` package contains data model. + + +geos.mesh.model.CellTypeCounts filter +-------------------------------------- + +.. automodule:: geos.mesh.model.CellTypeCounts + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/geos_mesh_docs/modules.rst b/docs/geos_mesh_docs/modules.rst index fa6a3558..4e13c711 100644 --- a/docs/geos_mesh_docs/modules.rst +++ b/docs/geos_mesh_docs/modules.rst @@ -11,4 +11,10 @@ GEOS Mesh tools io + model + + processing + + stats + utils \ No newline at end of file diff --git a/docs/geos_mesh_docs/processing.rst b/docs/geos_mesh_docs/processing.rst new file mode 100644 index 00000000..d79db6db --- /dev/null +++ b/docs/geos_mesh_docs/processing.rst @@ -0,0 +1,13 @@ +Processing filters +^^^^^^^^^^^^^^^^^^^ + +The `processing` module of `geos-mesh` package contains filters to process meshes. + + +geos.mesh.processing.SplitMesh filter +-------------------------------------- + +.. automodule:: geos.mesh.processing.SplitMesh + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/geos_mesh_docs/stats.rst b/docs/geos_mesh_docs/stats.rst new file mode 100644 index 00000000..50664514 --- /dev/null +++ b/docs/geos_mesh_docs/stats.rst @@ -0,0 +1,13 @@ +Mesh stats tools +^^^^^^^^^^^^^^^^ + +The `stats` module of `geos-mesh` package contains filter to compute statistics on meshes. + + +geos.mesh.stats.CellTypeCounter filter +-------------------------------------- + +.. automodule:: geos.mesh.stats.CellTypeCounter + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/geos_mesh_docs/utils.rst b/docs/geos_mesh_docs/utils.rst index 62c15c72..31f83c3f 100644 --- a/docs/geos_mesh_docs/utils.rst +++ b/docs/geos_mesh_docs/utils.rst @@ -1,7 +1,7 @@ Mesh utilities ^^^^^^^^^^^^^^^^ -The `utils` module of `geos-mesh` package contains different utilities methods for VTK meshes. +The `utils` module of `geos-mesh` package contains various utilities for VTK meshes. geos.mesh.utils.genericHelpers module diff --git a/geos-mesh/pyproject.toml b/geos-mesh/pyproject.toml index 1a184bf8..8a994b47 100644 --- a/geos-mesh/pyproject.toml +++ b/geos-mesh/pyproject.toml @@ -29,8 +29,9 @@ dependencies = [ "networkx >= 2.4", "tqdm >= 4.67", "numpy >= 2.2", + "pandas >= 2.2", "meshio >= 5.3", - "pandas", + "typing_extensions >= 4.12", ] [project.scripts] @@ -48,7 +49,8 @@ build = [ "build ~= 1.2" ] dev = [ - "mypy", + "mypy", + "ruff", "yapf", ] test = [ @@ -57,7 +59,7 @@ test = [ ] [tool.pytest.ini_options] -addopts="--import-mode=importlib" +addopts = "--import-mode=importlib" console_output_style = "count" pythonpath = ["src"] python_classes = "Test" @@ -66,7 +68,3 @@ python_functions = "test*" testpaths = ["tests"] norecursedirs = "bin" filterwarnings = [] - -[tool.coverage.run] -branch = true -source = ["src/geos"] \ No newline at end of file diff --git a/geos-mesh/src/geos/mesh/doctor/checks/generate_fractures.py b/geos-mesh/src/geos/mesh/doctor/checks/generate_fractures.py index 17198237..82e25b7b 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/generate_fractures.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/generate_fractures.py @@ -13,7 +13,8 @@ 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.utils.arrayHelpers import has_invalid_field +from geos.mesh.utils.arrayHelpers import has_array + from geos.mesh.utils.genericHelpers import to_vtk_id_list, vtk_iter from geos.mesh.io.vtkIO import VtkOutput, read_mesh, write_mesh """ @@ -558,7 +559,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_array( 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 ) diff --git a/geos-mesh/src/geos/mesh/io/vtkIO.py b/geos-mesh/src/geos/mesh/io/vtkIO.py index aa4e4015..1b93648a 100644 --- a/geos-mesh/src/geos/mesh/io/vtkIO.py +++ b/geos-mesh/src/geos/mesh/io/vtkIO.py @@ -160,6 +160,7 @@ def __write_vtu( mesh: vtkUnstructuredGrid, output: str, toBinary: bool = False def write_mesh( mesh: vtkPointSet, vtk_output: VtkOutput, canOverwrite: bool = False ) -> int: """Write mesh to disk. + Nothing is done if file already exists. Args: diff --git a/geos-mesh/src/geos/mesh/model/CellTypeCounts.py b/geos-mesh/src/geos/mesh/model/CellTypeCounts.py new file mode 100644 index 00000000..bab12085 --- /dev/null +++ b/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 ) -> 'CellTypeCounts': + """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 diff --git a/geos-mesh/src/geos/mesh/processing/SplitMesh.py b/geos-mesh/src/geos/mesh/processing/SplitMesh.py new file mode 100644 index 00000000..95a219bc --- /dev/null +++ b/geos-mesh/src/geos/mesh/processing/SplitMesh.py @@ -0,0 +1,475 @@ +# 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.util.vtkAlgorithm import VTKPythonAlgorithmBase +from vtkmodules.vtkCommonCore import ( + vtkPoints, + vtkIdTypeArray, + vtkDataArray, + vtkInformation, + vtkInformationVector, +) +from vtkmodules.vtkCommonDataModel import ( + vtkUnstructuredGrid, + vtkCellArray, + vtkCellData, + vtkCell, + vtkCellTypes, + VTK_TRIANGLE, + VTK_QUAD, + VTK_TETRA, + VTK_HEXAHEDRON, + VTK_PYRAMID, + VTK_WEDGE, + VTK_POLYHEDRON, + VTK_POLYGON, +) + +from vtkmodules.util.numpy_support import ( numpy_to_vtk, vtk_to_numpy ) + +from geos.mesh.stats.CellTypeCounter import CellTypeCounter +from geos.mesh.model.CellTypeCounts import CellTypeCounts + +__doc__ = """ +SplitMesh module is a vtk filter that split cells of a mesh composed of Tetrahedra, pyramids, and hexahedra. + +Filter input and output types are vtkUnstructuredGrid. + +To use the filter: + +.. code-block:: python + + from geos.mesh.processing.SplitMesh import SplitMesh + + # filter inputs + input :vtkUnstructuredGrid + + # instanciate the filter + filter :SplitMesh = SplitMesh() + # set input data object + filter.SetInputDataObject(input) + # do calculations + filter.Update() + # get output object + output :vtkUnstructuredGrid = filter.GetOutputDataObject(0) +""" + + +class SplitMesh( VTKPythonAlgorithmBase ): + + def __init__( self ) -> None: + """SplitMesh filter splits each cell using edge centers.""" + super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkUnstructuredGrid" ) + + self.inData: vtkUnstructuredGrid + self.cells: vtkCellArray + self.points: vtkPoints + self.originalId: vtkIdTypeArray + self.cellTypes: list[ int ] + + 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" ) + return 1 + + 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. + """ + self.inData = self.GetInputData( inInfoVec, 0, 0 ) + output: vtkUnstructuredGrid = self.GetOutputData( outInfoVec, 0 ) + + assert self.inData is not None, "Input mesh is undefined." + assert output is not None, "Output mesh is undefined." + + nb_cells: int = self.inData.GetNumberOfCells() + counts: CellTypeCounts = self._get_cell_counts() + nb_tet: int = counts.getTypeCount( VTK_TETRA ) + nb_pyr: int = counts.getTypeCount( VTK_PYRAMID ) + nb_hex: int = counts.getTypeCount( VTK_HEXAHEDRON ) + nb_triangles: int = counts.getTypeCount( VTK_TRIANGLE ) + nb_quad: int = counts.getTypeCount( VTK_QUAD ) + nb_polygon = counts.getTypeCount( VTK_POLYGON ) + nb_polyhedra = counts.getTypeCount( VTK_POLYHEDRON ) + assert counts.getTypeCount( VTK_WEDGE ) == 0, "Input mesh contains wedges that are not currently supported." + assert nb_polyhedra * nb_polygon == 0, "Input mesh is composed of both polygons and polyhedra, but it must contains only one of the two." + nbNewPoints: int = 0 + nbNewPoints = nb_hex * 19 + nb_tet * 6 + nb_pyr * 9 if nb_polyhedra > 0 else nb_triangles * 3 + nb_quad * 5 + nbNewCells: int = nb_hex * 8 + nb_tet * 8 + nb_pyr * 10 * nb_triangles * 4 + nb_quad * 4 + + self.points = vtkPoints() + self.points.DeepCopy( self.inData.GetPoints() ) + self.points.Resize( self.inData.GetNumberOfPoints() + nbNewPoints ) + + self.cells = vtkCellArray() + self.cells.AllocateExact( nbNewCells, 8 ) + self.originalId = vtkIdTypeArray() + self.originalId.SetName( "OriginalID" ) + self.originalId.Allocate( nbNewCells ) + self.cellTypes = [] + for c in range( nb_cells ): + cell: vtkCell = self.inData.GetCell( c ) + cellType: int = cell.GetCellType() + if cellType == VTK_HEXAHEDRON: + self._split_hexahedron( cell, c ) + elif cellType == VTK_TETRA: + self._split_tetrahedron( cell, c ) + elif cellType == VTK_PYRAMID: + self._split_pyramid( cell, c ) + elif cellType == VTK_TRIANGLE: + self._split_triangle( cell, c ) + elif cellType == VTK_QUAD: + self._split_quad( cell, c ) + else: + raise TypeError( f"Cell type {vtkCellTypes.GetClassNameFromTypeId(cellType)} is not supported." ) + # add points and cells + output.SetPoints( self.points ) + output.SetCells( self.cellTypes, self.cells ) + # add attribute saving original cell ids + cellArrays: vtkCellData = output.GetCellData() + assert cellArrays is not None, "Cell data is undefined." + cellArrays.AddArray( self.originalId ) + # transfer all cell arrays + self._transferCellArrays( output ) + return 1 + + def _get_cell_counts( self: Self ) -> CellTypeCounts: + """Get the number of cells of each type. + + Returns: + CellTypeCounts: cell type counts + """ + filter: CellTypeCounter = CellTypeCounter() + filter.SetInputDataObject( self.inData ) + filter.Update() + return filter.GetCellTypeCounts() + + def _addMidPoint( self: Self, ptA: int, ptB: int ) -> int: + """Add a point at the center of the edge defined by input point ids. + + Args: + ptA (int): first point Id + ptB (int): second point Id + + Returns: + int: inserted point Id + """ + ptACoor: npt.NDArray[ np.float64 ] = np.array( self.points.GetPoint( ptA ) ) + ptBCoor: npt.NDArray[ np.float64 ] = np.array( self.points.GetPoint( ptB ) ) + center: npt.NDArray[ np.float64 ] = ( ptACoor + ptBCoor ) / 2. + return self.points.InsertNextPoint( center[ 0 ], center[ 1 ], center[ 2 ] ) + + def _split_tetrahedron( self: Self, cell: vtkCell, index: int ) -> None: + r"""Split a tetrahedron. + + Let's suppose an input tetrahedron composed of nodes (0, 1, 2, 3), + the cell is splitted in 8 tetrahedra using edge centers. + + 2 + ,/|`\ + ,/ | `\ + ,6 '. `5 + ,/ 8 `\ + ,/ | `\ + 0--------4--'.--------1 + `\. | ,/ + `\. | ,9 + `7. '. ,/ + `\. |/ + `3 + + Args: + cell (vtkCell): cell to split + index (int): index of the cell + """ + pt0: int = cell.GetPointId( 0 ) + pt1: int = cell.GetPointId( 1 ) + pt2: int = cell.GetPointId( 2 ) + pt3: int = cell.GetPointId( 3 ) + pt4: int = self._addMidPoint( pt0, pt1 ) + pt5: int = self._addMidPoint( pt1, pt2 ) + pt6: int = self._addMidPoint( pt0, pt2 ) + pt7: int = self._addMidPoint( pt0, pt3 ) + pt8: int = self._addMidPoint( pt2, pt3 ) + pt9: int = self._addMidPoint( pt1, pt3 ) + + self.cells.InsertNextCell( 4, [ pt0, pt4, pt6, pt7 ] ) + self.cells.InsertNextCell( 4, [ pt7, pt9, pt8, pt3 ] ) + self.cells.InsertNextCell( 4, [ pt9, pt4, pt5, pt1 ] ) + self.cells.InsertNextCell( 4, [ pt5, pt6, pt8, pt2 ] ) + self.cells.InsertNextCell( 4, [ pt6, pt8, pt7, pt4 ] ) + self.cells.InsertNextCell( 4, [ pt4, pt8, pt7, pt9 ] ) + self.cells.InsertNextCell( 4, [ pt4, pt8, pt9, pt5 ] ) + self.cells.InsertNextCell( 4, [ pt5, pt4, pt8, pt6 ] ) + for _ in range( 8 ): + self.originalId.InsertNextValue( index ) + self.cellTypes.extend( [ VTK_TETRA ] * 8 ) + + def _split_pyramid( self: Self, cell: vtkCell, index: int ) -> None: + r"""Split a pyramid. + + Let's suppose an input pyramid composed of nodes (0, 1, 2, 3, 4), + the cell is splitted in 8 pyramids using edge centers. + + 4 + ,/|\ + ,/ .'|\ + ,/ | | \ + ,/ .' | `. + ,7 | 12 \ + ,/ .' | \ + ,/ 9 | 11 + 0--------6-.'----3 `. + `\ | `\ \ + `5 .'13 10 \ + `\ | `\ \ + `\.' `\` + 1--------8-------2 + + Args: + cell (vtkCell): cell to split + index (int): index of the cell + """ + pt0: int = cell.GetPointId( 0 ) + pt1: int = cell.GetPointId( 1 ) + pt2: int = cell.GetPointId( 2 ) + pt3: int = cell.GetPointId( 3 ) + pt4: int = cell.GetPointId( 4 ) + pt5: int = self._addMidPoint( pt0, pt1 ) + pt6: int = self._addMidPoint( pt0, pt3 ) + pt7: int = self._addMidPoint( pt0, pt4 ) + pt8: int = self._addMidPoint( pt1, pt2 ) + pt9: int = self._addMidPoint( pt1, pt4 ) + pt10: int = self._addMidPoint( pt2, pt3 ) + pt11: int = self._addMidPoint( pt2, pt4 ) + pt12: int = self._addMidPoint( pt3, pt4 ) + pt13: int = self._addMidPoint( pt5, pt10 ) + + self.cells.InsertNextCell( 5, [ pt5, pt1, pt8, pt13, pt9 ] ) + self.cells.InsertNextCell( 5, [ pt13, pt8, pt2, pt10, pt11 ] ) + self.cells.InsertNextCell( 5, [ pt3, pt6, pt13, pt10, pt12 ] ) + self.cells.InsertNextCell( 5, [ pt6, pt0, pt5, pt13, pt7 ] ) + self.cells.InsertNextCell( 5, [ pt12, pt7, pt9, pt11, pt4 ] ) + self.cells.InsertNextCell( 5, [ pt11, pt9, pt7, pt12, pt13 ] ) + + self.cells.InsertNextCell( 4, [ pt7, pt9, pt5, pt13 ] ) + self.cells.InsertNextCell( 4, [ pt9, pt11, pt8, pt13 ] ) + self.cells.InsertNextCell( 4, [ pt11, pt12, pt10, pt13 ] ) + self.cells.InsertNextCell( 4, [ pt12, pt7, pt6, pt13 ] ) + for _ in range( 10 ): + self.originalId.InsertNextValue( index ) + self.cellTypes.extend( [ VTK_PYRAMID ] * 8 ) + + def _split_hexahedron( self: Self, cell: vtkCell, index: int ) -> None: + r"""Split a hexahedron. + + Let's suppose an input hexahedron composed of nodes (0, 1, 2, 3, 4, 5, 6, 7), + the cell is splitted in 8 hexahedra using edge centers. + + 3----13----2 + |\ |\ + |15 24 | 14 + 9 \ 20 11 \ + | 7----19+---6 + |22 | 26 | 23| + 0---+-8----1 | + \ 17 25 \ 18 + 10| 21 12| + \| \| + 4----16----5 + + Args: + cell (vtkCell): cell to split + index (int): index of the cell + """ + pt0: int = cell.GetPointId( 0 ) + pt1: int = cell.GetPointId( 1 ) + pt2: int = cell.GetPointId( 2 ) + pt3: int = cell.GetPointId( 3 ) + pt4: int = cell.GetPointId( 4 ) + pt5: int = cell.GetPointId( 5 ) + pt6: int = cell.GetPointId( 6 ) + pt7: int = cell.GetPointId( 7 ) + pt8: int = self._addMidPoint( pt0, pt1 ) + pt9: int = self._addMidPoint( pt0, pt3 ) + pt10: int = self._addMidPoint( pt0, pt4 ) + pt11: int = self._addMidPoint( pt1, pt2 ) + pt12: int = self._addMidPoint( pt1, pt5 ) + pt13: int = self._addMidPoint( pt2, pt3 ) + pt14: int = self._addMidPoint( pt2, pt6 ) + pt15: int = self._addMidPoint( pt3, pt7 ) + pt16: int = self._addMidPoint( pt4, pt5 ) + pt17: int = self._addMidPoint( pt4, pt7 ) + pt18: int = self._addMidPoint( pt5, pt6 ) + pt19: int = self._addMidPoint( pt6, pt7 ) + pt20: int = self._addMidPoint( pt9, pt11 ) + pt21: int = self._addMidPoint( pt10, pt12 ) + pt22: int = self._addMidPoint( pt9, pt17 ) + pt23: int = self._addMidPoint( pt11, pt18 ) + pt24: int = self._addMidPoint( pt14, pt15 ) + pt25: int = self._addMidPoint( pt17, pt18 ) + pt26: int = self._addMidPoint( pt22, pt23 ) + + self.cells.InsertNextCell( 8, [ pt10, pt21, pt26, pt22, pt4, pt16, pt25, pt17 ] ) + self.cells.InsertNextCell( 8, [ pt21, pt12, pt23, pt26, pt16, pt5, pt18, pt25 ] ) + self.cells.InsertNextCell( 8, [ pt0, pt8, pt20, pt9, pt10, pt21, pt26, pt22 ] ) + self.cells.InsertNextCell( 8, [ pt8, pt1, pt11, pt20, pt21, pt12, pt23, pt26 ] ) + self.cells.InsertNextCell( 8, [ pt22, pt26, pt24, pt15, pt17, pt25, pt19, pt7 ] ) + self.cells.InsertNextCell( 8, [ pt26, pt23, pt14, pt24, pt25, pt18, pt6, pt19 ] ) + self.cells.InsertNextCell( 8, [ pt9, pt20, pt13, pt3, pt22, pt26, pt24, pt15 ] ) + self.cells.InsertNextCell( 8, [ pt20, pt11, pt2, pt13, pt26, pt23, pt14, pt24 ] ) + for _ in range( 8 ): + self.originalId.InsertNextValue( index ) + self.cellTypes.extend( [ VTK_HEXAHEDRON ] * 8 ) + + def _split_triangle( self: Self, cell: vtkCell, index: int ) -> None: + r"""Split a triangle. + + Let's suppose an input triangle composed of nodes (0, 1, 2), + the cell is splitted in 3 triangles using edge centers. + + 2 + |\ + | \ + 5 4 + | \ + | \ + 0-----3----1 + + Args: + cell (vtkCell): cell to split + index (int): index of the cell + """ + pt0: int = cell.GetPointId( 0 ) + pt1: int = cell.GetPointId( 1 ) + pt2: int = cell.GetPointId( 2 ) + pt3: int = self._addMidPoint( pt0, pt1 ) + pt4: int = self._addMidPoint( pt1, pt2 ) + pt5: int = self._addMidPoint( pt0, pt2 ) + + self.cells.InsertNextCell( 3, [ pt0, pt3, pt5 ] ) + self.cells.InsertNextCell( 3, [ pt3, pt1, pt4 ] ) + self.cells.InsertNextCell( 3, [ pt5, pt4, pt2 ] ) + self.cells.InsertNextCell( 3, [ pt3, pt4, pt5 ] ) + for _ in range( 4 ): + self.originalId.InsertNextValue( index ) + self.cellTypes.extend( [ VTK_TRIANGLE ] * 4 ) + + def _split_quad( self: Self, cell: vtkCell, index: int ) -> None: + r"""Split a quad. + + Let's suppose an input quad composed of nodes (0, 1, 2, 3), + the cell is splitted in 4 quads using edge centers. + + 3-----6-----2 + | | + | | + 7 8 5 + | | + | | + 0-----4-----1 + + Args: + cell (vtkCell): cell to split + index (int): index of the cell + """ + pt0: int = cell.GetPointId( 0 ) + pt1: int = cell.GetPointId( 1 ) + pt2: int = cell.GetPointId( 2 ) + pt3: int = cell.GetPointId( 3 ) + pt4: int = self._addMidPoint( pt0, pt1 ) + pt5: int = self._addMidPoint( pt1, pt2 ) + pt6: int = self._addMidPoint( pt2, pt3 ) + pt7: int = self._addMidPoint( pt3, pt0 ) + pt8: int = self._addMidPoint( pt7, pt5 ) + + self.cells.InsertNextCell( 4, [ pt0, pt4, pt8, pt7 ] ) + self.cells.InsertNextCell( 4, [ pt4, pt1, pt5, pt8 ] ) + self.cells.InsertNextCell( 4, [ pt8, pt5, pt2, pt6 ] ) + self.cells.InsertNextCell( 4, [ pt7, pt8, pt6, pt3 ] ) + for _ in range( 4 ): + self.originalId.InsertNextValue( index ) + self.cellTypes.extend( [ VTK_QUAD ] * 4 ) + + def _transferCellArrays( self: Self, splittedMesh: vtkUnstructuredGrid ) -> bool: + """Transfer arrays from input mesh to splitted mesh. + + Args: + splittedMesh (vtkUnstructuredGrid): splitted mesh + + Returns: + bool: True if arrays were successfully transfered. + """ + cellDataSplitted: vtkCellData = splittedMesh.GetCellData() + assert cellDataSplitted is not None, "Cell data of splitted mesh should be defined." + cellData: vtkCellData = self.inData.GetCellData() + assert cellData is not None, "Cell data of input mesh should be defined." + # for each array of input mesh + for i in range( cellData.GetNumberOfArrays() ): + array: vtkDataArray = cellData.GetArray( i ) + assert array is not None, "Array should be defined." + npArray: npt.NDArray[ np.float64 ] = vtk_to_numpy( array ) + # get number of components + dims: tuple[ int, ...] = npArray.shape + ny: int = 1 if len( dims ) == 1 else dims[ 1 ] + # create new array with nb cells from splitted mesh and number of components from array to copy + newNpArray: npt.NDArray[ np.float64 ] = np.full( ( splittedMesh.GetNumberOfCells(), ny ), np.nan ) + # for each cell, copy the values from input mesh + for c in range( splittedMesh.GetNumberOfCells() ): + idParent: int = int( self.originalId.GetTuple1( c ) ) + newNpArray[ c ] = npArray[ idParent ] + # set array the splitted mesh + newArray: vtkDataArray = numpy_to_vtk( newNpArray ) + newArray.SetName( array.GetName() ) + cellDataSplitted.AddArray( newArray ) + cellDataSplitted.Modified() + splittedMesh.Modified() + return True diff --git a/geos-mesh/src/geos/mesh/stats/CellTypeCounter.py b/geos-mesh/src/geos/mesh/stats/CellTypeCounter.py new file mode 100644 index 00000000..749432d7 --- /dev/null +++ b/geos-mesh/src/geos/mesh/stats/CellTypeCounter.py @@ -0,0 +1,120 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Antoine Mazuyer, Martin Lemay +from typing_extensions import Self +from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase +from vtkmodules.vtkCommonCore import ( + vtkInformation, + vtkInformationVector, + vtkIntArray, +) +from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkCell, vtkTable, vtkCellTypes, VTK_VERTEX, + VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON ) + +from geos.mesh.model.CellTypeCounts import CellTypeCounts + +__doc__ = """ +CellTypeCounter module is a vtk filter that computes cell type counts. + +Filter input is a vtkUnstructuredGrid, output is a vtkTable + +To use the filter: + +.. code-block:: python + + from geos.mesh.stats.CellTypeCounter import CellTypeCounter + + # filter inputs + input :vtkUnstructuredGrid + + # instanciate the filter + filter :CellTypeCounter = CellTypeCounter() + # set input data object + filter.SetInputDataObject(input) + # do calculations + filter.Update() + # get counts + counts :CellTypeCounts = filter.GetCellTypeCounts() +""" + + +class CellTypeCounter( VTKPythonAlgorithmBase ): + + def __init__( self ) -> None: + """CellTypeCounter filter computes mesh stats.""" + super().__init__( nInputPorts=1, nOutputPorts=1, inputType="vtkUnstructuredGrid", outputType="vtkTable" ) + self.counts: CellTypeCounts + + 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" ) + return 1 + + 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 = self.GetInputData( inInfoVec, 0, 0 ) + outTable: vtkTable = vtkTable.GetData( outInfoVec, 0 ) + assert inData is not None, "Input mesh is undefined." + assert outTable is not None, "Output table is undefined." + + # compute cell type counts + self.counts = CellTypeCounts() + self.counts.setTypeCount( VTK_VERTEX, inData.GetNumberOfPoints() ) + for i in range( inData.GetNumberOfCells() ): + cell: vtkCell = inData.GetCell( i ) + self.counts.addType( cell.GetCellType() ) + + # create output table + # first reset output table + outTable.RemoveAllRows() + outTable.RemoveAllColumns() + outTable.SetNumberOfRows( 1 ) + + # create columns per types + for cellType in self.getAllCellTypes(): + array: vtkIntArray = vtkIntArray() + array.SetName( vtkCellTypes.GetClassNameFromTypeId( cellType ) ) + array.SetNumberOfComponents( 1 ) + array.SetNumberOfValues( 1 ) + array.SetValue( 0, self.counts.getTypeCount( cellType ) ) + outTable.AddColumn( array ) + return 1 + + def GetCellTypeCounts( self: Self ) -> CellTypeCounts: + """Get CellTypeCounts object. + + Returns: + CellTypeCounts: CellTypeCounts object. + """ + return self.counts + + def getAllCellTypes( self: Self ) -> tuple[ int, ...]: + """Get all cell type ids managed by CellTypeCount class. + + Returns: + tuple[int,...]: tuple containg cell type ids. + """ + return ( VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON ) diff --git a/geos-mesh/src/geos/mesh/utils/arrayHelpers.py b/geos-mesh/src/geos/mesh/utils/arrayHelpers.py index 2ca957fa..3139d67f 100644 --- a/geos-mesh/src/geos/mesh/utils/arrayHelpers.py +++ b/geos-mesh/src/geos/mesh/utils/arrayHelpers.py @@ -26,48 +26,38 @@ """ -def has_invalid_field( mesh: vtkUnstructuredGrid, invalid_fields: list[ str ] ) -> bool: - """Checks if a mesh contains at least a data arrays within its cell, field or point data - having a certain name. If so, returns True, else False. +def has_array( mesh: vtkUnstructuredGrid, array_names: list[ str ] ) -> bool: + """Checks if input mesh contains at least one of input data arrays. Args: mesh (vtkUnstructuredGrid): An unstructured mesh. - invalid_fields (list[str]): Field name of an array in any data from the data. + array_names (list[str]): List of array names. Returns: - bool: True if one field found, else False. + bool: True if at least one array is found, else False. """ # Check the cell data fields - cell_data = mesh.GetCellData() - for i in range( cell_data.GetNumberOfArrays() ): - if cell_data.GetArrayName( i ) in invalid_fields: - logging.error( f"The mesh contains an invalid cell field name '{cell_data.GetArrayName( i )}'." ) - return True - # Check the field data fields - field_data = mesh.GetFieldData() - for i in range( field_data.GetNumberOfArrays() ): - if field_data.GetArrayName( i ) in invalid_fields: - logging.error( f"The mesh contains an invalid field name '{field_data.GetArrayName( i )}'." ) - return True - # Check the point data fields - point_data = mesh.GetPointData() - for i in range( point_data.GetNumberOfArrays() ): - if point_data.GetArrayName( i ) in invalid_fields: - logging.error( f"The mesh contains an invalid point field name '{point_data.GetArrayName( i )}'." ) - return True + data: vtkFieldData | None + for data in ( mesh.GetCellData(), mesh.GetFieldData(), mesh.GetPointData() ): + if data is None: + continue # type: ignore[unreachable] + for arrayName in array_names: + if data.HasArray( arrayName ): + logging.error( f"The mesh contains the array named '{arrayName}'." ) + return True return False def getFieldType( data: vtkFieldData ) -> str: - """A vtk grid can contain 3 types of field data: + """Returns whether the data is "vtkFieldData", "vtkCellData" or "vtkPointData". + + A vtk mesh can contain 3 types of field data: - vtkFieldData (parent class) - vtkCellData (inheritance of vtkFieldData) - vtkPointData (inheritance of vtkFieldData) - The goal is to return whether the data is "vtkFieldData", "vtkCellData" or "vtkPointData". - Args: - data (vtkFieldData) + data (vtkFieldData): vtk field data Returns: str: "vtkFieldData", "vtkCellData" or "vtkPointData" @@ -86,7 +76,7 @@ def getArrayNames( data: vtkFieldData ) -> list[ str ]: """Get the names of all arrays stored in a "vtkFieldData", "vtkCellData" or "vtkPointData". Args: - data (vtkFieldData) + data (vtkFieldData): vtk field data Returns: list[ str ]: The array names in the order that they are stored in the field data. @@ -100,8 +90,9 @@ def getArrayByName( data: vtkFieldData, name: str ) -> Optional[ vtkDataArray ]: """Get the vtkDataArray corresponding to the given name. Args: - data (vtkFieldData) - name (str) + data (vtkFieldData): vtk field data + name (str): array name + Returns: Optional[ vtkDataArray ]: The vtkDataArray associated with the name given. None if not found. @@ -116,8 +107,9 @@ def getCopyArrayByName( data: vtkFieldData, name: str ) -> Optional[ vtkDataArra """Get the copy of a vtkDataArray corresponding to the given name. Args: - data (vtkFieldData) - name (str) + data (vtkFieldData): vtk field data + name (str): array name + Returns: Optional[ vtkDataArray ]: The copy of the vtkDataArray associated with the name given. None if not found. @@ -132,7 +124,8 @@ def getNumpyGlobalIdsArray( data: Union[ vtkCellData, vtkPointData ] ) -> Option """Get a numpy array of the GlobalIds. Args: - data (Union[ vtkCellData, vtkPointData ]) + data (Union[ vtkCellData, vtkPointData ]): Cell or point array. + Returns: Optional[ npt.NDArray[ np.int64 ] ]: The numpy array of GlobalIds. @@ -144,26 +137,25 @@ def getNumpyGlobalIdsArray( data: Union[ vtkCellData, vtkPointData ] ) -> Option return vtk_to_numpy( global_ids ) -def getNumpyArrayByName( data: vtkFieldData, name: str, sorted: bool = False ) -> Optional[ npt.NDArray ]: +def getNumpyArrayByName( data: vtkCellData | vtkPointData, name: str, sorted: bool = False ) -> Optional[ npt.NDArray ]: """Get the numpy array of a given vtkDataArray found by its name. + If sorted is selected, this allows the option to reorder the values wrt GlobalIds. If not GlobalIds was found, no reordering will be perform. Args: - data (vtkFieldData) - name (str) + data (vtkCellData | vtkPointData): vtk field data. + name (str): Array name to sort sorted (bool, optional): Sort the output array with the help of GlobalIds. Defaults to False. Returns: - Optional[ npt.NDArray ] + Optional[ npt.NDArray ]: Sorted array """ dataArray: Optional[ vtkDataArray ] = getArrayByName( data, name ) if dataArray is not None: - arr: Optional[ npt.NDArray ] = vtk_to_numpy( dataArray ) - if sorted: - fieldType: str = getFieldType( data ) - if fieldType in [ "vtkCellData", "vtkPointData" ]: - sortArrayByGlobalIds( data, arr ) + arr: npt.NDArray[ np.float64 ] = vtk_to_numpy( dataArray ) + if sorted and ( data.IsA( "vtkCellData" ) or data.IsA( "vtkPointData" ) ): + sortArrayByGlobalIds( data, arr ) return arr return None @@ -491,7 +483,8 @@ def getComponentNamesDataSet( dataSet: vtkDataSet, attributeName: str, onPoints: """ array: vtkDoubleArray = getVtkArrayInObject( dataSet, attributeName, onPoints ) - componentNames: list[ str ] = list() + componentNames: list[ str ] = [] + if array.GetNumberOfComponents() > 1: componentNames += [ array.GetComponentName( i ) for i in range( array.GetNumberOfComponents() ) ] return tuple( componentNames ) @@ -658,12 +651,12 @@ def computeCellCenterCoordinates( mesh: vtkDataSet ) -> vtkDataArray: return pts.GetData() -def sortArrayByGlobalIds( data: Union[ vtkCellData, vtkFieldData ], arr: npt.NDArray[ np.int64 ] ) -> None: - """Sort an array following global Ids +def sortArrayByGlobalIds( data: Union[ vtkCellData, vtkPointData ], arr: npt.NDArray[ np.float64 ] ) -> None: + """Sort an array following global Ids. Args: data (vtkFieldData): Global Ids array - arr (npt.NDArray[ np.int64 ]): Array to sort + arr (npt.NDArray[ np.float64 ]): Array to sort """ globalids: Optional[ npt.NDArray[ np.int64 ] ] = getNumpyGlobalIdsArray( data ) if globalids is not None: diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py index 27005395..2faf092d 100644 --- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py +++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py @@ -1,10 +1,14 @@ # SPDX-License-Identifier: Apache-2.0 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. # SPDX-FileContributor: Martin Lemay, Paloma Martinez -from typing import Any, Iterator, List -from vtkmodules.vtkCommonCore import vtkIdList -from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkPolyData, vtkPlane +import numpy as np +import numpy.typing as npt +from typing import Iterator, List, Sequence, Any, Union +from vtkmodules.util.numpy_support import numpy_to_vtk +from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints, reference +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkMultiBlockDataSet, vtkPolyData, vtkDataSet, vtkDataObject, vtkPlane, vtkCellTypes, vtkIncrementalOctreePointLocator from vtkmodules.vtkFiltersCore import vtk3DLinearGridPlaneCutter +from geos.mesh.utils.multiblockHelpers import ( getBlockElementIndexesFlatten, getBlockFromFlatIndex ) __doc__ = """ Generic VTK utilities. @@ -32,21 +36,19 @@ def to_vtk_id_list( data: List[ int ] ) -> vtkIdList: return result -def vtk_iter( vtkContainer ) -> Iterator[ Any ]: - """ - Utility function transforming a vtk "container" (e.g. vtkIdList) into an iterable to be used for building built-ins - python containers. +def vtk_iter( vtkContainer: vtkIdList | vtkCellTypes ) -> Iterator[ Any ]: + """Utility function transforming a vtk "container" into an iterable. Args: - vtkContainer: A vtk container + vtkContainer (vtkIdList | vtkCellTypes): A vtk container Returns: - The iterator + Iterator[ Any ]: The iterator """ - if hasattr( vtkContainer, "GetNumberOfIds" ): + if isinstance( vtkContainer, vtkIdList ): for i in range( vtkContainer.GetNumberOfIds() ): yield vtkContainer.GetId( i ) - elif hasattr( vtkContainer, "GetNumberOfTypes" ): + elif isinstance( vtkContainer, vtkCellTypes ): for i in range( vtkContainer.GetNumberOfTypes() ): yield vtkContainer.GetCellType( i ) @@ -81,3 +83,187 @@ def extractSurfaceFromElevation( mesh: vtkUnstructuredGrid, elevation: float ) - cutter.SetInterpolateAttributes( True ) cutter.Update() return cutter.GetOutputDataObject( 0 ) + + +def getBounds( + input: Union[ vtkUnstructuredGrid, + vtkMultiBlockDataSet ] ) -> tuple[ float, float, float, float, float, float ]: + """Get bounds of either single of composite data set. + + Args: + input (Union[vtkUnstructuredGrid, vtkMultiBlockDataSet]): input mesh + + Returns: + tuple[float, float, float, float, float, float]: tuple containing + bounds (xmin, xmax, ymin, ymax, zmin, zmax) + + """ + if isinstance( input, vtkMultiBlockDataSet ): + return getMultiBlockBounds( input ) + else: + return getMonoBlockBounds( input ) + + +def getMonoBlockBounds( input: vtkUnstructuredGrid, ) -> tuple[ float, float, float, float, float, float ]: + """Get boundary box extrema coordinates for a vtkUnstructuredGrid. + + Args: + input (vtkMultiBlockDataSet): input single block mesh + + Returns: + tuple[float, float, float, float, float, float]: tuple containing + bounds (xmin, xmax, ymin, ymax, zmin, zmax) + + """ + return input.GetBounds() + + +def getMultiBlockBounds( input: vtkMultiBlockDataSet, ) -> tuple[ float, float, float, float, float, float ]: + """Get boundary box extrema coordinates for a vtkMultiBlockDataSet. + + Args: + input (vtkMultiBlockDataSet): input multiblock mesh + + Returns: + tuple[float, float, float, float, float, float]: bounds. + + """ + xmin, ymin, zmin = 3 * [ np.inf ] + xmax, ymax, zmax = 3 * [ -1.0 * np.inf ] + blockIndexes: list[ int ] = getBlockElementIndexesFlatten( input ) + for blockIndex in blockIndexes: + block0: vtkDataObject = getBlockFromFlatIndex( input, blockIndex ) + assert block0 is not None, "Mesh is undefined." + block: vtkDataSet = vtkDataSet.SafeDownCast( block0 ) + bounds: tuple[ float, float, float, float, float, float ] = block.GetBounds() + xmin = bounds[ 0 ] if bounds[ 0 ] < xmin else xmin + xmax = bounds[ 1 ] if bounds[ 1 ] > xmax else xmax + ymin = bounds[ 2 ] if bounds[ 2 ] < ymin else ymin + ymax = bounds[ 3 ] if bounds[ 3 ] > ymax else ymax + zmin = bounds[ 4 ] if bounds[ 4 ] < zmin else zmin + zmax = bounds[ 5 ] if bounds[ 5 ] > zmax else zmax + return xmin, xmax, ymin, ymax, zmin, zmax + + +def getBoundsFromPointCoords( cellPtsCoord: list[ npt.NDArray[ np.float64 ] ] ) -> Sequence[ float ]: + """Compute bounding box coordinates of the list of points. + + Args: + cellPtsCoord (list[npt.NDArray[np.float64]]): list of points + + Returns: + Sequence[float]: bounding box coordinates (xmin, xmax, ymin, ymax, zmin, zmax) + """ + bounds: list[ float ] = [ + np.inf, + -np.inf, + np.inf, + -np.inf, + np.inf, + -np.inf, + ] + for ptsCoords in cellPtsCoord: + mins: npt.NDArray[ np.float64 ] = np.min( ptsCoords, axis=0 ) + maxs: npt.NDArray[ np.float64 ] = np.max( ptsCoords, axis=0 ) + for i in range( 3 ): + bounds[ 2 * i ] = float( min( bounds[ 2 * i ], mins[ i ] ) ) + bounds[ 2 * i + 1 ] = float( max( bounds[ 2 * i + 1 ], maxs[ i ] ) ) + return bounds + + +def createSingleCellMesh( cellType: int, ptsCoord: npt.NDArray[ np.float64 ] ) -> vtkUnstructuredGrid: + """Create a mesh that consists of a single cell. + + Args: + cellType (int): cell type + ptsCoord (1DArray[np.float64]): cell point coordinates + + Returns: + vtkUnstructuredGrid: output mesh + """ + nbPoints: int = ptsCoord.shape[ 0 ] + points: npt.NDArray[ np.float64 ] = np.vstack( ( ptsCoord, ) ) + # Convert points to vtkPoints object + vtkpts: vtkPoints = vtkPoints() + vtkpts.SetData( numpy_to_vtk( points ) ) + + # create cells from point ids + cellsID: vtkIdList = vtkIdList() + for j in range( nbPoints ): + cellsID.InsertNextId( j ) + + # add cell to mesh + mesh: vtkUnstructuredGrid = vtkUnstructuredGrid() + mesh.SetPoints( vtkpts ) + mesh.Allocate( 1 ) + mesh.InsertNextCell( cellType, cellsID ) + return mesh + + +def createMultiCellMesh( cellTypes: list[ int ], + cellPtsCoord: list[ npt.NDArray[ np.float64 ] ], + sharePoints: bool = True ) -> vtkUnstructuredGrid: + """Create a mesh that consists of multiple cells. + + .. WARNING:: the mesh is not check for conformity. + + Args: + cellTypes (list[int]): cell type + cellPtsCoord (list[1DArray[np.float64]]): list of cell point coordinates + sharePoints (bool): if True, cells share points, else a new point is created for each cell vertex + + Returns: + vtkUnstructuredGrid: output mesh + """ + assert len( cellPtsCoord ) == len( cellTypes ), "The lists of cell types of point coordinates must be of same size." + nbCells: int = len( cellPtsCoord ) + mesh: vtkUnstructuredGrid = vtkUnstructuredGrid() + points: vtkPoints + cellVertexMapAll: list[ tuple[ int, ...] ] + points, cellVertexMapAll = createVertices( cellPtsCoord, sharePoints ) + assert len( cellVertexMapAll ) == len( + cellTypes ), "The lists of cell types of cell point ids must be of same size." + mesh.SetPoints( points ) + mesh.Allocate( nbCells ) + # create mesh cells + for cellType, ptsId in zip( cellTypes, cellVertexMapAll, strict=True ): + # create cells from point ids + cellsID: vtkIdList = vtkIdList() + for ptId in ptsId: + cellsID.InsertNextId( ptId ) + mesh.InsertNextCell( cellType, cellsID ) + return mesh + + +def createVertices( cellPtsCoord: list[ npt.NDArray[ np.float64 ] ], + shared: bool = True ) -> tuple[ vtkPoints, list[ tuple[ int, ...] ] ]: + """Create vertices from cell point coordinates list. + + Args: + cellPtsCoord (list[npt.NDArray[np.float64]]): list of cell point coordinates + shared (bool, optional): If True, collocated points are merged. Defaults to True. + + Returns: + tuple[vtkPoints, list[tuple[int, ...]]]: tuple containing points and the + map of cell point ids + """ + # get point bounds + bounds: Sequence[ float ] = getBoundsFromPointCoords( cellPtsCoord ) + points: vtkPoints = vtkPoints() + # use point locator to check for colocated points + pointsLocator = vtkIncrementalOctreePointLocator() + pointsLocator.InitPointInsertion( points, bounds ) + cellVertexMapAll: list[ tuple[ int, ...] ] = [] + ptId: reference = reference( 0 ) + ptsCoords: npt.NDArray[ np.float64 ] + for ptsCoords in cellPtsCoord: + cellVertexMap: list[ int ] = [] + pt: npt.NDArray[ np.float64 ] # 1DArray + for pt in ptsCoords: + if shared: + pointsLocator.InsertUniquePoint( pt.tolist(), ptId ) # type: ignore[arg-type] + else: + pointsLocator.InsertPointWithoutChecking( pt.tolist(), ptId, 1 ) # type: ignore[arg-type] + cellVertexMap += [ ptId.get() ] # type: ignore + cellVertexMapAll += [ tuple( cellVertexMap ) ] + return points, cellVertexMapAll diff --git a/geos-mesh/tests/data/hexa_cell.csv b/geos-mesh/tests/data/hexa_cell.csv new file mode 100644 index 00000000..741fa4f7 --- /dev/null +++ b/geos-mesh/tests/data/hexa_cell.csv @@ -0,0 +1,8 @@ +0.0, 0.0, 0.0 +1.0, 0.0, 0.0 +1.0, 1.0, 0.0 +0.0, 1.0, 0.0 +0.0, 0.0, 1.0 +1.0, 0.0, 1.0 +1.0, 1.0, 1.0 +0.0, 1.0, 1.0 \ No newline at end of file diff --git a/geos-mesh/tests/data/hexa_mesh.csv b/geos-mesh/tests/data/hexa_mesh.csv new file mode 100644 index 00000000..cc55f562 --- /dev/null +++ b/geos-mesh/tests/data/hexa_mesh.csv @@ -0,0 +1,64 @@ +0.0,0.0,0.5 +0.5,0.0,0.5 +0.5,0.5,0.5 +0.0,0.5,0.5 +0.0,0.0,1.0 +0.5,0.0,1.0 +0.5,0.5,1.0 +0.0,0.5,1.0 +0.5,0.0,0.5 +1.0,0.0,0.5 +1.0,0.5,0.5 +0.5,0.5,0.5 +0.5,0.0,1.0 +1.0,0.0,1.0 +1.0,0.5,1.0 +0.5,0.5,1.0 +0.0,0.0,0.0 +0.5,0.0,0.0 +0.5,0.5,0.0 +0.0,0.5,0.0 +0.0,0.0,0.5 +0.5,0.0,0.5 +0.5,0.5,0.5 +0.0,0.5,0.5 +0.5,0.0,0.0 +1.0,0.0,0.0 +1.0,0.5,0.0 +0.5,0.5,0.0 +0.5,0.0,0.5 +1.0,0.0,0.5 +1.0,0.5,0.5 +0.5,0.5,0.5 +0.0,0.5,0.5 +0.5,0.5,0.5 +0.5,1.0,0.5 +0.0,1.0,0.5 +0.0,0.5,1.0 +0.5,0.5,1.0 +0.5,1.0,1.0 +0.0,1.0,1.0 +0.5,0.5,0.5 +1.0,0.5,0.5 +1.0,1.0,0.5 +0.5,1.0,0.5 +0.5,0.5,1.0 +1.0,0.5,1.0 +1.0,1.0,1.0 +0.5,1.0,1.0 +0.0,0.5,0.0 +0.5,0.5,0.0 +0.5,1.0,0.0 +0.0,1.0,0.0 +0.0,0.5,0.5 +0.5,0.5,0.5 +0.5,1.0,0.5 +0.0,1.0,0.5 +0.5,0.5,0.0 +1.0,0.5,0.0 +1.0,1.0,0.0 +0.5,1.0,0.0 +0.5,0.5,0.5 +1.0,0.5,0.5 +1.0,1.0,0.5 +0.5,1.0,0.5 diff --git a/geos-mesh/tests/data/pyramid_cell.csv b/geos-mesh/tests/data/pyramid_cell.csv new file mode 100644 index 00000000..864deaf4 --- /dev/null +++ b/geos-mesh/tests/data/pyramid_cell.csv @@ -0,0 +1,5 @@ +0.0, 0.0, 0.0 +1.0, 0.0, 0.0 +1.0, 1.0, 0.0 +0.0, 1.0, 0.0 +0.5, 0.5, 1.0 \ No newline at end of file diff --git a/geos-mesh/tests/data/pyramid_mesh.csv b/geos-mesh/tests/data/pyramid_mesh.csv new file mode 100644 index 00000000..c435d2e4 --- /dev/null +++ b/geos-mesh/tests/data/pyramid_mesh.csv @@ -0,0 +1,46 @@ +0.5,0.0,0.0 +1.0,0.0,0.0 +1.0,0.5,0.0 +0.5,0.5,0.0 +0.8,0.2,0.5 +0.5,0.5,0.0 +1.0,0.5,0.0 +1.0,1.0,0.0 +0.5,1.0,0.0 +0.8,0.8,0.5 +0.0,1.0,0.0 +0.0,0.5,0.0 +0.5,0.5,0.0 +0.5,1.0,0.0 +0.2,0.8,0.5 +0.0,0.5,0.0 +0.0,0.0,0.0 +0.5,0.0,0.0 +0.5,0.5,0.0 +0.2,0.2,0.5 +0.2,0.8,0.5 +0.2,0.2,0.5 +0.8,0.2,0.5 +0.8,0.8,0.5 +0.5,0.5,1.0 +0.8,0.8,0.5 +0.8,0.2,0.5 +0.2,0.2,0.5 +0.2,0.8,0.5 +0.5,0.5,0.0 +0.2,0.2,0.5 +0.8,0.2,0.5 +0.5,0.0,0.0 +0.5,0.5,0.0 +0.8,0.2,0.5 +0.8,0.8,0.5 +1.0,0.5,0.0 +0.5,0.5,0.0 +0.8,0.8,0.5 +0.2,0.8,0.5 +0.5,1.0,0.0 +0.5,0.5,0.0 +0.2,0.8,0.5 +0.2,0.2,0.5 +0.0,0.5,0.0 +0.5,0.5,0.0 diff --git a/geos-mesh/tests/data/quad_cell.csv b/geos-mesh/tests/data/quad_cell.csv new file mode 100644 index 00000000..ffca9522 --- /dev/null +++ b/geos-mesh/tests/data/quad_cell.csv @@ -0,0 +1,4 @@ +0.0, 0.0, 0.0 +1.0, 0.0, 0.0 +1.0, 1.0, 0.0 +0.0, 1.0, 0.0 \ No newline at end of file diff --git a/geos-mesh/tests/data/tetra_cell.csv b/geos-mesh/tests/data/tetra_cell.csv new file mode 100644 index 00000000..38b971aa --- /dev/null +++ b/geos-mesh/tests/data/tetra_cell.csv @@ -0,0 +1,4 @@ +0.0, 0.0, 0.0 +1.0, 0.0, 0.0 +0.0, 0.0, 1.0 +0.0, 1.0, 0.0 \ No newline at end of file diff --git a/geos-mesh/tests/data/tetra_mesh.csv b/geos-mesh/tests/data/tetra_mesh.csv new file mode 100644 index 00000000..2f3414b4 --- /dev/null +++ b/geos-mesh/tests/data/tetra_mesh.csv @@ -0,0 +1,32 @@ +0.0,0.0,0.0 +0.5,0.0,0.0 +0.0,0.0,0.5 +0.0,0.5,0.0 +0.0,0.5,0.0 +0.5,0.5,0.0 +0.0,0.5,0.5 +0.0,1.0,0.0 +0.5,0.5,0.0 +0.5,0.0,0.0 +0.5,0.0,0.5 +1.0,0.0,0.0 +0.5,0.0,0.5 +0.0,0.0,0.5 +0.0,0.5,0.5 +0.0,0.0,1.0 +0.0,0.0,0.5 +0.0,0.5,0.5 +0.0,0.5,0.0 +0.5,0.0,0.0 +0.5,0.0,0.0 +0.0,0.5,0.5 +0.0,0.5,0.0 +0.5,0.5,0.0 +0.5,0.0,0.0 +0.0,0.5,0.5 +0.5,0.5,0.0 +0.5,0.0,0.5 +0.5,0.0,0.5 +0.5,0.0,0.0 +0.0,0.5,0.5 +0.0,0.0,0.5 diff --git a/geos-mesh/tests/data/triangle_cell.csv b/geos-mesh/tests/data/triangle_cell.csv new file mode 100644 index 00000000..b8b70c27 --- /dev/null +++ b/geos-mesh/tests/data/triangle_cell.csv @@ -0,0 +1,3 @@ +0.0, 0.0, 0.0 +1.0, 0.0, 0.0 +0.0, 1.0, 0.0 \ No newline at end of file diff --git a/geos-mesh/tests/test_CellTypeCounter.py b/geos-mesh/tests/test_CellTypeCounter.py new file mode 100644 index 00000000..6da50c5c --- /dev/null +++ b/geos-mesh/tests/test_CellTypeCounter.py @@ -0,0 +1,159 @@ +# SPDX-FileContributor: Martin Lemay +# SPDX-License-Identifier: Apache 2.0 +# ruff: noqa: E402 # disable Module level import not at top of file +import os +from dataclasses import dataclass +import numpy as np +import numpy.typing as npt +import pytest +from typing import ( + Iterator, ) + +from geos.mesh.utils.genericHelpers import createSingleCellMesh, createMultiCellMesh +from geos.mesh.stats.CellTypeCounter import CellTypeCounter +from geos.mesh.model.CellTypeCounts import CellTypeCounts + +from vtkmodules.vtkCommonDataModel import ( + vtkUnstructuredGrid, + vtkCellTypes, + vtkCell, + VTK_TRIANGLE, + VTK_QUAD, + VTK_TETRA, + VTK_VERTEX, + VTK_POLYHEDRON, + VTK_POLYGON, + VTK_PYRAMID, + VTK_HEXAHEDRON, + VTK_WEDGE, +) + +data_root: str = os.path.join( os.path.dirname( os.path.abspath( __file__ ) ), "data" ) + +filename_all: tuple[ str, ...] = ( "triangle_cell.csv", "quad_cell.csv", "tetra_cell.csv", "pyramid_cell.csv", + "hexa_cell.csv" ) +cellType_all: tuple[ int, ...] = ( VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_HEXAHEDRON ) + +filename_all2: tuple[ str, ...] = ( "tetra_mesh.csv", "hexa_mesh.csv" ) +cellType_all2: tuple[ int, ...] = ( VTK_TETRA, VTK_HEXAHEDRON ) +nbPtsCell_all2: tuple[ int ] = ( 4, 8 ) + + +@dataclass( frozen=True ) +class TestCase: + """Test case.""" + __test__ = False + #: mesh + mesh: vtkUnstructuredGrid + + +def __generate_test_data_single_cell() -> Iterator[ TestCase ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: iterator on test cases + """ + for cellType, filename in zip( cellType_all, filename_all, strict=True ): + ptsCoord: npt.NDArray[ np.float64 ] = np.loadtxt( os.path.join( data_root, filename ), + dtype=float, + delimiter=',' ) + mesh: vtkUnstructuredGrid = createSingleCellMesh( cellType, ptsCoord ) + yield TestCase( mesh ) + + +ids: list[ str ] = [ vtkCellTypes.GetClassNameFromTypeId( cellType ) for cellType in cellType_all ] + + +@pytest.mark.parametrize( "test_case", __generate_test_data_single_cell(), ids=ids ) +def test_CellTypeCounter_single( test_case: TestCase ) -> None: + """Test of CellTypeCounter filter. + + Args: + test_case (TestCase): test case + """ + filter: CellTypeCounter = CellTypeCounter() + filter.SetInputDataObject( test_case.mesh ) + filter.Update() + countsObs: CellTypeCounts = filter.GetCellTypeCounts() + assert countsObs is not None, "CellTypeCounts is undefined" + + assert countsObs.getTypeCount( VTK_VERTEX ) == test_case.mesh.GetNumberOfPoints( + ), f"Number of vertices should be {test_case.mesh.GetNumberOfPoints()}" + + # compute counts for each type of cell + elementTypes: tuple[ int ] = ( VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_HEXAHEDRON, VTK_WEDGE ) + counts: npt.NDArray[ np.int64 ] = np.zeros( len( elementTypes ) ) + for i in range( test_case.mesh.GetNumberOfCells() ): + cell: vtkCell = test_case.mesh.GetCell( i ) + index: int = elementTypes.index( cell.GetCellType() ) + counts[ index ] += 1 + # check cell type counts + for i, elementType in enumerate( elementTypes ): + assert int( + countsObs.getTypeCount( elementType ) + ) == counts[ i ], f"The number of {vtkCellTypes.GetClassNameFromTypeId(elementType)} should be {counts[i]}." + + nbPolygon: int = counts[ 0 ] + counts[ 1 ] + nbPolyhedra: int = np.sum( counts[ 2: ] ) + assert int( countsObs.getTypeCount( VTK_POLYGON ) ) == nbPolygon, f"The number of faces should be {nbPolygon}." + assert int( + countsObs.getTypeCount( VTK_POLYHEDRON ) ) == nbPolyhedra, f"The number of polyhedra should be {nbPolyhedra}." + + +def __generate_test_data_multi_cell() -> Iterator[ TestCase ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: iterator on test cases + """ + for cellType, filename, nbPtsCell in zip( cellType_all2, filename_all2, nbPtsCell_all2, strict=True ): + ptsCoords: npt.NDArray[ np.float64 ] = np.loadtxt( os.path.join( data_root, filename ), + dtype=float, + delimiter=',' ) + # split array to get a list of coordinates per cell + cellPtsCoords: list[ npt.NDArray[ np.float64 ] ] = [ + ptsCoords[ i:i + nbPtsCell ] for i in range( 0, ptsCoords.shape[ 0 ], nbPtsCell ) + ] + nbCells: int = int( ptsCoords.shape[ 0 ] / nbPtsCell ) + cellTypes = nbCells * [ cellType ] + mesh: vtkUnstructuredGrid = createMultiCellMesh( cellTypes, cellPtsCoords, True ) + yield TestCase( mesh ) + + +ids2: list[ str ] = [ os.path.splitext( name )[ 0 ] for name in filename_all2 ] + + +@pytest.mark.parametrize( "test_case", __generate_test_data_multi_cell(), ids=ids2 ) +def test_CellTypeCounter_multi( test_case: TestCase ) -> None: + """Test of CellTypeCounter filter. + + Args: + test_case (TestCase): test case + """ + filter: CellTypeCounter = CellTypeCounter() + filter.SetInputDataObject( test_case.mesh ) + filter.Update() + countsObs: CellTypeCounts = filter.GetCellTypeCounts() + assert countsObs is not None, "CellTypeCounts is undefined" + + assert countsObs.getTypeCount( VTK_VERTEX ) == test_case.mesh.GetNumberOfPoints( + ), f"Number of vertices should be {test_case.mesh.GetNumberOfPoints()}" + + # compute counts for each type of cell + elementTypes: tuple[ int ] = ( VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_HEXAHEDRON, VTK_WEDGE ) + counts: npt.NDArray[ np.int64 ] = np.zeros( len( elementTypes ), dtype=int ) + for i in range( test_case.mesh.GetNumberOfCells() ): + cell: vtkCell = test_case.mesh.GetCell( i ) + index: int = elementTypes.index( cell.GetCellType() ) + counts[ index ] += 1 + # check cell type counts + for i, elementType in enumerate( elementTypes ): + assert int( + countsObs.getTypeCount( elementType ) + ) == counts[ i ], f"The number of {vtkCellTypes.GetClassNameFromTypeId(elementType)} should be {counts[i]}." + + nbPolygon: int = counts[ 0 ] + counts[ 1 ] + nbPolyhedra: int = np.sum( counts[ 2: ] ) + assert int( countsObs.getTypeCount( VTK_POLYGON ) ) == nbPolygon, f"The number of faces should be {nbPolygon}." + assert int( + countsObs.getTypeCount( VTK_POLYHEDRON ) ) == nbPolyhedra, f"The number of polyhedra should be {nbPolyhedra}." diff --git a/geos-mesh/tests/test_CellTypeCounts.py b/geos-mesh/tests/test_CellTypeCounts.py new file mode 100644 index 00000000..32d3c7ea --- /dev/null +++ b/geos-mesh/tests/test_CellTypeCounts.py @@ -0,0 +1,220 @@ +# SPDX-FileContributor: Martin Lemay +# SPDX-License-Identifier: Apache 2.0 +# ruff: noqa: E402 # disable Module level import not at top of file +from dataclasses import dataclass +import pytest +from typing import ( + Iterator, ) + +from geos.mesh.model.CellTypeCounts import CellTypeCounts + +from vtkmodules.vtkCommonDataModel import ( vtkCellTypes, VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, + VTK_HEXAHEDRON, VTK_WEDGE, VTK_VERTEX ) + +# inputs +nbVertex_all: tuple[ int ] = ( 3, 4, 5, 8, 10, 20 ) +nbTri_all: tuple[ int ] = ( 1, 0, 3, 0, 0, 4 ) +nbQuad_all: tuple[ int ] = ( 0, 1, 0, 6, 0, 3 ) +nbTetra_all: tuple[ int ] = ( 0, 0, 1, 0, 4, 0 ) +nbPyr_all: tuple[ int ] = ( 0, 0, 0, 0, 0, 4 ) +nbWed_all: tuple[ int ] = ( 0, 0, 0, 0, 0, 2 ) +nbHexa_all: tuple[ int ] = ( 0, 0, 0, 1, 0, 5 ) + + +@dataclass( frozen=True ) +class TestCase: + """Test case.""" + __test__ = False + nbVertex: tuple[ int ] + nbTri: tuple[ int ] + nbQuad: tuple[ int ] + nbTetra: tuple[ int ] + nbPyr: tuple[ int ] + nbWed: tuple[ int ] + nbHexa: tuple[ int ] + + +def __generate_test_data() -> Iterator[ TestCase ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: iterator on test cases + """ + for nbVertex, nbTri, nbQuad, nbTetra, nbPyr, nbWed, nbHexa in zip( nbVertex_all, + nbTri_all, + nbQuad_all, + nbTetra_all, + nbPyr_all, + nbWed_all, + nbHexa_all, + strict=True ): + yield TestCase( nbVertex, nbTri, nbQuad, nbTetra, nbPyr, nbWed, nbHexa ) + + +def __get_expected_counts( + nbVertex: int, + nbTri: int, + nbQuad: int, + nbTetra: int, + nbPyr: int, + nbWed: int, + nbHexa: int, +) -> str: + nbFaces: int = nbTri + nbQuad + nbPolyhedre: int = nbTetra + nbPyr + nbHexa + nbWed + countsExp: str = "" + countsExp += "| | |\n" + countsExp += "| - | - |\n" + countsExp += f"| **Total Number of Vertices** | {int(nbVertex):12} |\n" + countsExp += f"| **Total Number of Polygon** | {int(nbFaces):12} |\n" + countsExp += f"| **Total Number of Polyhedron** | {int(nbPolyhedre):12} |\n" + countsExp += f"| **Total Number of Cells** | {int(nbPolyhedre+nbFaces):12} |\n" + countsExp += "| - | - |\n" + for cellType, nb in zip( ( + VTK_TRIANGLE, + VTK_QUAD, + ), ( + nbTri, + nbQuad, + ), strict=True ): + countsExp += f"| **Total Number of {vtkCellTypes.GetClassNameFromTypeId(cellType):<13}** | {int(nb):12} |\n" + for cellType, nb in zip( ( VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON ), ( nbTetra, nbPyr, nbWed, nbHexa ), + strict=True ): + countsExp += f"| **Total Number of {vtkCellTypes.GetClassNameFromTypeId(cellType):<13}** | {int(nb):12} |\n" + return countsExp + + +def test_CellTypeCounts_init() -> None: + """Test of CellTypeCounts . + + Args: + test_case (TestCase): test case + """ + counts: CellTypeCounts = CellTypeCounts() + assert counts.getTypeCount( VTK_VERTEX ) == 0, "Number of vertices must be 0" + assert counts.getTypeCount( VTK_TRIANGLE ) == 0, "Number of triangles must be 0" + assert counts.getTypeCount( VTK_QUAD ) == 0, "Number of quads must be 0" + assert counts.getTypeCount( VTK_TETRA ) == 0, "Number of tetrahedra must be 0" + assert counts.getTypeCount( VTK_PYRAMID ) == 0, "Number of pyramids must be 0" + assert counts.getTypeCount( VTK_WEDGE ) == 0, "Number of wedges must be 0" + assert counts.getTypeCount( VTK_HEXAHEDRON ) == 0, "Number of hexahedra must be 0" + + +@pytest.mark.parametrize( "test_case", __generate_test_data() ) +def test_CellTypeCounts_addType( test_case: TestCase ) -> None: + """Test of CellTypeCounts . + + Args: + test_case (TestCase): test case + """ + counts: CellTypeCounts = CellTypeCounts() + for _ in range( test_case.nbVertex ): + counts.addType( VTK_VERTEX ) + for _ in range( test_case.nbTri ): + counts.addType( VTK_TRIANGLE ) + for _ in range( test_case.nbQuad ): + counts.addType( VTK_QUAD ) + for _ in range( test_case.nbTetra ): + counts.addType( VTK_TETRA ) + for _ in range( test_case.nbPyr ): + counts.addType( VTK_PYRAMID ) + for _ in range( test_case.nbWed ): + counts.addType( VTK_WEDGE ) + for _ in range( test_case.nbHexa ): + counts.addType( VTK_HEXAHEDRON ) + + assert counts.getTypeCount( VTK_VERTEX ) == test_case.nbVertex, f"Number of vertices must be {test_case.nbVertex}" + assert counts.getTypeCount( VTK_TRIANGLE ) == test_case.nbTri, f"Number of triangles must be {test_case.nbTri}" + assert counts.getTypeCount( VTK_QUAD ) == test_case.nbQuad, f"Number of quads must be {test_case.nbQuad}" + assert counts.getTypeCount( VTK_TETRA ) == test_case.nbTetra, f"Number of tetrahedra must be {test_case.nbTetra}" + assert counts.getTypeCount( VTK_PYRAMID ) == test_case.nbPyr, f"Number of pyramids must be {test_case.nbPyr}" + assert counts.getTypeCount( VTK_WEDGE ) == test_case.nbWed, f"Number of wedges must be {test_case.nbWed}" + assert counts.getTypeCount( VTK_HEXAHEDRON ) == test_case.nbHexa, f"Number of hexahedra must be {test_case.nbHexa}" + + +@pytest.mark.parametrize( "test_case", __generate_test_data() ) +def test_CellTypeCounts_setCount( test_case: TestCase ) -> None: + """Test of CellTypeCounts . + + Args: + test_case (TestCase): test case + """ + counts: CellTypeCounts = CellTypeCounts() + counts.setTypeCount( VTK_VERTEX, test_case.nbVertex ) + counts.setTypeCount( VTK_TRIANGLE, test_case.nbTri ) + counts.setTypeCount( VTK_QUAD, test_case.nbQuad ) + counts.setTypeCount( VTK_TETRA, test_case.nbTetra ) + counts.setTypeCount( VTK_PYRAMID, test_case.nbPyr ) + counts.setTypeCount( VTK_WEDGE, test_case.nbWed ) + counts.setTypeCount( VTK_HEXAHEDRON, test_case.nbHexa ) + + assert counts.getTypeCount( VTK_VERTEX ) == test_case.nbVertex, f"Number of vertices must be {test_case.nbVertex}" + assert counts.getTypeCount( VTK_TRIANGLE ) == test_case.nbTri, f"Number of triangles must be {test_case.nbTri}" + assert counts.getTypeCount( VTK_QUAD ) == test_case.nbQuad, f"Number of quads must be {test_case.nbQuad}" + assert counts.getTypeCount( VTK_TETRA ) == test_case.nbTetra, f"Number of tetrahedra must be {test_case.nbTetra}" + assert counts.getTypeCount( VTK_PYRAMID ) == test_case.nbPyr, f"Number of pyramids must be {test_case.nbPyr}" + assert counts.getTypeCount( VTK_WEDGE ) == test_case.nbWed, f"Number of wedges must be {test_case.nbWed}" + assert counts.getTypeCount( VTK_HEXAHEDRON ) == test_case.nbHexa, f"Number of hexahedra must be {test_case.nbHexa}" + + +@pytest.mark.parametrize( "test_case", __generate_test_data() ) +def test_CellTypeCounts_add( test_case: TestCase ) -> None: + """Test of CellTypeCounts . + + Args: + test_case (TestCase): test case + """ + counts1: CellTypeCounts = CellTypeCounts() + counts1.setTypeCount( VTK_VERTEX, test_case.nbVertex ) + counts1.setTypeCount( VTK_TRIANGLE, test_case.nbTri ) + counts1.setTypeCount( VTK_QUAD, test_case.nbQuad ) + counts1.setTypeCount( VTK_TETRA, test_case.nbTetra ) + counts1.setTypeCount( VTK_PYRAMID, test_case.nbPyr ) + counts1.setTypeCount( VTK_WEDGE, test_case.nbWed ) + counts1.setTypeCount( VTK_HEXAHEDRON, test_case.nbHexa ) + + counts2: CellTypeCounts = CellTypeCounts() + counts2.setTypeCount( VTK_VERTEX, test_case.nbVertex ) + counts2.setTypeCount( VTK_TRIANGLE, test_case.nbTri ) + counts2.setTypeCount( VTK_QUAD, test_case.nbQuad ) + counts2.setTypeCount( VTK_TETRA, test_case.nbTetra ) + counts2.setTypeCount( VTK_PYRAMID, test_case.nbPyr ) + counts2.setTypeCount( VTK_WEDGE, test_case.nbWed ) + counts2.setTypeCount( VTK_HEXAHEDRON, test_case.nbHexa ) + + newcounts: CellTypeCounts = counts1 + counts2 + assert newcounts.getTypeCount( VTK_VERTEX ) == int( + 2 * test_case.nbVertex ), f"Number of vertices must be {int(2 * test_case.nbVertex)}" + assert newcounts.getTypeCount( VTK_TRIANGLE ) == int( + 2 * test_case.nbTri ), f"Number of triangles must be {int(2 * test_case.nbTri)}" + assert newcounts.getTypeCount( VTK_QUAD ) == int( + 2 * test_case.nbQuad ), f"Number of quads must be {int(2 * test_case.nbQuad)}" + assert newcounts.getTypeCount( VTK_TETRA ) == int( + 2 * test_case.nbTetra ), f"Number of tetrahedra must be {int(2 * test_case.nbTetra)}" + assert newcounts.getTypeCount( VTK_PYRAMID ) == int( + 2 * test_case.nbPyr ), f"Number of pyramids must be {int(2 * test_case.nbPyr)}" + assert newcounts.getTypeCount( VTK_WEDGE ) == int( + 2 * test_case.nbWed ), f"Number of wedges must be {int(2 * test_case.nbWed)}" + assert newcounts.getTypeCount( VTK_HEXAHEDRON ) == int( + 2 * test_case.nbHexa ), f"Number of hexahedra must be {int(2 * test_case.nbHexa)}" + + +@pytest.mark.parametrize( "test_case", __generate_test_data() ) +def test_CellTypeCounts_print( test_case: TestCase ) -> None: + """Test of CellTypeCounts . + + Args: + test_case (TestCase): test case + """ + counts: CellTypeCounts = CellTypeCounts() + counts.setTypeCount( VTK_VERTEX, test_case.nbVertex ) + counts.setTypeCount( VTK_TRIANGLE, test_case.nbTri ) + counts.setTypeCount( VTK_QUAD, test_case.nbQuad ) + counts.setTypeCount( VTK_TETRA, test_case.nbTetra ) + counts.setTypeCount( VTK_PYRAMID, test_case.nbPyr ) + counts.setTypeCount( VTK_WEDGE, test_case.nbWed ) + counts.setTypeCount( VTK_HEXAHEDRON, test_case.nbHexa ) + line: str = counts.print() + lineExp: str = __get_expected_counts( test_case.nbVertex, test_case.nbTri, test_case.nbQuad, test_case.nbTetra, + test_case.nbPyr, test_case.nbWed, test_case.nbHexa ) + assert line == lineExp, "Output counts string differs from expected value." diff --git a/geos-mesh/tests/test_SplitMesh.py b/geos-mesh/tests/test_SplitMesh.py new file mode 100644 index 00000000..2a73cb58 --- /dev/null +++ b/geos-mesh/tests/test_SplitMesh.py @@ -0,0 +1,205 @@ +# SPDX-FileContributor: Martin Lemay +# SPDX-License-Identifier: Apache 2.0 +# ruff: noqa: E402 # disable Module level import not at top of file +import os +from dataclasses import dataclass +import numpy as np +import numpy.typing as npt +import pytest +from typing import ( + Iterator, ) + +from geos.mesh.utils.genericHelpers import createSingleCellMesh +from geos.mesh.processing.SplitMesh import SplitMesh + +from vtkmodules.util.numpy_support import vtk_to_numpy + +from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkCellArray, vtkCellData, vtkCellTypes, VTK_TRIANGLE, + VTK_QUAD, VTK_TETRA, VTK_HEXAHEDRON, VTK_PYRAMID ) + +from vtkmodules.vtkCommonCore import ( + vtkPoints, + vtkIdList, + vtkDataArray, +) + +data_root: str = os.path.join( os.path.dirname( os.path.abspath( __file__ ) ), "data" ) + +############################################################### +# create single tetra mesh # +############################################################### +tetra_cell_type: int = VTK_TETRA +tetra_path = "tetra_cell.csv" + +# expected results +tetra_points_out: npt.NDArray[ np.float64 ] = np.array( + [ [ 0., 0., 0. ], [ 1., 0., 0. ], [ 0., 0., 1. ], [ 0., 1., 0. ], [ 0.5, 0., 0. ], [ 0.5, 0., 0.5 ], + [ 0., 0., 0.5 ], [ 0., 0.5, 0. ], [ 0., 0.5, 0.5 ], [ 0.5, 0.5, 0. ] ], np.float64 ) +tetra_cells_out: list[ list[ int ] ] = [ [ 0, 4, 6, 7 ], [ 7, 9, 8, 3 ], [ 9, 4, 5, 1 ], [ 5, 6, 8, 2 ], [ 6, 8, 7, 4 ], + [ 4, 8, 7, 9 ], [ 4, 8, 9, 5 ], [ 5, 4, 8, 6 ] ] + +############################################################### +# create single hexa mesh # +############################################################### +hexa_cell_type: int = VTK_HEXAHEDRON +hexa_path = "hexa_cell.csv" + +# expected results +hexa_points_out: npt.NDArray[ np.float64 ] = np.array( + [ [ 0.0, 0.0, 0.0 ], [ 1.0, 0.0, 0.0 ], [ 1.0, 1.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 1.0 ], [ 1.0, 0.0, 1.0 ], + [ 1.0, 1.0, 1.0 ], [ 0.0, 1.0, 1.0 ], [ 0.5, 0.0, 0.0 ], [ 0.0, 0.5, 0.0 ], [ 0.0, 0.0, 0.5 ], [ 1.0, 0.5, 0.0 ], + [ 1.0, 0.0, 0.5 ], [ 0.5, 1.0, 0.0 ], [ 1.0, 1.0, 0.5 ], [ 0.0, 1.0, 0.5 ], [ 0.5, 0.0, 1.0 ], [ 0.0, 0.5, 1.0 ], + [ 1.0, 0.5, 1.0 ], [ 0.5, 1.0, 1.0 ], [ 0.5, 0.5, 0.0 ], [ 0.5, 0.0, 0.5 ], [ 0.0, 0.5, 0.5 ], [ 1.0, 0.5, 0.5 ], + [ 0.5, 1.0, 0.5 ], [ 0.5, 0.5, 1.0 ], [ 0.5, 0.5, 0.5 ] ], np.float64 ) +hexa_cells_out: list[ list[ int ] ] = [ [ 10, 21, 26, 22, 4, 16, 25, 17 ], [ 21, 12, 23, 26, 16, 5, 18, 25 ], + [ 0, 8, 20, 9, 10, 21, 26, 22 ], [ 8, 1, 11, 20, 21, 12, 23, 26 ], + [ 22, 26, 24, 15, 17, 25, 19, 7 ], [ 26, 23, 14, 24, 25, 18, 6, 19 ], + [ 9, 20, 13, 3, 22, 26, 24, 15 ], [ 20, 11, 2, 13, 26, 23, 14, 24 ] ] + +############################################################### +# create single pyramid mesh # +############################################################### +pyramid_cell_type: int = VTK_PYRAMID +pyramid_path = "pyramid_cell.csv" + +# expected results +pyramid_points_out: npt.NDArray[ np.float64 ] = np.array( + [ [ 0.0, 0.0, 0.0 ], [ 1.0, 0.0, 0.0 ], [ 1.0, 1.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.5, 0.5, 1.0 ], [ 0.5, 0.0, 0.0 ], + [ 0.0, 0.5, 0.0 ], [ 0.25, 0.25, 0.5 ], [ 1.0, 0.5, 0.0 ], [ 0.75, 0.25, 0.5 ], [ 0.5, 1.0, 0.0 ], + [ 0.75, 0.75, 0.5 ], [ 0.25, 0.75, 0.5 ], [ 0.5, 0.5, 0.0 ] ], np.float64 ) +pyramid_cells_out: list[ list[ int ] ] = [ [ 5, 1, 8, 13, 9 ], [ 13, 8, 2, 10, 11 ], [ 3, 6, 13, 10, 12 ], + [ 6, 0, 5, 13, 7 ], [ 12, 7, 9, 11, 4 ], [ 11, 9, 7, 12, 13 ], + [ 7, 9, 5, 13 ], [ 9, 11, 8, 13 ], [ 11, 12, 10, 13 ], [ 12, 7, 6, 13 ] ] + +############################################################### +# create single triangle mesh # +############################################################### +triangle_cell_type: int = VTK_TRIANGLE +triangle_path = "triangle_cell.csv" + +# expected results +triangle_points_out: npt.NDArray[ np.float64 ] = np.array( [ [ 0.0, 0.0, 0.0 ], [ 1.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0 ], + [ 0.5, 0.0, 0.0 ], [ 0.5, 0.5, 0.0 ], [ 0.0, 0.5, 0.0 ] ], + np.float64 ) +triangle_cells_out: list[ list[ int ] ] = [ [ 0, 3, 5 ], [ 3, 1, 4 ], [ 5, 4, 2 ], [ 3, 4, 5 ] ] + +############################################################### +# create single quad mesh # +############################################################### +quad_cell_type: int = VTK_QUAD +quad_path = "quad_cell.csv" + +# expected results +quad_points_out: npt.NDArray[ np.float64 ] = np.array( + [ [ 0.0, 0.0, 0.0 ], [ 1.0, 0.0, 0.0 ], [ 1.0, 1.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.5, 0.0, 0.0 ], [ 1.0, 0.5, 0.0 ], + [ 0.5, 1.0, 0.0 ], [ 0.0, 0.5, 0.0 ], [ 0.5, 0.5, 0.0 ] ], np.float64 ) +quad_cells_out: list[ list[ int ] ] = [ [ 0, 4, 8, 7 ], [ 4, 1, 5, 8 ], [ 8, 5, 2, 6 ], [ 7, 8, 6, 3 ] ] + +############################################################### +# create multi cell mesh # +############################################################### +# TODO: add tests cases composed of multi-cell meshes of various types + +data_filename_all = ( tetra_path, hexa_path, pyramid_path, triangle_path, quad_path ) +cell_types_all = ( tetra_cell_type, hexa_cell_type, pyramid_cell_type, triangle_cell_type, quad_cell_type ) +points_out_all = ( tetra_points_out, hexa_points_out, pyramid_points_out, triangle_points_out, quad_points_out ) +cells_out_all = ( tetra_cells_out, hexa_cells_out, pyramid_cells_out, triangle_cells_out, quad_cells_out ) + + +@dataclass( frozen=True ) +class TestCase: + """Test case.""" + __test__ = False + #: VTK cell type + cellType: int + #: mesh + mesh: vtkUnstructuredGrid + #: expected new point coordinates + pointsExp: npt.NDArray[ np.float64 ] + #: expected new cell point ids + cellsExp: list[ int ] + + +def __generate_split_mesh_test_data() -> Iterator[ TestCase ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: iterator on test cases + """ + for cellType, data_path, pointsExp, cellsExp in zip( cell_types_all, + data_filename_all, + points_out_all, + cells_out_all, + strict=True ): + ptsCoord: npt.NDArray[ np.float64 ] = np.loadtxt( os.path.join( data_root, data_path ), + dtype=float, + delimiter=',' ) + mesh: vtkUnstructuredGrid = createSingleCellMesh( cellType, ptsCoord ) + yield TestCase( cellType, mesh, pointsExp, cellsExp ) + + +ids = [ vtkCellTypes.GetClassNameFromTypeId( cellType ) for cellType in cell_types_all ] + + +@pytest.mark.parametrize( "test_case", __generate_split_mesh_test_data(), ids=ids ) +def test_single_cell_split( test_case: TestCase ) -> None: + """Test of SplitMesh filter with meshes composed of a single cell. + + Args: + test_case (TestCase): test case + """ + cellTypeName: str = vtkCellTypes.GetClassNameFromTypeId( test_case.cellType ) + filter: SplitMesh = SplitMesh() + filter.SetInputDataObject( test_case.mesh ) + filter.Update() + output: vtkUnstructuredGrid = filter.GetOutputDataObject( 0 ) + assert output is not None, "Output mesh is undefined." + pointsOut: vtkPoints = output.GetPoints() + assert pointsOut is not None, "Points from output mesh are undefined." + assert pointsOut.GetNumberOfPoints( + ) == test_case.pointsExp.shape[ 0 ], f"Number of points is expected to be {test_case.pointsExp.shape[0]}." + pointCoords: npt.NDArray[ np.float64 ] = vtk_to_numpy( pointsOut.GetData() ) + print( "Points coords: ", cellTypeName, pointCoords.tolist() ) + assert np.array_equal( pointCoords.ravel(), test_case.pointsExp.ravel() ), "Points coordinates mesh are wrong." + + cellsOut: vtkCellArray = output.GetCells() + typesArray0: npt.NDArray[ np.int64 ] = vtk_to_numpy( output.GetDistinctCellTypesArray() ) + print( "typesArray0", cellTypeName, typesArray0 ) + + assert cellsOut is not None, "Cells from output mesh are undefined." + assert cellsOut.GetNumberOfCells() == len( + test_case.cellsExp ), f"Number of cells is expected to be {len(test_case.cellsExp)}." + # check cell types + types: vtkCellTypes = vtkCellTypes() + output.GetCellTypes( types ) + assert types is not None, "Cell types must be defined" + typesArray: npt.NDArray[ np.int64 ] = vtk_to_numpy( types.GetCellTypesArray() ) + + print( "typesArray", cellTypeName, typesArray ) + assert ( typesArray.size == 1 ) and ( typesArray[ 0 ] == test_case.cellType ), f"All cells must be {cellTypeName}" + + for i in range( cellsOut.GetNumberOfCells() ): + ptIds = vtkIdList() + cellsOut.GetCellAtId( i, ptIds ) + cellsOutObs: list[ int ] = [ ptIds.GetId( j ) for j in range( ptIds.GetNumberOfIds() ) ] + nbPtsExp: int = len( test_case.cellsExp[ i ] ) + print( "cell type", cellTypeName, i, vtkCellTypes.GetClassNameFromTypeId( types.GetCellType( i ) ) ) + print( "cellsOutObs: ", cellTypeName, i, cellsOutObs ) + assert ptIds is not None, "Point ids must be defined" + assert ptIds.GetNumberOfIds() == nbPtsExp, f"Cells must be defined by {nbPtsExp} points." + assert cellsOutObs == test_case.cellsExp[ i ], "Cell point ids are wrong." + + # test originalId array was created + cellData: vtkCellData = output.GetCellData() + assert cellData is not None, "Cell data should be defined." + array: vtkDataArray = cellData.GetArray( "OriginalID" ) + assert array is not None, "OriginalID array should be defined." + + # test other arrays were transferred + cellDataInput: vtkCellData = test_case.mesh.GetCellData() + assert cellDataInput is not None, "Cell data from input mesh should be defined." + nbArrayInput: int = cellDataInput.GetNumberOfArrays() + nbArraySplited: int = cellData.GetNumberOfArrays() + assert nbArraySplited == nbArrayInput + 1, f"Number of arrays should be {nbArrayInput + 1}" + + #assert False diff --git a/geos-mesh/tests/test_collocated_nodes.py b/geos-mesh/tests/test_collocated_nodes.py index 2b74e30f..ecc36dbb 100644 --- a/geos-mesh/tests/test_collocated_nodes.py +++ b/geos-mesh/tests/test_collocated_nodes.py @@ -6,8 +6,7 @@ def get_points() -> Iterator[ Tuple[ vtkPoints, int ] ]: - """ - Generates the data for the cases. + """Generates the data for the cases. One case has two nodes at the exact same position. The other has two differente nodes :return: Generator to (vtk points, number of expected duplicated locations) diff --git a/geos-mesh/tests/test_genericHelpers.py b/geos-mesh/tests/test_genericHelpers.py new file mode 100644 index 00000000..f77bc6f2 --- /dev/null +++ b/geos-mesh/tests/test_genericHelpers.py @@ -0,0 +1,188 @@ +# SPDX-FileContributor: Martin Lemay +# SPDX-License-Identifier: Apache 2.0 +# ruff: noqa: E402 # disable Module level import not at top of file +import os +from dataclasses import dataclass +import numpy as np +import numpy.typing as npt +import pytest +from typing import ( + Iterator, ) + +from geos.mesh.utils.genericHelpers import getBoundsFromPointCoords, createVertices, createMultiCellMesh + +from vtkmodules.util.numpy_support import vtk_to_numpy + +from vtkmodules.vtkCommonDataModel import ( + vtkUnstructuredGrid, + vtkCellArray, + vtkCellTypes, + VTK_TETRA, + VTK_HEXAHEDRON, +) + +from vtkmodules.vtkCommonCore import ( + vtkPoints, + vtkIdList, +) + +# TODO: add case whith various cell types + +# inputs +data_root: str = os.path.join( os.path.dirname( os.path.abspath( __file__ ) ), "data" ) +data_filename_all: tuple[ str, ...] = ( "tetra_mesh.csv", "hexa_mesh.csv" ) +celltypes_all: tuple[ int ] = ( VTK_TETRA, VTK_HEXAHEDRON ) +nbPtsCell_all: tuple[ int ] = ( 4, 8 ) + +# expected results if shared vertices +hexa_points_out: npt.NDArray[ np.float64 ] = np.array( + [ [ 0.0, 0.0, 0.5 ], [ 0.5, 0.0, 0.5 ], [ 0.5, 0.5, 0.5 ], [ 0.0, 0.5, 0.5 ], [ 0.0, 0.0, 1.0 ], [ 0.5, 0.0, 1.0 ], + [ 0.5, 0.5, 1.0 ], [ 0.0, 0.5, 1.0 ], [ 1.0, 0.0, 0.5 ], [ 1.0, 0.5, 0.5 ], [ 1.0, 0.0, 1.0 ], [ 1.0, 0.5, 1.0 ], + [ 0.0, 0.0, 0.0 ], [ 0.5, 0.0, 0.0 ], [ 0.5, 0.5, 0.0 ], [ 0.0, 0.5, 0.0 ], [ 1.0, 0.0, 0.0 ], [ 1.0, 0.5, 0.0 ], + [ 0.5, 1.0, 0.5 ], [ 0.0, 1.0, 0.5 ], [ 0.5, 1.0, 1.0 ], [ 0.0, 1.0, 1.0 ], [ 1.0, 1.0, 0.5 ], [ 1.0, 1.0, 1.0 ], + [ 0.5, 1.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 1.0, 1.0, 0.0 ] ], np.float64 ) +tetra_points_out: npt.NDArray[ np.float64 ] = np.array( + [ [ 0.0, 0.0, 0.0 ], [ 0.5, 0.0, 0.0 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.5, 0.0 ], [ 0.5, 0.5, 0.0 ], [ 0.0, 0.5, 0.5 ], + [ 0.0, 1.0, 0.0 ], [ 0.5, 0.0, 0.5 ], [ 1.0, 0.0, 0.0 ], [ 0.0, 0.0, 1.0 ] ], np.float64 ) +points_out_all = ( tetra_points_out, hexa_points_out ) + +tetra_cellPtsIdsExp = [ ( 0, 1, 2, 3 ), ( 3, 4, 5, 6 ), ( 4, 1, 7, 8 ), ( 7, 2, 5, 9 ), ( 2, 5, 3, 1 ), ( 1, 5, 3, 4 ), + ( 1, 5, 4, 7 ), ( 7, 1, 5, 2 ) ] +hexa_cellPtsIdsExp = [ ( 0, 1, 2, 3, 4, 5, 6, 7 ), ( 1, 8, 9, 2, 5, 10, 11, 6 ), ( 12, 13, 14, 15, 0, 1, 2, 3 ), + ( 13, 16, 17, 14, 1, 8, 9, 2 ), ( 3, 2, 18, 19, 7, 6, 20, 21 ), ( 2, 9, 22, 18, 6, 11, 23, 20 ), + ( 15, 14, 24, 25, 3, 2, 18, 19 ), ( 14, 17, 26, 24, 2, 9, 22, 18 ) ] +cellPtsIdsExp_all = ( tetra_cellPtsIdsExp, hexa_cellPtsIdsExp ) + + +@dataclass( frozen=True ) +class TestCase: + """Test case.""" + __test__ = False + #: cell types + cellTypes: list[ int ] + #: cell point coordinates + cellPtsCoords: list[ npt.NDArray[ np.float64 ] ] + #: share or unshare vertices + share: bool + #: expected points + pointsExp: npt.NDArray[ np.float64 ] + #: expected cell point ids + cellPtsIdsExp: tuple[ tuple[ int ] ] + + +def __generate_test_data() -> Iterator[ TestCase ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: iterator on test cases + """ + for path, celltype, nbPtsCell, pointsExp0, cellPtsIdsExp0 in zip( data_filename_all, + celltypes_all, + nbPtsCell_all, + points_out_all, + cellPtsIdsExp_all, + strict=True ): + # all points coordinates + ptsCoords: npt.NDArray[ np.float64 ] = np.loadtxt( os.path.join( data_root, path ), dtype=float, delimiter=',' ) + # split array to get a list of coordinates per cell + cellPtsCoords: list[ npt.NDArray[ np.float64 ] ] = [ + ptsCoords[ i:i + nbPtsCell ] for i in range( 0, ptsCoords.shape[ 0 ], nbPtsCell ) + ] + nbCells: int = int( ptsCoords.shape[ 0 ] / nbPtsCell ) + cellTypes = nbCells * [ celltype ] + for shared in ( False, True ): + pointsExp: npt.NDArray[ np.float64 ] = pointsExp0 if shared else ptsCoords + cellPtsIdsExp = cellPtsIdsExp0 if shared else [ + tuple( range( i * nbPtsCell, ( i + 1 ) * nbPtsCell, 1 ) ) for i in range( nbCells ) + ] + yield TestCase( cellTypes, cellPtsCoords, shared, pointsExp, cellPtsIdsExp ) + + +ids: list[ str ] = [ + os.path.splitext( name )[ 0 ] + f"_{shared}]" for name in data_filename_all for shared in ( False, True ) +] + + +@pytest.mark.parametrize( "test_case", __generate_test_data(), ids=ids ) +def test_createVertices( test_case: TestCase ) -> None: + """Test of createVertices method. + + Args: + test_case (TestCase): test case + """ + pointsOut: vtkPoints + cellPtsIds: list[ tuple[ int, ...] ] + pointsOut, cellPtsIds = createVertices( test_case.cellPtsCoords, test_case.share ) + assert pointsOut is not None, "Output points is undefined." + assert cellPtsIds is not None, "Output cell point map is undefined." + nbPtsExp: int = test_case.pointsExp.shape[ 0 ] + assert pointsOut.GetNumberOfPoints() == nbPtsExp, f"Number of points is expected to be {nbPtsExp}." + pointCoords: npt.NDArray[ np.float64 ] = vtk_to_numpy( pointsOut.GetData() ) + print( "Points coords Obs: ", pointCoords.tolist() ) + assert np.array_equal( pointCoords, test_case.pointsExp ), "Points coordinates are wrong." + print( "Cell points coords: ", cellPtsIds ) + assert cellPtsIds == test_case.cellPtsIdsExp, f"Cell point Ids are expected to be {test_case.cellPtsIdsExp}" + + +ids: list[ str ] = [ + os.path.splitext( name )[ 0 ] + f"_{shared}]" for name in data_filename_all for shared in ( False, True ) +] + + +@pytest.mark.parametrize( "test_case", __generate_test_data(), ids=ids ) +def test_createMultiCellMesh( test_case: TestCase ) -> None: + """Test of createMultiCellMesh method. + + Args: + test_case (TestCase): test case + """ + output: vtkUnstructuredGrid = createMultiCellMesh( test_case.cellTypes, test_case.cellPtsCoords, test_case.share ) + assert output is not None, "Output mesh is undefined." + + # tests on points + pointsOut: vtkPoints = output.GetPoints() + assert pointsOut is not None, "Output points is undefined." + nbPtsExp: int = test_case.pointsExp.shape[ 0 ] + assert pointsOut.GetNumberOfPoints() == nbPtsExp, f"Number of points is expected to be {nbPtsExp}." + pointCoords: npt.NDArray[ np.float64 ] = vtk_to_numpy( pointsOut.GetData() ) + print( "Points coords Obs: ", pointCoords.tolist() ) + assert np.array_equal( pointCoords, test_case.pointsExp ), "Points coordinates are wrong." + + # tests on cells + cellsOut: vtkCellArray = output.GetCells() + assert cellsOut is not None, "Cells from output mesh are undefined." + nbCells: int = len( test_case.cellPtsCoords ) + assert cellsOut.GetNumberOfCells() == nbCells, f"Number of cells is expected to be {nbCells}." + + # check cell types + types: vtkCellTypes = vtkCellTypes() + output.GetCellTypes( types ) + assert types is not None, "Cell types must be defined" + typesArray: npt.NDArray[ np.int64 ] = vtk_to_numpy( types.GetCellTypesArray() ) + print( "typesArray.size ", typesArray.size ) + assert ( typesArray.size == 1 ) and ( typesArray[ 0 ] == test_case.cellTypes[ 0 ] ), "Cell types are wrong" + + for cellId in range( output.GetNumberOfCells() ): + ptIds = vtkIdList() + cellsOut.GetCellAtId( cellId, ptIds ) + cellsOutObs: tuple[ int ] = tuple( [ ptIds.GetId( j ) for j in range( ptIds.GetNumberOfIds() ) ] ) + print( "cellsOutObs: ", cellsOutObs ) + nbCellPts: int = len( test_case.cellPtsIdsExp[ cellId ] ) + assert ptIds is not None, "Point ids must be defined" + assert ptIds.GetNumberOfIds() == nbCellPts, f"Cells must be defined by {nbCellPts} points." + assert cellsOutObs == test_case.cellPtsIdsExp[ cellId ], "Cell point ids are wrong." + + +def test_getBoundsFromPointCoords() -> None: + """Test of getBoundsFromPointCoords method.""" + # input + cellPtsCoord: list[ npt.NDArray[ np.float64 ] ] = [ + np.array( [ [ 5, 4, 3 ], [ 1, 8, 4 ], [ 2, 5, 7 ] ], dtype=float ), + np.array( [ [ 1, 4, 6 ], [ 2, 7, 9 ], [ 4, 5, 6 ] ], dtype=float ), + np.array( [ [ 3, 7, 8 ], [ 5, 7, 3 ], [ 4, 7, 3 ] ], dtype=float ), + np.array( [ [ 1, 7, 2 ], [ 0, 1, 2 ], [ 2, 3, 7 ] ], dtype=float ), + ] + # expected output + boundsExp: list[ float ] = [ 0., 5., 1., 8., 2., 9. ] + boundsObs: list[ float ] = getBoundsFromPointCoords( cellPtsCoord ) + assert boundsExp == boundsObs, f"Expected bounds are {boundsExp}." diff --git a/geos-mesh/tests/test_genericHelpers_createSingleCellMesh.py b/geos-mesh/tests/test_genericHelpers_createSingleCellMesh.py new file mode 100644 index 00000000..d058c537 --- /dev/null +++ b/geos-mesh/tests/test_genericHelpers_createSingleCellMesh.py @@ -0,0 +1,98 @@ +# SPDX-FileContributor: Martin Lemay +# SPDX-License-Identifier: Apache 2.0 +# ruff: noqa: E402 # disable Module level import not at top of file +import os +from dataclasses import dataclass +import numpy as np +import numpy.typing as npt +import pytest +from typing import ( + Iterator, ) + +from geos.mesh.utils.genericHelpers import createSingleCellMesh + +from vtkmodules.util.numpy_support import vtk_to_numpy + +from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkCellArray, vtkCellTypes, VTK_TRIANGLE, VTK_QUAD, + VTK_TETRA, VTK_HEXAHEDRON, VTK_PYRAMID ) + +from vtkmodules.vtkCommonCore import ( + vtkPoints, + vtkIdList, +) + +# inputs +data_root: str = os.path.join( os.path.dirname( os.path.abspath( __file__ ) ), "data" ) + +data_filename_all: tuple[ str, ...] = ( "triangle_cell.csv", "quad_cell.csv", "tetra_cell.csv", "pyramid_cell.csv", + "hexa_cell.csv" ) +cell_type_all: tuple[ int, ...] = ( VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_HEXAHEDRON ) + + +@dataclass( frozen=True ) +class TestCase: + """Test case.""" + __test__ = False + #: VTK cell type + cellType: int + #: cell point coordinates + cellPoints: npt.NDArray[ np.float64 ] + + +def __generate_test_data() -> Iterator[ TestCase ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: iterator on test cases + """ + for cellType, path in zip( cell_type_all, data_filename_all, strict=True ): + cell: npt.NDArray[ np.float64 ] = np.loadtxt( os.path.join( data_root, path ), dtype=float, delimiter=',' ) + yield TestCase( cellType, cell ) + + +ids: list[ str ] = [ vtkCellTypes.GetClassNameFromTypeId( cellType ) for cellType in cell_type_all ] + + +@pytest.mark.parametrize( "test_case", __generate_test_data(), ids=ids ) +def test_createSingleCellMesh( test_case: TestCase ) -> None: + """Test of createSingleCellMesh method. + + Args: + test_case (TestCase): test case + """ + cellTypeName: str = vtkCellTypes.GetClassNameFromTypeId( test_case.cellType ) + output: vtkUnstructuredGrid = createSingleCellMesh( test_case.cellType, test_case.cellPoints ) + + assert output is not None, "Output mesh is undefined." + pointsOut: vtkPoints = output.GetPoints() + nbPtsExp: int = len( test_case.cellPoints ) + assert pointsOut is not None, "Points from output mesh are undefined." + assert pointsOut.GetNumberOfPoints() == nbPtsExp, f"Number of points is expected to be {nbPtsExp}." + pointCoords: npt.NDArray[ np.float64 ] = vtk_to_numpy( pointsOut.GetData() ) + print( "Points coords: ", cellTypeName, pointCoords.tolist() ) + assert np.array_equal( pointCoords.ravel(), test_case.cellPoints.ravel() ), "Points coordinates are wrong." + + cellsOut: vtkCellArray = output.GetCells() + typesArray0: npt.NDArray[ np.int64 ] = vtk_to_numpy( output.GetDistinctCellTypesArray() ) + print( "typesArray0", cellTypeName, typesArray0 ) + + assert cellsOut is not None, "Cells from output mesh are undefined." + assert cellsOut.GetNumberOfCells() == 1, "Number of cells is expected to be 1." + # check cell types + types: vtkCellTypes = vtkCellTypes() + output.GetCellTypes( types ) + assert types is not None, "Cell types must be defined" + typesArray: npt.NDArray[ np.int64 ] = vtk_to_numpy( types.GetCellTypesArray() ) + + print( "typesArray", cellTypeName, typesArray ) + assert ( typesArray.size == 1 ) and ( typesArray[ 0 ] == test_case.cellType ), f"Cell must be {cellTypeName}" + + ptIds = vtkIdList() + cellsOut.GetCellAtId( 0, ptIds ) + cellsOutObs: list[ int ] = [ ptIds.GetId( j ) for j in range( ptIds.GetNumberOfIds() ) ] + + print( "cell type", cellTypeName, vtkCellTypes.GetClassNameFromTypeId( types.GetCellType( 0 ) ) ) + print( "cellsOutObs: ", cellTypeName, cellsOutObs ) + assert ptIds is not None, "Point ids must be defined" + assert ptIds.GetNumberOfIds() == nbPtsExp, f"Cells must be defined by {nbPtsExp} points." + assert cellsOutObs == list( range( nbPtsExp ) ), "Cell point ids are wrong." diff --git a/geos-mesh/tests/test_supported_elements.py b/geos-mesh/tests/test_supported_elements.py index bb0213c6..5e87fb38 100644 --- a/geos-mesh/tests/test_supported_elements.py +++ b/geos-mesh/tests/test_supported_elements.py @@ -11,8 +11,7 @@ # TODO Update this test to have access to another meshTests file @pytest.mark.parametrize( "base_name", ( "supportedElements.vtk", "supportedElementsAsVTKPolyhedra.vtk" ) ) def test_supported_elements( base_name ) -> None: - """ - Testing that the supported elements are properly detected as supported! + """Testing that the supported elements are properly detected as supported! :param base_name: Supported elements are provided as standard elements or polyhedron elements. """ ... @@ -25,8 +24,7 @@ def test_supported_elements( base_name ) -> None: def make_dodecahedron() -> Tuple[ vtkPoints, vtkIdList ]: - """ - Returns the points and faces for a dodecahedron. + """Returns the points and faces for a dodecahedron. This code was adapted from an official vtk example. :return: The tuple of points and faces (as vtk instances). """ @@ -81,8 +79,7 @@ def make_dodecahedron() -> Tuple[ vtkPoints, vtkIdList ]: # TODO make this test work def test_dodecahedron() -> None: - """ - Tests whether a dodecahedron is support by GEOS or not. + """Tests whether a dodecahedron is support by GEOS or not. """ points, faces = make_dodecahedron() mesh = vtkUnstructuredGrid() diff --git a/geos-pv/src/PVplugins/PVCellTypeCounter.py b/geos-pv/src/PVplugins/PVCellTypeCounter.py new file mode 100644 index 00000000..fb465d86 --- /dev/null +++ b/geos-pv/src/PVplugins/PVCellTypeCounter.py @@ -0,0 +1,159 @@ +# 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 sys +from pathlib import Path +from typing_extensions import Self +from typing import Optional + +from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] + VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy ) + +from vtkmodules.vtkCommonCore import ( + vtkInformation, + vtkInformationVector, +) +from vtkmodules.vtkCommonDataModel import ( + vtkPointSet, + vtkTable, +) + +# update sys.path to load all GEOS Python Package dependencies +geos_pv_path: Path = Path( __file__ ).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.stats.CellTypeCounter import CellTypeCounter +from geos.mesh.model.CellTypeCounts import CellTypeCounts + +__doc__ = """ +Compute cell type counts. Counts are dumped in to Output message window and can be exporter in a file. + +To use it: + +* Load the module in Paraview: Tools>Manage Plugins...>Load new>PVCellTypeCounter. +* Select the input mesh. +* Apply the filter. + +""" + + +@smproxy.filter( name="PVCellTypeCounter", label="Cell Type Counter" ) +@smhint.xml( '' ) +@smproperty.input( name="Input", port_index=0 ) +@smdomain.datatype( + dataTypes=[ "vtkUnstructuredGrid" ], + composite_data_supported=True, +) +class PVCellTypeCounter( VTKPythonAlgorithmBase ): + + def __init__( self: Self ) -> None: + """Merge collocated points.""" + super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkTable" ) + + self._filename: Optional[ str ] = None + self._saveToFile: bool = True + # used to concatenate results if vtkMultiBlockDataSet + self._countsAll: CellTypeCounts = CellTypeCounts() + + @smproperty.intvector( + name="SetSaveToFile", + label="Save to file", + default_values=0, + panel_visibility="default", + ) + @smdomain.xml( """ + + + Specify if mesh statistics are dumped into a file. + + """ ) + def SetSaveToFile( self: Self, saveToFile: bool ) -> None: + """Setter to save the stats into a file. + + Args: + saveToFile (bool): if True, a file will be saved. + """ + if self._saveToFile != saveToFile: + self._saveToFile = saveToFile + self.Modified() + + @smproperty.stringvector( name="FilePath", label="File Path" ) + @smdomain.xml( """ + + Output file path. + + + + + """ ) + def SetFileName( self: Self, fname: str ) -> None: + """Specify filename for the filter to write. + + Args: + fname (str): file path + """ + if self._filename != fname: + self._filename = fname + self.Modified() + + @smproperty.xml( """ + + + + + + + """ ) + def d09GroupAdvancedOutputParameters( self: Self ) -> None: + """Organize groups.""" + self.Modified() + + def RequestData( + self: Self, + request: vtkInformation, # noqa: F841 + inInfoVec: list[ vtkInformationVector ], + outInfoVec: vtkInformationVector, + ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestData. + + Args: + request (vtkInformation): request + inInfoVec (list[vtkInformationVector]): input objects + outInfoVec (vtkInformationVector): output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + inputMesh: vtkPointSet = self.GetInputData( inInfoVec, 0, 0 ) + outputTable: vtkTable = vtkTable.GetData( outInfoVec, 0 ) + assert inputMesh is not None, "Input server mesh is null." + assert outputTable is not None, "Output pipeline is null." + + filter: CellTypeCounter = CellTypeCounter() + filter.SetInputDataObject( inputMesh ) + filter.Update() + outputTable.ShallowCopy( filter.GetOutputDataObject( 0 ) ) + + # print counts in Output Messages view + counts: CellTypeCounts = filter.GetCellTypeCounts() + print( counts.print() ) + + self._countsAll += counts + # save to file if asked + if self._saveToFile and self._filename is not None: + try: + with open( self._filename, 'w' ) as fout: + fout.write( self._countsAll.print() ) + print( f"File {self._filename} was successfully written." ) + except Exception as e: + print( "Error while exporting the file due to:" ) + print( str( e ) ) + return 1 diff --git a/geos-pv/src/PVplugins/PVSplitMesh.py b/geos-pv/src/PVplugins/PVSplitMesh.py new file mode 100644 index 00000000..673f790e --- /dev/null +++ b/geos-pv/src/PVplugins/PVSplitMesh.py @@ -0,0 +1,68 @@ +# 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 sys +from pathlib import Path +from typing_extensions import Self + +from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] + smdomain, smhint, smproperty, smproxy, +) + +from vtkmodules.vtkCommonDataModel import ( + vtkPointSet, ) + +# update sys.path to load all GEOS Python Package dependencies +geos_pv_path: Path = Path( __file__ ).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.SplitMesh import SplitMesh +from geos.pv.utils.AbstractPVPluginVtkWrapper import AbstractPVPluginVtkWrapper + +__doc__ = """ +Split each cell of input mesh to smaller cells. + +Output mesh is of same type as input mesh. If input mesh is a composite mesh, the plugin split cells of each part independently. + +To use it: + +* Load the module in Paraview: Tools>Manage Plugins...>Load new>PVSplitMesh. +* Select the input mesh. +* Apply the filter. + +""" + + +@smproxy.filter( name="PVSplitMesh", label="Split Mesh" ) +@smhint.xml( '' ) +@smproperty.input( name="Input", port_index=0 ) +@smdomain.datatype( + dataTypes=[ "vtkPointSet" ], + composite_data_supported=True, +) +class PVSplitMesh( AbstractPVPluginVtkWrapper ): + + def __init__( self: Self ) -> None: + """Split mesh cells.""" + super().__init__() + + def applyVtkFilter( + self: Self, + input: vtkPointSet, + ) -> vtkPointSet: + """Apply vtk filter. + + Args: + input (vtkPointSet): input mesh + + Returns: + vtkPointSet: output mesh + """ + filter: SplitMesh = SplitMesh() + filter.SetInputDataObject( input ) + filter.Update() + return filter.GetOutputDataObject( 0 ) diff --git a/geos-pv/src/geos/pv/utils/AbstractPVPluginVtkWrapper.py b/geos-pv/src/geos/pv/utils/AbstractPVPluginVtkWrapper.py new file mode 100644 index 00000000..8ae8c27f --- /dev/null +++ b/geos-pv/src/geos/pv/utils/AbstractPVPluginVtkWrapper.py @@ -0,0 +1,97 @@ +# 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 Any +from typing_extensions import Self + +from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] + VTKPythonAlgorithmBase, ) + +from vtkmodules.vtkCommonCore import ( + vtkInformation, + vtkInformationVector, +) + +__doc__ = """ +AbstractPVPluginVtkWrapper module defines the parent Paraview plugin from which inheritates PV plugins that directly wrap a vtk filter. + +To use it, make children PV plugins inherited from AbstractPVPluginVtkWrapper. Output mesh is of same type as input mesh. +If output type needs to be specified, this must be done in the child class. +""" + + +class AbstractPVPluginVtkWrapper( VTKPythonAlgorithmBase ): + + def __init__( self: Self ) -> None: + """Abstract Paraview Plugin class.""" + super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkPointSet" ) + + 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 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: + inputMesh: Any = self.GetInputData( inInfoVec, 0, 0 ) + outputMesh: Any = self.GetOutputData( outInfoVec, 0 ) + assert inputMesh is not None, "Input server mesh is null." + assert outputMesh is not None, "Output pipeline is null." + + tmpMesh = self.applyVtkFilter( inputMesh ) + assert tmpMesh is not None, "Output mesh is null." + outputMesh.ShallowCopy( tmpMesh ) + print( "Filter was successfully applied." ) + except ( AssertionError, Exception ) as e: + print( f"Filter failed due to: {e}" ) + return 0 + return 1 + + def applyVtkFilter( + self: Self, + input: Any, + ) -> Any: + """Apply vtk filter. + + Args: + input (Any): input object + + Returns: + Any: output mesh + """ + pass