diff --git a/README.md b/README.md index 4dc52d71..e2c8e11b 100644 --- a/README.md +++ b/README.md @@ -145,7 +145,14 @@ If you would like to contribute to GEOS Python packages, please respect the foll 1. Create a new branch named from this template: `[CONTRIBUTOR]/[TYPE]/[TITLE]` where CONTRIBUTOR is the name of the contributor, TYPE is the type of contribution among 'feature', 'refactor', 'doc', 'ci', TITLE is a short title for the branch. 2. Add your code trying to integrate into the current code architecture. -3. Push the branch, open a new PR respecting naming [semantics](https://gist.github.com/joshbuchea/6f47e86d2510bce28f8e7f42ae84c716), and add reviewers +3. Run mypy, ruff and yapf in this order + ``` + python -m pip install --upgrade mypy ruff yapf + python -m mypy --config-file ./.mypy.ini --check-untyped-defs ./ + python -m ruff check --fix --config .ruff.toml ./ + python -m yapf -r -i --style .style.yapf ./ + ``` +4. Push the branch, open a new PR respecting naming [semantics](https://gist.github.com/joshbuchea/6f47e86d2510bce28f8e7f42ae84c716), and add reviewers If you do not have the rights to push the code and open new PRs, consider opening a new issue to explain what you want to do and ask for the dev rights. diff --git a/docs/geos_mesh_docs/model.rst b/docs/geos_mesh_docs/model.rst index ead57332..49e07753 100644 --- a/docs/geos_mesh_docs/model.rst +++ b/docs/geos_mesh_docs/model.rst @@ -4,10 +4,18 @@ Model The `model` module of `geos-mesh` package contains data model. -geos.mesh.model.CellTypeCounts filter +geos.mesh.model.CellTypeCounts module -------------------------------------- .. automodule:: geos.mesh.model.CellTypeCounts + :members: + :undoc-members: + :show-inheritance: + +geos.mesh.model.QualityMetricSummary module +-------------------------------------------- + +.. automodule:: geos.mesh.model.QualityMetricSummary :members: :undoc-members: :show-inheritance: \ No newline at end of file diff --git a/docs/geos_mesh_docs/processing.rst b/docs/geos_mesh_docs/processing.rst index d79db6db..8a448e31 100644 --- a/docs/geos_mesh_docs/processing.rst +++ b/docs/geos_mesh_docs/processing.rst @@ -4,6 +4,14 @@ Processing filters The `processing` module of `geos-mesh` package contains filters to process meshes. +geos.mesh.processing.meshQualityMetricHelpers module +----------------------------------------------------- + +.. automodule:: geos.mesh.processing.meshQualityMetricHelpers + :members: + :undoc-members: + :show-inheritance: + geos.mesh.processing.SplitMesh filter -------------------------------------- diff --git a/docs/geos_mesh_docs/stats.rst b/docs/geos_mesh_docs/stats.rst index 50664514..77e9eee6 100644 --- a/docs/geos_mesh_docs/stats.rst +++ b/docs/geos_mesh_docs/stats.rst @@ -4,10 +4,18 @@ Mesh stats tools The `stats` module of `geos-mesh` package contains filter to compute statistics on meshes. -geos.mesh.stats.CellTypeCounter filter --------------------------------------- +geos.mesh.stats.CellTypeCounterEnhanced filter +----------------------------------------------- -.. automodule:: geos.mesh.stats.CellTypeCounter +.. automodule:: geos.mesh.stats.CellTypeCounterEnhanced + :members: + :undoc-members: + :show-inheritance: + +geos.mesh.stats.MeshQualityEnhanced filter +------------------------------------------- + +.. automodule:: geos.mesh.stats.MeshQualityEnhanced :members: :undoc-members: :show-inheritance: \ No newline at end of file diff --git a/geos-mesh/pyproject.toml b/geos-mesh/pyproject.toml index 8a994b47..3bad4c2b 100644 --- a/geos-mesh/pyproject.toml +++ b/geos-mesh/pyproject.toml @@ -32,6 +32,7 @@ dependencies = [ "pandas >= 2.2", "meshio >= 5.3", "typing_extensions >= 4.12", + "geos-utils @ file:./geos-utils", ] [project.scripts] diff --git a/geos-mesh/src/geos/mesh/doctor/checks/collocated_nodes.py b/geos-mesh/src/geos/mesh/doctor/checks/collocated_nodes.py index 74cbbe8c..ec753908 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/collocated_nodes.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/collocated_nodes.py @@ -20,7 +20,7 @@ class Result: def __check( mesh, options: Options ) -> Result: - points = mesh.GetPoints() + points: vtkPoints = mesh.GetPoints() locator = vtkIncrementalOctreePointLocator() locator.SetTolerance( options.tolerance ) diff --git a/geos-mesh/src/geos/mesh/model/CellTypeCounts.py b/geos-mesh/src/geos/mesh/model/CellTypeCounts.py index bab12085..ba9b7692 100644 --- a/geos-mesh/src/geos/mesh/model/CellTypeCounts.py +++ b/geos-mesh/src/geos/mesh/model/CellTypeCounts.py @@ -10,6 +10,21 @@ __doc__ = """ CellTypeCounts stores the number of elements of each type. + +To use the filter: + +.. code-block:: python + + from geos.mesh.model.CellTypeCounts import CellTypeCounts + + counts: CellTypeCounts = CellTypeCounts() + + # set data + counts.addType( cellType ) + counts.setTypeCount( cellType, count ) + + # get data + count: int counts.getTypeCount( cellType ) """ @@ -43,6 +58,14 @@ def __add__( self: Self, other: Self ) -> 'CellTypeCounts': newCounts._counts = self._counts + other._counts return newCounts + def getCounts( self: Self ) -> npt.NDArray[ np.int64 ]: + """Get all counts. + + Returns: + npt.NDArray[ np.int64 ]: counts + """ + return self._counts + def addType( self: Self, cellType: int ) -> None: """Increment the number of cell of input type. @@ -74,6 +97,10 @@ def getTypeCount( self: Self, cellType: int ) -> int: """ return int( self._counts[ cellType ] ) + def reset( self: Self ) -> None: + """Reset counts.""" + self._counts = np.zeros( VTK_NUMBER_OF_CELL_TYPES, dtype=float ) + def _updateGeneralCounts( self: Self, cellType: int, count: int ) -> None: """Update generic type counters. diff --git a/geos-mesh/src/geos/mesh/model/QualityMetricSummary.py b/geos-mesh/src/geos/mesh/model/QualityMetricSummary.py new file mode 100644 index 00000000..eff54e15 --- /dev/null +++ b/geos-mesh/src/geos/mesh/model/QualityMetricSummary.py @@ -0,0 +1,586 @@ +# 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 +import pandas as pd +from enum import Enum +from typing import Any +from typing_extensions import Self, Iterable +from packaging.version import Version +import matplotlib as mpl +import matplotlib.cm as cm +import matplotlib.pyplot as plt +from matplotlib.patches import Rectangle +from matplotlib.figure import Figure +from matplotlib.axes import Axes +from vtkmodules.vtkCommonDataModel import ( vtkCellTypes, VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, + VTK_HEXAHEDRON, VTK_WEDGE, VTK_POLYGON, VTK_POLYHEDRON ) +from geos.mesh.processing.meshQualityMetricHelpers import ( QUALITY_METRIC_OTHER_START_INDEX, getAllCellTypesExtended, + getQualityMeasureNameFromIndex, getQualityMetricFromIndex, + MeshQualityMetricEnum, CellQualityMetricEnum, + VtkCellQualityMetricEnum, CellQualityMetricAdditionalEnum, + QualityMetricOtherEnum, QualityRange ) +from geos.mesh.model.CellTypeCounts import CellTypeCounts + +__doc__ = """ +QualityMetricSummary stores the statistics of mesh quality metrics. + +To use QualityMetricSummary: + +.. code-block:: python + + from geos.mesh.model.QualityMetricSummary import QualityMetricSummary, StatTypes + + qualityMetricSummary: QualityMetricSummary = QualityMetricSummary() + # set data + qualityMetricSummary.setCellTypeCounts(counts) + qualityMetricSummary.setCellStatValueFromStatMetricAndCellType(cellMetricIndex, cellType, statType, value)) + qualityMetricSummary.setOtherStatValueFromMetric(metricIndex, statType, value) + + # get stats + count: int = qualityMetricSummary.getCellTypeCountsOfCellType(cellType) + value: float = qualityMetricSummary.getCellStatValueFromStatMetricAndCellType(cellMetricIndex, cellType, statType) + subSetStats: pd.DataFrame = stats.getStatsFromMetric(cellMetricIndex) + stats: pd.DataFrame = stats.getAllCellStats() + + # output figure + fig: Figure = stats.plotSummaryFigure() +""" + + +class StatTypes( Enum ): + MEAN = ( 0, "Mean", float, np.nanmean ) + STD_DEV = ( 1, "StdDev", float, np.nanstd ) + MIN = ( 2, "Min", float, np.nanmin ) + MAX = ( 3, "Max", float, np.nanmax ) + COUNT = ( 4, "Count", int, lambda v: np.count_nonzero( np.isfinite( v ) ) ) + + def getIndex( self: Self ) -> int: + """Get stat index. + + Returns: + int: Index + """ + return self.value[ 0 ] + + def getString( self: Self ) -> str: + """Get stat name. + + Returns: + str: Name + """ + return self.value[ 1 ] + + def getType( self: Self ) -> object: + """Get stat type. + + Returns: + object: Type + """ + return self.value[ 2 ] + + def compute( self: Self, array: Iterable[ float ] ) -> int | float: + """Compute statistics using function. + + Args: + array (Iterable[float]): Input array + + Returns: + int | float: Output stat + """ + return self.value[ 3 ]( array ) + + @staticmethod + def getNameFromIndex( index: int ) -> str: + """Get stat name from index. + + Args: + index (int): Index + + Returns: + str: Name + """ + return list( StatTypes )[ index ].getString() + + @staticmethod + def getIndexFromName( name: str ) -> int: + """Get stat index from name. + + Args: + name (str): Name + + Returns: + int: Index + """ + for stat in list( StatTypes ): + if stat.getString() == name: + return stat.getIndex() + return -1 + + @staticmethod + def getTypeFromIndex( index: int ) -> object: + """Get stat type from index. + + Args: + index (int): Index + + Returns: + object: Type + """ + return list( StatTypes )[ index ].getType() + + +_RANGE_COLORS: tuple[ str, str, str ] = ( + 'lightcoral', + 'sandybrown', + 'palegreen', +) + + +class QualityMetricSummary(): + + _LEVELS: tuple[ str, str ] = ( "MetricIndex", "CellType" ) + _CELL_TYPES_PLOT: tuple[ int, ...] = ( VTK_TRIANGLE, VTK_QUAD, VTK_POLYGON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, + VTK_HEXAHEDRON, VTK_POLYHEDRON ) + _CELL_TYPES_NAME: list[ str ] = [ + vtkCellTypes.GetClassNameFromTypeId( cellType ).removeprefix( "vtk" ) for cellType in _CELL_TYPES_PLOT + ] + + def __init__( self: Self ) -> None: + """CellTypeCounts stores the number of cells of each type.""" + #: stores for each cell type, and metric type, the stats + self._counts: CellTypeCounts = CellTypeCounts() + self._cellStats: pd.DataFrame + self._meshOtherStats: pd.DataFrame + self._initStats() + + def _initStats( self: Self ) -> None: + """Initialize self._cellStats dataframe.""" + rows: list[ int ] = [ statType.getIndex() for statType in list( StatTypes ) ] + nb_rows: int = len( rows ) + cellTypes: list[ int ] = getAllCellTypesExtended() + indexes = [ ( metric.getMetricIndex(), cellType ) + for metric in ( list( VtkCellQualityMetricEnum ) + list( CellQualityMetricAdditionalEnum ) ) + for cellType in cellTypes if metric.isApplicableToCellType( cellType ) ] + df_indexes: pd.MultiIndex = pd.MultiIndex.from_tuples( ( indexes ), names=self._LEVELS ) + nb_col: int = df_indexes.size + self._cellStats = pd.DataFrame( np.full( ( nb_rows, nb_col ), np.nan ), columns=df_indexes, index=rows ) + + columns = [ metric.getMetricIndex() for metric in list( QualityMetricOtherEnum ) ] + self._meshOtherStats = pd.DataFrame( np.full( ( nb_rows, len( columns ) ), np.nan ), + columns=columns, + index=rows ) + + def setCellTypeCounts( self: Self, counts: CellTypeCounts ) -> None: + """Set cell type counts. + + Args: + counts (CellTypeCounts): CellTypeCounts instance + """ + self._counts = counts + + def getCellTypeCountsObject( self: Self ) -> int: + """Get cell type counts. + + Returns: + int: Number of cell + """ + return self._counts + + def getCellTypeCountsOfCellType( self: Self, cellType: int ) -> int: + """Get cell type counts. + + Returns: + int: Number of cell + """ + return self._counts.getTypeCount( cellType ) + + def isCellStatsValidForMetricAndCellType( + self: Self, + metricIndex: int, + cellType: int, + ) -> bool: + """Returns True if input quality metric applies to input cell type. + + Args: + metricIndex (int): Metric index + cellType (int): Cell type index + + Returns: + bool: True if input quality metric applies + """ + return bool( np.any( np.isfinite( self.getStatsFromMetricAndCellType( metricIndex, cellType ) ) ) ) + + def getAllCellStats( self: Self ) -> pd.DataFrame: + """Get all cell stats including nan values. + + Returns: + pd.DataFrame: Stats + """ + return self._cellStats + + def getAllValidCellStats( self: Self ) -> pd.DataFrame: + """Get all valid cell stats. + + Returns: + pd.DataFrame: Stats + """ + return self._cellStats.dropna( axis=1 ) + + def getAllValidOtherMetricStats( self: Self ) -> pd.DataFrame: + """Get all valid other metric stats. + + Returns: + pd.DataFrame: Stats + """ + print( self._meshOtherStats.head() ) + return self._meshOtherStats.dropna( axis=1 ) + + def getCellStatValueFromStatMetricAndCellType( + self: Self, + metricIndex: int, + cellType: int, + statType: StatTypes, + ) -> float: + """Get cell stat value for the given metric and cell types. + + Args: + metricIndex (int): Metric index + cellType (int): Cell type index + statType (StatTypes): Stat number + + Returns: + float: Stats value + """ + if ( metricIndex, cellType ) not in self._cellStats.columns: + raise IndexError( f"Index ({metricIndex}, {cellType}) not in QualityMetricSummary stats" ) + return self._cellStats[ ( metricIndex, cellType ) ][ statType.getIndex() ] + + def getStatsFromMetricAndCellType( self: Self, metricIndex: int, cellType: int ) -> pd.Series: + """Get stats for the given metric and cell types. + + Args: + metricIndex (int): Metric index + cellType (int): Cell type index + + Returns: + pd.Series: Stats + """ + if ( metricIndex, cellType ) not in self._cellStats.columns: + raise IndexError( f"Index ({metricIndex}, {cellType}) not in QualityMetricSummary stats" ) + return self._cellStats[ ( metricIndex, cellType ) ] + + def getStatsFromMetric( + self: Self, + metricIndex: int, + ) -> pd.DataFrame: + """Get stats for the given metric index. + + Args: + metricIndex (int): Metric index + + Returns: + pd.DataFrame: Stats + """ + if metricIndex < QUALITY_METRIC_OTHER_START_INDEX: + return self._cellStats.xs( metricIndex, level=self._LEVELS[ 0 ], axis=1 ) + else: + return self._meshOtherStats[ metricIndex ] + + def setOtherStatValueFromMetric( self: Self, metricIndex: int, statType: StatTypes, value: int | float ) -> None: + """Set other stat value for the given metric. + + Args: + metricIndex (int): Metric index + statType (StatTypes): Stat number + value (int | float): Value + """ + if metricIndex not in self._meshOtherStats.columns: + raise IndexError( f"Index {metricIndex} not in QualityMetricSummary meshOtherStats" ) + self._meshOtherStats.loc[ statType.getIndex(), metricIndex ] = value + + def getCellStatsFromCellType( + self: Self, + cellType: int, + ) -> pd.DataFrame: + """Get cell stats for the given cell type. + + Args: + cellType (int): Cell type index + + Returns: + pd.DataFrame: Stats + """ + return self._cellStats.xs( cellType, level=self._LEVELS[ 1 ], axis=1 ) + + def setCellStatValueFromStatMetricAndCellType( self: Self, metricIndex: int, cellType: int, statType: StatTypes, + value: int | float ) -> None: + """Set cell stats for the given metric and cell types. + + Args: + metricIndex (int): Metric index + cellType (int): Cell type index + statType (StatTypes): Stat number + value (int | float): Value + """ + if ( metricIndex, cellType ) not in self._cellStats.columns: + raise IndexError( f"Index ({metricIndex}, {cellType}) not in QualityMetricSummary stats" ) + self._cellStats.loc[ statType.getIndex(), ( metricIndex, cellType ) ] = value + + def getComputedCellMetricIndexes( self: Self ) -> list[ Any ]: + """Get the list of index of computed cell quality metrics. + + Returns: + tuple[int]: List of metrics index + """ + validCellStats: pd.DataFrame = self.getAllValidCellStats() + columns: list[ int ] = validCellStats.columns.get_level_values( 0 ).to_list() + return np.unique( columns ).tolist() + + def getComputedOtherMetricIndexes( self: Self ) -> list[ Any ]: + """Get the list of index of computed other quality metrics. + + Returns: + tuple[int]: List of metrics index + """ + validOtherStats: pd.DataFrame = self.getAllValidOtherMetricStats() + columns: list[ int ] = [ validOtherStats.columns.to_list() ] + return np.unique( columns ).tolist() + + def getAllComputedMetricIndexes( self: Self ) -> list[ Any ]: + """Get the list of index of all computed metrics. + + Returns: + tuple[int]: List of metrics index + """ + return self.getComputedCellMetricIndexes() + self.getComputedOtherMetricIndexes() + + def plotSummaryFigure( self: Self ) -> Figure: + """Plot quality metric summary figure. + + Returns: + plt.figure: Output Figure + """ + computedCellMetrics: list[ int ] = self.getComputedCellMetricIndexes() + computedOtherMetrics: list[ int ] = self.getComputedOtherMetricIndexes() + # compute layout + nbAxes: int = len( computedCellMetrics ) + len( computedOtherMetrics ) + 1 + ncols: int = 3 + nrows: int = 1 + # 3 columns for these number of axes, else 4 columns + if nbAxes not in ( 1, 2, 3, 5, 6, 9 ): + ncols = 4 + nrows = nbAxes // ncols + if nbAxes % ncols > 0: + nrows += 1 + figSize = ( ncols * 3, nrows * 4 ) + fig, axes = plt.subplots( nrows, ncols, figsize=figSize, tight_layout=True ) + axesFlat = axes.flatten() + # index of current axes + currentAxIndex: int = 0 + + # plot cell type counts + self._plotCellTypeCounts( axesFlat[ 0 ] ) + currentAxIndex += 1 + + # plot other mesh quality stats + ax: Axes + if len( computedOtherMetrics ) > 0: + ax = axesFlat[ currentAxIndex ] + self._plotOtherMetricStats( ax ) + currentAxIndex += 1 + # plot cell quality metrics + for metricIndex in computedCellMetrics: + ax = axesFlat[ currentAxIndex ] + self._plotCellMetricStats( ax, metricIndex ) + currentAxIndex += 1 + + # remove unused axes + for ax in axesFlat[ currentAxIndex: ]: + ax.remove() + return fig + + def _plotCellTypeCounts( self: Self, ax: Axes ) -> None: + """Plot cell type counts. + + Args: + ax (Axes): Axes object + """ + xticks: npt.NDArray[ np.int64 ] = np.arange( len( self._CELL_TYPES_PLOT ), dtype=int ) + toplot: list[ int ] = [ self._counts.getTypeCount( cellType ) for cellType in self._CELL_TYPES_PLOT ] + p = ax.bar( self._CELL_TYPES_NAME, toplot ) + # bar_label only available for matplotlib version >= 3.3 + if Version( mpl.__version__ ) >= Version( "3.3" ): + plt.bar_label( p, label_type='center', rotation=90, padding=5 ) + ax.set_xticks( xticks ) + ax.set_xticklabels( self._CELL_TYPES_NAME, rotation=30, ha="right" ) + ax.set_xlabel( "Cell types" ) + ax.set_title( "Cell Type Counts" ) + + def _plotOtherMetricStats( self: Self, ax0: Axes ) -> None: + """Plot other metric stats. + + Args: + ax0 (Axes): Axes object + metricIndex (int): Metric index + """ + # order of cell types in each axes + computedMetrics: list[ int ] = self.getComputedOtherMetricIndexes() + # get data to plot + maxs: pd.Series = self._meshOtherStats.loc[ StatTypes.MAX.getIndex(), computedMetrics ] + mins: pd.Series = self._meshOtherStats.loc[ StatTypes.MIN.getIndex(), computedMetrics ] + means: pd.Series = self._meshOtherStats.loc[ StatTypes.MEAN.getIndex(), computedMetrics ] + stdDev: pd.Series = self._meshOtherStats.loc[ StatTypes.STD_DEV.getIndex(), computedMetrics ] + xticks: npt.NDArray[ np.int64 ] = np.arange( means.index.size, dtype=int ) + xtickslabels = [ getQualityMeasureNameFromIndex( metricIndex ) for metricIndex in computedMetrics ] + # define colors + cmap: mpl.colors.Colormap = cm.get_cmap( 'plasma' ) + colors = cmap( np.linspace( 0, 1, len( computedMetrics ) ) ) + + # min max rectangle width + recWidth: float = 0.5 + xtick: float = 0.0 + ax: Axes + for k, metricIndex in enumerate( computedMetrics ): + ax = ax0 if k == 0 else ax0.twinx() + color = colors[ k ] + # add rectangle from min to max + x: float = xtick - recWidth / 2.0 + y: float = mins[ metricIndex ] + recHeight: float = maxs[ metricIndex ] - mins[ metricIndex ] + ax.add_patch( Rectangle( ( x, y ), recWidth, recHeight, edgecolor=color, fill=False, lw=1 ) ) + + # plot mean and error bars for std dev + ax.errorbar( k, means[ metricIndex ], yerr=stdDev[ metricIndex ], fmt='-o', color=color ) + xtick += 1.0 + + # set y axis color + ax.yaxis.label.set_color( color ) + ax.tick_params( axis='y', colors=color ) + + # set x tick names + ax0.set_xticks( xticks ) + ax0.set_xticklabels( xtickslabels, rotation=30, ha="right" ) + ax0.set_xlabel( "Mesh Quality Metric" ) + ax0.set_title( "Other Mesh Quality Metric" ) + + def _plotCellMetricStats( self: Self, ax: Axes, metricIndex: int ) -> None: + """Plot cell metric stats. + + Args: + ax (Axes): Axes object + metricIndex (int): Metric index + """ + # get data to plot + maxs: pd.Series = self._cellStats.loc[ StatTypes.MAX.getIndex(), metricIndex ] + mins: pd.Series = self._cellStats.loc[ StatTypes.MIN.getIndex(), metricIndex ] + means: pd.Series = self._cellStats.loc[ StatTypes.MEAN.getIndex(), metricIndex ] + xticks: npt.NDArray[ np.int64 ] = np.arange( means.index.size, dtype=int ) + stdDev: pd.Series = self._cellStats.loc[ StatTypes.STD_DEV.getIndex(), metricIndex ] + + # order of cell types in each axes + xtickslabels: list[ str ] = [] + # min max rectangle width + recWidth: float = 0.5 + # range rectangle width + rangeRecWidth: float = 1.8 * recWidth + ylim0: float = mins.max() + ylim1: float = maxs.min() + xtick: float = 0.0 + for k, cellType in enumerate( self._CELL_TYPES_PLOT ): + if cellType in means.index: + xtickslabels += [ self._CELL_TYPES_NAME[ k ] ] + # add quality metric range + ( ylim0, ylim1 ) = self._plotRangePatch( ax, metricIndex, cellType, ylim0, ylim1, xtick, rangeRecWidth ) + # add rectangle from min to max + x: float = xtick - recWidth / 2.0 + y: float = mins[ cellType ] + recHeight: float = maxs[ cellType ] - mins[ cellType ] + ax.add_patch( Rectangle( ( x, y ), recWidth, recHeight, edgecolor='black', fill=False, lw=1 ) ) + # plot mean and error bars for std dev + ax.errorbar( xtick, means[ cellType ], yerr=stdDev[ cellType ], fmt='-o', color='k' ) + xtick += 1.0 + + # set y axis limits + ax.set_ylim( 0.1 * ylim0, 1.1 * ylim1 ) + # set x tick names + ax.set_xticks( xticks ) #, xtickslabels, rotation=30, ha="right") + ax.set_xticklabels( xtickslabels, rotation=30, ha="right" ) + ax.set_xlabel( "Cell types" ) + ax.set_title( f"{getQualityMeasureNameFromIndex(metricIndex)}" ) + + def _plotRangePatch( self: Self, ax: Axes, metricIndex: int, cellType: int, ylim0: float, ylim1: float, + xtick: float, rangeRecWidth: float ) -> tuple[ float, float ]: + """Plot quality metric ranges. + + Args: + ax (Axes): Axes object + metricIndex (int): Metric index + cellType (int): Cell type index + ylim0 (float): Min y + ylim1 (float): Max y + xtick (float): Abscissa + rangeRecWidth (float): Patch width + + Returns: + tuple[float, float]: Tuple containing min y and max y + """ + try: + metric: MeshQualityMetricEnum = getQualityMetricFromIndex( metricIndex ) + assert isinstance( metric, CellQualityMetricEnum ), "Mesh quality metric is of wrong type." + # add quality range patches if relevant + qualityRange: QualityRange | None = metric.getQualityRange( cellType ) + if qualityRange is not None: + ( ylim0, ylim1 ) = self._plotQualityRange( ax, qualityRange, xtick - rangeRecWidth / 2.0, + ( ylim0, ylim1 ), rangeRecWidth ) + else: + # add white patch for tick alignment + ax.add_patch( + Rectangle( + ( xtick - rangeRecWidth / 2.0, 0. ), + rangeRecWidth, + 1.0, + facecolor='w', + fill=True, + ) ) + except AssertionError as e: + print( "Cell quality metric range cannot be displayed due to: ", e ) + return ( ylim0, ylim1 ) + + def _plotQualityRange( self: Self, ax: Axes, qualityRange: QualityRange, x: float, ylim: tuple[ float, float ], + rangeRecWidth: float ) -> tuple[ float, float ]: + """Plot quality range patches. + + Args: + ax (Axes): Axes object + qualityRange (QualityRange): Quality ranges to plot + x (float): Origin abscissa of the patches + ylim (tuple[float, float]): Y limits for updates + rangeRecWidth (float): Patch width + + Returns: + tuple[float, float]: Y limits for updates + """ + ylim0: float = ylim[ 0 ] + ylim1: float = ylim[ 1 ] + for k, ( vmin, vmax ) in enumerate( + ( qualityRange.fullRange, qualityRange.normalRange, qualityRange.acceptableRange ) ): + if not np.isfinite( vmin ): + vmin = -1e12 + else: + ylim0 = min( ylim0, vmin ) + if not np.isfinite( vmax ): + vmax = 1e12 + else: + ylim1 = max( ylim1, vmax ) + y: float = vmin + recHeight = vmax - vmin + ax.add_patch( Rectangle( + ( x, y ), + rangeRecWidth, + recHeight, + facecolor=_RANGE_COLORS[ k ], + fill=True, + ) ) + return ( ylim0, ylim1 ) diff --git a/geos-mesh/src/geos/mesh/model/__init__.py b/geos-mesh/src/geos/mesh/model/__init__.py new file mode 100644 index 00000000..b7db2541 --- /dev/null +++ b/geos-mesh/src/geos/mesh/model/__init__.py @@ -0,0 +1 @@ +# Empty diff --git a/geos-mesh/src/geos/mesh/processing/SplitMesh.py b/geos-mesh/src/geos/mesh/processing/SplitMesh.py index 95a219bc..0275774e 100644 --- a/geos-mesh/src/geos/mesh/processing/SplitMesh.py +++ b/geos-mesh/src/geos/mesh/processing/SplitMesh.py @@ -30,7 +30,7 @@ from vtkmodules.util.numpy_support import ( numpy_to_vtk, vtk_to_numpy ) -from geos.mesh.stats.CellTypeCounter import CellTypeCounter +from geos.mesh.stats.CellTypeCounterEnhanced import CellTypeCounterEnhanced from geos.mesh.model.CellTypeCounts import CellTypeCounts __doc__ = """ @@ -187,10 +187,10 @@ def _get_cell_counts( self: Self ) -> CellTypeCounts: Returns: CellTypeCounts: cell type counts """ - filter: CellTypeCounter = CellTypeCounter() + filter: CellTypeCounterEnhanced = CellTypeCounterEnhanced() filter.SetInputDataObject( self.inData ) filter.Update() - return filter.GetCellTypeCounts() + return filter.GetCellTypeCountsObject() def _addMidPoint( self: Self, ptA: int, ptB: int ) -> int: """Add a point at the center of the edge defined by input point ids. diff --git a/geos-mesh/src/geos/mesh/processing/__init__.py b/geos-mesh/src/geos/mesh/processing/__init__.py new file mode 100644 index 00000000..b7db2541 --- /dev/null +++ b/geos-mesh/src/geos/mesh/processing/__init__.py @@ -0,0 +1 @@ +# Empty diff --git a/geos-mesh/src/geos/mesh/processing/meshQualityMetricHelpers.py b/geos-mesh/src/geos/mesh/processing/meshQualityMetricHelpers.py new file mode 100644 index 00000000..992f00ea --- /dev/null +++ b/geos-mesh/src/geos/mesh/processing/meshQualityMetricHelpers.py @@ -0,0 +1,690 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Antoine Mazuyer, Martin Lemay +from dataclasses import dataclass +import numpy as np +from typing import Optional +from typing_extensions import Self +from enum import Enum +from vtkmodules.vtkFiltersVerdict import vtkMeshQuality +from vtkmodules.vtkCommonDataModel import ( VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_HEXAHEDRON, VTK_WEDGE, + VTK_POLYGON, VTK_POLYHEDRON ) + +__doc__ = """ +Helpers for MeshQuality metrics. +""" + +#: start index of additional cell quality metrics +CELL_QUALITY_METRIC_ADDITIONAL_START_INDEX: int = 100 +#: start index of other mesh quality metrics +QUALITY_METRIC_OTHER_START_INDEX: int = 200 + + +@dataclass( frozen=True ) +class QualityRange(): + """Defines metric quality ranges.""" + acceptableRange: tuple[ float, float ] + normalRange: tuple[ float, float ] + fullRange: tuple[ float, float ] + + +class MeshQualityMetricEnum( Enum ): + + def __init__( + self: Self, + metricIndex: int, + name: str, + ) -> None: + """Define the enumeration to add attributes to mesh quality measures. + + Args: + metricIndex (int): Index of QualityMeasureTypes + name (str): Name of the metric + applicableToCellTypes (tuple[bool, ...]): Tuple defining for each cell type if the + metric is applicable. + """ + self.metricIndex: int = int( metricIndex ) + self.metricName: str = name + + def getMetricIndex( self: Self ) -> int: + """Get metric index. + + Returns: + int: Metric index + """ + return self.metricIndex + + def getMetricName( self: Self ) -> str: + """Get metric name. + + Returns: + str: Metric name + """ + return self.metricName + + +class CellQualityMetricEnum( MeshQualityMetricEnum ): + + def __init__( + self: Self, + metricIndex: int, + name: str, + applicableToCellTypes: tuple[ bool, ...], + qualityRanges: tuple[ QualityRange | None, ...], + ) -> None: + """Define the enumeration to add attributes to mesh quality measures. + + Args: + metricIndex (int): Index of QualityMeasureTypes + name (str): Name of the metric + applicableToCellTypes (tuple[bool, ...]): Tuple defining for each cell type if the + metric is applicable. + qualityRanges (tuple[QualityRange | None,...]): Quality range limits for each cell type + starting from best to worst quality. + """ + super().__init__( metricIndex, name ) + self.applicableToCellTypes: tuple[ bool, ...] = applicableToCellTypes + self.qualityRanges: tuple[ QualityRange | None, ...] = qualityRanges + + def getApplicableCellTypes( self: Self ) -> set[ int ]: + """Get the list of cell type indexes the metric applies to. + + Returns: + set[int]: Set of cell type indexes + """ + cellTypes = set() + for i, cellType in enumerate( getAllCellTypes() ): + if self.applicableToCellTypes[ i ]: + cellTypes.add( cellType ) + return cellTypes + + def isApplicableToCellType( self: Self, cellType: int ) -> bool: + """Return True if the metric is applicable to input cell type, False otherwise. + + Args: + cellType (int): Cell type index + + Returns: + bool: True if the metric is applicable + """ + assert cellType in getAllCellTypesExtended(), f"Cell type index {cellType} not in supported cell types." + cellTypes: list[ int ] = [ + cellType, + ] + if cellType == VTK_POLYGON: + cellTypes = getPolygonCellTypes() + if cellType == VTK_POLYHEDRON: + cellTypes = getPolyhedronCellTypes() + for cellType in cellTypes: + cellTypeIndex: int = getAllCellTypes().index( cellType ) + if not self.applicableToCellTypes[ cellTypeIndex ]: + return False + return True + + def getQualityRange( self: Self, cellType: int ) -> Optional[ QualityRange ]: + """Get quality range for input cell type. + + Args: + cellType (int): Cell type index + + Returns: + tuple[float, float, float, float, bool] | None: Quality range from best to worst. Last element + yields True if the range is symmetrical to negative values. + """ + if cellType not in getAllCellTypes(): + return None + cellTypeIndex: int = getAllCellTypes().index( cellType ) + return self.qualityRanges[ cellTypeIndex ] + + +#: from https://vtk.org/doc/nightly/html/vtkMeshQuality_8h_source.html +class VtkCellQualityMetricEnum( CellQualityMetricEnum ): + """Cell quality metric enumeration. + + The order of boolean is the same as getAllCellTypes method: + (VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON) + + .. CAUTION:: The order of the enum must follow the one of vtkMeshQuality.QualityMeasureTypes. + + """ + #: ratio of cell longest and shortest edge lengths + EDGE_RATIO = ( vtkMeshQuality.QualityMeasureTypes.EDGE_RATIO, "Edge Ratio", ( True, True, True, False, True, True ), + ( + QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), ( 1.0, np.inf ) ), + QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), ( 1.0, np.inf ) ), + QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), ( 1.0, np.inf ) ), + None, + QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), ( 1.0, np.inf ) ), + QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), ( 1.0, np.inf ) ), + ) ) + #: ratio of the maximum edge length to the inradius (for simplicial elements but adapted for quads). + #: may be adapted for polyhedron other than tet by splitting the polyhedron in tet and computing the + #: max of aspect ratio of all tet + ASPECT_RATIO = ( + vtkMeshQuality.QualityMeasureTypes.ASPECT_RATIO, + "Aspect Ratio", + ( True, True, True, False, False, False ), + ( QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), + ( 1.0, np.inf ) ), QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), ( 1.0, np.inf ) ), + QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), ( 1.0, np.inf ) ), None, None, None ), + ) + #: ratio between the radius of the inscribed circle/sphere to the radius of the circum-circle/sphere + #: normalized so that the ratio yields 1 for equilateral cell + RADIUS_RATIO = ( + vtkMeshQuality.QualityMeasureTypes.RADIUS_RATIO, + "Radius Ratio", + ( True, True, True, False, False, False ), + ( QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), + ( 1.0, np.inf ) ), QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), ( 1.0, np.inf ) ), + QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), ( 1.0, np.inf ) ), None, None, None ), + ) + #: sum of the edge lengths squared divided by the area for triangles. Adapted for Tetrahedron. + #: normalized so that equal to 1 when the element is equilateral triangle + ASPECT_FROBENIUS = ( + vtkMeshQuality.QualityMeasureTypes.ASPECT_FROBENIUS, + "Aspect Frobenius", + ( True, False, True, False, False, False ), + ( QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), + ( 1.0, np.inf ) ), None, QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), + ( 1.0, np.inf ) ), None, None, None ), + ) + #: median of Aspect Frobenius over all triangles of quads / tetrahedra of hexahedron + MEDIAN_ASPECT_FROBENIUS = ( vtkMeshQuality.QualityMeasureTypes.MED_ASPECT_FROBENIUS, "Med Aspect Frobenius", + ( False, True, False, False, False, True ), + ( None, QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), ( 1.0, np.inf ) ), None, None, None, + QualityRange( ( 1.0, 3.0 ), ( 1.0, 3.0 ), ( 9.0, np.inf ) ) ) ) + #: maximum of Aspect Frobenius over all triangles of Quads / tetrahedra of hexahedron + MAXIMUM_ASPECT_FROBENIUS = ( vtkMeshQuality.QualityMeasureTypes.MAX_ASPECT_FROBENIUS, "Maximum Aspect Frobenius", + ( False, True, False, False, True, True ), + ( None, QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), ( 1.0, np.inf ) ), None, None, + QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), ( 1.0, np.inf ) ), + QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), ( 1.0, np.inf ) ) ) ) + #: minimum angle between two neighboring edges for polygons / faces for tetrahedron. + MINIMUM_ANGLE = ( vtkMeshQuality.QualityMeasureTypes.MIN_ANGLE, "Minimum Angle (°)", ( True, True, True, False, + False, False ), + ( QualityRange( ( 30.0, 60.0 ), ( 0.0, 60.0 ), + ( 0.0, 360.0 ) ), QualityRange( ( 45.0, 90.0 ), ( 0.0, 90.0 ), ( 0.0, 360. ) ), + QualityRange( ( 40.0, 180. / np.pi * np.arccos( 1 / 3 ) ), + ( 0.0, 180. / np.pi * np.arccos( 1 / 3 ) ), ( 0.0, 360.0 ) ), None, None, None ) ) + #: the smallest ratio of the height of a vertex above its opposing triangle to + #: the longest edge of that opposing triangle across all vertices of the tetrahedron + COLLAPSE_RATIO = ( vtkMeshQuality.QualityMeasureTypes.COLLAPSE_RATIO, "Collapse Ratio", ( False, False, True, False, + False, False ), + ( None, None, QualityRange( ( 0.1, 1.0 ), ( 0.0, np.inf ), + ( 0.0, np.inf ) ), None, None, None ) ) + #: maximum angle between two neighboring edges for polygons / faces for tetrahedron. + MAXIMUM_ANGLE = ( vtkMeshQuality.QualityMeasureTypes.MAX_ANGLE, "Maximum Angle (°)", + ( True, True, False, False, False, False ), ( QualityRange( + ( 60., 90.0 ), ( 60.0, 180.0 ), + ( 0.0, 180.0 ) ), QualityRange( ( 90.0, 135.0 ), ( 90.0, 360.0 ), + ( 0.0, 360. ) ), None, None, None, None ) ) + #: condition number of the weighted Jacobian matrix. + CONDITION = ( + vtkMeshQuality.QualityMeasureTypes.CONDITION, + "Condition", + ( True, True, True, False, True, True ), + ( QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), + ( 1.0, np.inf ) ), QualityRange( ( 1.0, 4.0 ), ( 1.0, 12.0 ), ( 1.0, np.inf ) ), + QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), + ( 1.0, np.inf ) ), None, QualityRange( ( 1.0, 4.0 ), ( 1.0, 12.0 ), ( 1.0, np.inf ) ), + QualityRange( ( 1.0, 4.0 ), ( 1.0, 12.0 ), ( 1.0, np.inf ) ) ), + ) + #: Jacobian divided by the product of the lengths of the longest edges + #: normalized so that a unit equilateral triangle has value 1. + SCALED_JACOBIAN = ( + vtkMeshQuality.QualityMeasureTypes.SCALED_JACOBIAN, + "Scaled Jacobian", + ( True, True, True, True, True, True ), + ( + QualityRange( ( 0.5, 2.0 * np.sqrt( 3 ) / 3.0 ), ( -2.0 * np.sqrt( 3 ) / 3.0, 2.0 * np.sqrt( 3 ) / 3.0 ), + ( -np.inf, np.inf ) ), + QualityRange( ( 0.30, 1.0 ), ( -1.0, 1.0 ), ( -1.0, np.inf ) ), + QualityRange( ( 0.5, 0.5 * np.sqrt( 2 ) ), ( -0.5 * np.sqrt( 2 ), 0.5 * np.sqrt( 2 ) ), + ( -np.inf, np.inf ) ), + QualityRange( ( 0.50, 1.0 ), ( -1.0, 1.0 ), ( -1.0, np.inf ) ), + QualityRange( ( 0.50, 1.0 ), ( -1.0, 1.0 ), ( -1.0, np.inf ) ), + QualityRange( ( 0.50, 1.0 ), ( -1.0, 1.0 ), ( -1.0, np.inf ) ), + ), + ) + #: same as Scaled Jacobian + SHEAR = ( vtkMeshQuality.QualityMeasureTypes.SHEAR, "Shear", ( False, True, False, False, False, True ), + ( None, QualityRange( ( 0.3, 0.6 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), None, None, None, + QualityRange( ( 0.3, 0.6 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ) ) ) + #: the minimum of the ratio of cell area/volume to the average area/volume of an ensemble of cells and its inverse. + RELATIVE_SIZE_SQUARED = ( + vtkMeshQuality.QualityMeasureTypes.RELATIVE_SIZE_SQUARED, + "Relative Size Squared", + ( True, True, True, False, False, True ), + ( QualityRange( ( 0.25, 0.50 ), ( 0.0, 1.0 ), + ( 0.0, 1.0 ) ), QualityRange( ( 0.30, 0.6 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.30, 0.50 ), ( 0.0, 1.0 ), + ( 0.0, 1.0 ) ), None, None, QualityRange( ( 0.50, 1.0 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ) ), + ) + #: inverse of Condition (Polygons) / Jacobian ratio + SHAPE = ( + vtkMeshQuality.QualityMeasureTypes.SHAPE, + "Shape", + ( True, True, True, True, True, True ), + ( + QualityRange( ( 0.25, 0.50 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.30, 0.60 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.30, 0.60 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.30, 0.60 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.30, 0.60 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.30, 0.60 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + ), + ) + #: relative size squared times shape + SHAPE_AND_SIZE = ( + vtkMeshQuality.QualityMeasureTypes.SHAPE_AND_SIZE, + "Shape And Size", + ( True, True, True, False, False, True ), + ( + QualityRange( ( 0.25, 0.5 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.20, 0.4 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.20, 0.4 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.20, 0.4 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.20, 0.4 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.20, 0.4 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + ), + ) + #: measure of how well-behaved the mapping from parameter space to world coordinates is. ratio of the + #: minimum of Jacobian determinant to cell area/volume + DISTORTION = ( + vtkMeshQuality.QualityMeasureTypes.DISTORTION, + "Distortion", + ( True, True, True, False, True, True ), + ( + QualityRange( ( 0.5, 1.0 ), ( 0.0, 1.0 ), ( -np.inf, np.inf ) ), + QualityRange( ( 0.5, 1.0 ), ( 0.0, 1.0 ), ( -np.inf, np.inf ) ), + QualityRange( ( 0.5, 1.0 ), ( 0.0, 1.0 ), ( -np.inf, np.inf ) ), + None, + QualityRange( ( 0.5, 1.0 ), ( 0.0, 1.0 ), ( -np.inf, np.inf ) ), + QualityRange( ( 0.5, 1.0 ), ( 0.0, 1.0 ), ( -np.inf, np.inf ) ), + ), + ) + #: maximum of edge ratio over all triangles of the cell + MAXIMUM_EDGE_RATIO = ( vtkMeshQuality.QualityMeasureTypes.MAX_EDGE_RATIO, "Maximum Edge Ratio", + ( False, True, False, False, False, True ), + ( None, QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), ( 1.0, np.inf ) ), None, None, None, + QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), ( 1.0, np.inf ) ) ) ) + #: measures the angle (absolute value of the cosine) between the principal axes. + SKEW = ( vtkMeshQuality.QualityMeasureTypes.SKEW, "Skew", ( False, True, False, False, False, True ), + ( None, QualityRange( ( 0.5, 1.0 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), None, None, None, + QualityRange( ( 0.0, 0.5 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ) ) ) + #: maximum ratio of cross derivative magnitude to principal axis magnitude + TAPER = ( vtkMeshQuality.QualityMeasureTypes.TAPER, "Taper", ( False, True, False, False, False, True ), + ( None, QualityRange( ( 0.0, 0.7 ), ( 0.0, 2.0 ), ( 0.0, np.inf ) ), None, None, None, + QualityRange( ( 0.0, 0.5 ), ( 0.0, 1.5 ), ( 0.0, np.inf ) ) ) ) + #: polyhedron volume + VOLUME = ( + vtkMeshQuality.QualityMeasureTypes.VOLUME, + "Volume (m3)", + ( False, False, True, True, True, True ), + ( None, None, QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), ( -np.inf, np.inf ) ), + QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), + ( -np.inf, np.inf ) ), QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), ( -np.inf, np.inf ) ), + QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), ( -np.inf, np.inf ) ) ), + ) + #: ratio of minimum edge length to longest diagonal length + STRETCH = ( vtkMeshQuality.QualityMeasureTypes.STRETCH, "Stretch", ( False, True, False, False, False, True ), + ( None, QualityRange( ( 0.25, 0.5 ), ( 0.0, 1.0 ), ( 0.0, np.inf ) ), None, None, None, + QualityRange( ( 0.25, 0.5 ), ( 0.0, 1.0 ), ( 0.0, np.inf ) ) ) ) + #: ratio of the minimum diagonal length to the maximum diagonal length + DIAGONAL = ( vtkMeshQuality.QualityMeasureTypes.DIAGONAL, "Diagonal", ( False, False, False, False, False, True ), ( + None, + None, + None, + None, + None, + QualityRange( ( 0.65, 1.0 ), ( 0.0, 1.0 ), ( 0.0, np.inf ) ), + ) ) + # ratio of cell volume to the gradient (divergence?) of cell volumes at the cell. Acceptable range is application dependent. + DIMENSION = ( vtkMeshQuality.QualityMeasureTypes.DIMENSION, "Dimension (m)", + ( False, False, False, False, False, True ), ( + None, + None, + None, + None, + None, + QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), ( 0.0, np.inf ) ), + ) ) + #: measures the maximum deviation of the metric tensor at the corners of the quadrilateral. Maximum of oddy for hexahedron. + ODDY = ( vtkMeshQuality.QualityMeasureTypes.ODDY, "Oddy", ( False, True, False, False, False, True ), ( + None, + QualityRange( ( 0.0, 0.5 ), ( 0.0, 1.5 ), ( 0.0, np.inf ) ), + None, + None, + None, + QualityRange( ( 0.0, 0.5 ), ( 0.0, 1.5 ), ( 0.0, np.inf ) ), + ) ) + #: relative size squared times shear + SHEAR_AND_SIZE = ( vtkMeshQuality.QualityMeasureTypes.SHEAR_AND_SIZE, "Shear And Size", ( False, True, False, False, + False, True ), + ( None, QualityRange( ( 0.2, 0.4 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), None, None, None, + QualityRange( ( 0.2, 0.4 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ) ) ) + #: minimum determinant of the Jacobian matrix evaluated at each corner and the center of the element + JACOBIAN = ( vtkMeshQuality.QualityMeasureTypes.JACOBIAN, "Jacobian", ( False, True, True, True, True, True ), ( + None, + QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), ( -np.inf, np.inf ) ), + QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), ( -np.inf, np.inf ) ), + QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), ( -np.inf, np.inf ) ), + QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), ( -np.inf, np.inf ) ), + QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), ( -np.inf, np.inf ) ), + ) ) + #: the cosine of the minimum dihedral angle formed by planes intersecting in diagonals (to the fourth power) + WARPAGE = ( vtkMeshQuality.QualityMeasureTypes.WARPAGE, "Warpage", ( False, True, False, False, False, False ), + ( None, QualityRange( ( 0.0, 0.7 ), ( 0.0, 2.0 ), ( 0.0, np.inf ) ), None, None, None, None ) ) + #: ratio of root-mean-square edge length to volume. + #: normalizing the metric to a value of 1 for equilateral tetrahedra + ASPECT_GAMMA = ( vtkMeshQuality.QualityMeasureTypes.ASPECT_GAMMA, "Aspect Gamma", ( False, False, True, False, + False, False ), + ( None, None, QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), ( 0.0, np.inf ) ), None, None, None ) ) + #: polygon area + AREA = ( + vtkMeshQuality.QualityMeasureTypes.AREA, + "Area (m2)", + ( True, True, False, False, False, False ), + ( QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), + ( -np.inf, np.inf ) ), QualityRange( ( 0.0, np.inf ), ( 0.0, np.inf ), + ( -np.inf, np.inf ) ), None, None, None, None ), + ) + #: maximum of ratio of angular deviation from ideal element + EQUIANGLE_SKEW = ( vtkMeshQuality.QualityMeasureTypes.EQUIANGLE_SKEW, "Equiangle Skew", + ( True, True, True, True, True, True ), ( + QualityRange( ( 0.0, 0.5 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.0, 0.5 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.0, 0.5 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.0, 0.5 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.0, 0.5 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + QualityRange( ( 0.0, 0.5 ), ( 0.0, 1.0 ), ( 0.0, 1.0 ) ), + ) ) + #: maximum of ratio of volume deviation from ideal element + EQUIVOLUME_SKEW = ( vtkMeshQuality.QualityMeasureTypes.EQUIVOLUME_SKEW, "Equivolume Skew", ( False, False, True, + False, False, False ), + ( None, None, QualityRange( ( 0.0, 0.3 ), ( 0.0, 0.9 ), ( 0.0, 1.0 ) ), None, None, None ) ) + #: maximum stretch over tetrahedra + MAXIMUM_STRETCH = ( vtkMeshQuality.QualityMeasureTypes.MAX_STRETCH, "Maximum Stretch", + ( False, False, False, False, True, False ), ( None, None, None, None, + QualityRange( ( 0.25, 0.5 ), ( 0.0, 1.0 ), + ( 0.0, np.inf ) ), None ) ) + #: mean of Aspect Frobenius over all triangles of Quads / tetrahedra of hexahedron + MEAN_ASPECT_FROBENIUS = ( vtkMeshQuality.QualityMeasureTypes.MEAN_ASPECT_FROBENIUS, "Mean Aspect Frobenius", + ( False, False, False, False, True, False ), ( None, None, None, None, + QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), + ( 1.0, np.inf ) ), None ) ) + #: ratio of tetrahedron volume over the volume of an equilateral tetrahedron with the same root mean squared edge length + MEAN_RATIO = ( vtkMeshQuality.QualityMeasureTypes.MEAN_RATIO, "Mean Ratio", + ( False, False, True, False, False, False ), ( None, None, + QualityRange( ( 0.0, 0.3 ), ( 0.0, 0.9 ), + ( 0.0, 1.0 ) ), None, None, None ) ) + #: ratio between the largest and smallest Jacobian determinant value + NODAL_JACOBIAN_RATIO = ( vtkMeshQuality.QualityMeasureTypes.NODAL_JACOBIAN_RATIO, "Nodal Jacobian Ratio", + ( False, False, False, False, False, True ), ( None, None, None, None, None, + QualityRange( ( 0.0, np.inf ), + ( 0.0, np.inf ), + ( 0.0, np.inf ) ) ) ) + #: ratio of the minimum sub-triangle inner radius to the outer triangle radius + NORMALIZED_INRADIUS = ( vtkMeshQuality.QualityMeasureTypes.NORMALIZED_INRADIUS, "Normalized Inradius", + ( True, False, True, False, False, False ), ( QualityRange( + ( 0.15, 0.5 ), ( -1.0, 1.0 ), + ( -1.0, 1.0 ) ), None, QualityRange( ( 0.15, 0.5 ), ( -1.0, 1.0 ), + ( -1.0, 1.0 ) ), None, None, None ) ) + #: measure used to quantify how far a cell deviates from orthogonality with respect to its face + #: maximum of sinus of the angle between the vector from polyhedron center and face center and face normal + #: yields 0 if vectors are parallel, 1 if they are orthogonal + SQUISH_INDEX = ( vtkMeshQuality.QualityMeasureTypes.SQUISH_INDEX, "Squish Index", + ( False, False, True, True, True, True ), ( None, None, + QualityRange( ( 0.0, 0.3 ), ( 0.0, 0.9 ), + ( 0.0, 1.0 ) ), None, None, None ) ) + #: no metric + NONE = ( vtkMeshQuality.QualityMeasureTypes.NONE, "None", ( False, False, False, False, False, False ), + ( None, None, None, None, None, None ) ) + + +class CellQualityMetricAdditionalEnum( CellQualityMetricEnum ): + """Additional cell quality metric enumeration. + + The order of boolean is the same as getAllCellTypes method: + (VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON) + + Metric index starts at 100 to prevent from conflicts with basic metrics and is incremented + in the order of the enumeration. + + """ + #: maximum of aspect ratio over all tets . + MAXIMUM_ASPECT_RATIO = ( + 100, + "Maximum Aspect Ratio", + ( False, False, False, True, True, True ), + ( QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), + ( 1.0, np.inf ) ), QualityRange( ( 1.0, 1.3 ), ( 1.0, 3.0 ), ( 1.0, np.inf ) ), + QualityRange( ( 1.0, 3.0 ), ( 1.0, 9.0 ), ( 1.0, np.inf ) ), None, None, None ), + ) + # other metrics can be defined the same way such as RADIUS_RATIO, EDGE_RATIO, ASPECT_GAMMA, EQUIVOLUME_SKEW, + # NORMALIZED_INRADIUS, (MAXIMUM_)ASPECT_FROBENIUS for pyramids, (MAXIMUM_)STRETCH for pyramids + + +class QualityMetricOtherEnum( MeshQualityMetricEnum ): + """Additional metrics that apply to the mesh, not to specific cell type. + + Metric index starts at 200 to prevent from conflicts with other metrics and is incremented + in the order of the enumeration. + """ + #: number of incident edges for each vertex + INCIDENT_VERTEX_COUNT = ( 200, "Incident Vertex Count" ) + + +def getAllCellTypesExtended() -> list[ int ]: + """Get all cell type ids. + + Returns: + tuple[int,...]: Tuple containing cell type ids. + """ + return getPolygonCellTypes() + getPolyhedronCellTypes() + [ VTK_POLYGON, VTK_POLYHEDRON ] + + +def getAllCellTypes() -> list[ int ]: + """Get all cell type ids. + + Returns: + tuple[int,...]: Tuple containing cell type ids. + """ + return getPolygonCellTypes() + getPolyhedronCellTypes() + + +def getChildrenCellTypes( parent: int ) -> list[ int ]: + """Get children cell type ids from parent id. + + Returns: + tuple[int,...]: Tuple containing cell type ids. + """ + if parent == VTK_POLYGON: + return getPolygonCellTypes() + elif parent == VTK_POLYHEDRON: + return getPolyhedronCellTypes() + else: + raise ValueError( f"Cell type {parent} is not supported." ) + + +def getPolygonCellTypes() -> list[ int ]: + """Get polygonal cell type ids. + + Returns: + tuple[int,...]: Tuple containing cell type ids. + """ + return [ VTK_TRIANGLE, VTK_QUAD ] + + +def getPolyhedronCellTypes() -> list[ int ]: + """Get polyhedra cell type ids. + + Returns: + tuple[int,...]: Tuple containing cell type ids. + """ + return [ VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON ] + + +def getQualityMetricFromIndex( metricIndex: int ) -> Optional[ MeshQualityMetricEnum ]: + """Get quality metric from its index. + + Args: + metricIndex (int): Metric index + + Raises: + IndexError: Metric index is out of range + + Returns: + MeshQualityMetricEnum | None: Quality metric + """ + if metricIndex < CELL_QUALITY_METRIC_ADDITIONAL_START_INDEX: + return list( VtkCellQualityMetricEnum )[ metricIndex ] + elif metricIndex < QUALITY_METRIC_OTHER_START_INDEX: + return list( CellQualityMetricAdditionalEnum )[ metricIndex - CELL_QUALITY_METRIC_ADDITIONAL_START_INDEX ] + elif metricIndex < QUALITY_METRIC_OTHER_START_INDEX + len( list( QualityMetricOtherEnum ) ): + return list( QualityMetricOtherEnum )[ metricIndex - QUALITY_METRIC_OTHER_START_INDEX ] + return None + + +def getQualityMeasureNameFromIndex( metricIndex: int ) -> Optional[ str ]: + """Get quality metric name from index. + + Args: + metricIndex (int): Index of quality measure + + Returns: + str | None: Name of quality measure. Returns None if metricIndex is undefined. + """ + metric = getQualityMetricFromIndex( metricIndex ) + if metric is None: + return None + return metric.getMetricName() + + +def getQualityMeasureIndexFromName( name: str ) -> int: + """Get quality metric index from name. + + Args: + name (str): Name of quality measure + + Returns: + int: Index of quality measure + """ + for metric in list( VtkCellQualityMetricEnum ) + list( CellQualityMetricAdditionalEnum ) + list( + QualityMetricOtherEnum ): + if metric.getMetricName() == name: + return metric.getMetricIndex() + return int( vtkMeshQuality.QualityMeasureTypes.NONE ) + + +def getCellQualityMeasureFromCellType( cellType: int ) -> set[ int ]: + """Get the indexes of mesh quality metrics defined for triangles. + + Returns: + set[int]: Set of possible indexes. + """ + if cellType not in getAllCellTypesExtended(): + raise ValueError( f"Cell type {cellType} not in supported cell types {getAllCellTypesExtended()}." ) + return { + metric.metricIndex + for metric in list( VtkCellQualityMetricEnum ) + list( CellQualityMetricAdditionalEnum ) + if metric.isApplicableToCellType( cellType ) + } + + +def getTriangleQualityMeasure() -> set[ int ]: + """Get the indexes of mesh quality metrics defined for triangles. + + Returns: + set[int]: Set of possible indexes. + """ + return getCellQualityMeasureFromCellType( VTK_TRIANGLE ) + + +def getQuadQualityMeasure() -> set[ int ]: + """Get the indexes of mesh quality metrics defined for quads. + + Returns: + set[int]: Set of possible indexes. + """ + return getCellQualityMeasureFromCellType( VTK_QUAD ) + + +def getCommonPolygonQualityMeasure() -> set[ int ]: + """Get the indexes of mesh quality metrics defined for both triangles and quads. + + Returns: + set[int]: Set of possible indexes. + """ + triangleMetrics: set[ int ] = getTriangleQualityMeasure() + quadMetrics: set[ int ] = getQuadQualityMeasure() + return triangleMetrics.intersection( quadMetrics ) + + +def getTetQualityMeasure() -> set[ int ]: + """Get the indexes of mesh quality metrics defined for quads. + + Returns: + set[int]: Set of possible indexes. + """ + return getCellQualityMeasureFromCellType( VTK_TETRA ) + + +def getPyramidQualityMeasure() -> set[ int ]: + """Get the indexes of mesh quality metrics defined for quads. + + Returns: + set[int]: Set of possible indexes. + """ + return getCellQualityMeasureFromCellType( VTK_PYRAMID ) + + +def getWedgeQualityMeasure() -> set[ int ]: + """Get the indexes of mesh quality metrics defined for quads. + + Returns: + set[int]: Set of possible indexes. + """ + return getCellQualityMeasureFromCellType( VTK_WEDGE ) + + +def getHexQualityMeasure() -> set[ int ]: + """Get the indexes of mesh quality metrics defined for quads. + + Returns: + set[int]: Set of possible indexes. + """ + return getCellQualityMeasureFromCellType( VTK_HEXAHEDRON ) + + +def getCommonPolyhedraQualityMeasure() -> set[ int ]: + """Get the indexes of mesh quality metrics defined for both triangles and quads. + + Returns: + set[int]: Set of possible indexes. + """ + tetMetrics: set[ int ] = getTetQualityMeasure() + pyrMetrics: set[ int ] = getPyramidQualityMeasure() + wedgeMetrics: set[ int ] = getWedgeQualityMeasure() + hexMetrics: set[ int ] = getHexQualityMeasure() + return tetMetrics.intersection( pyrMetrics ).intersection( wedgeMetrics ).intersection( hexMetrics ) + + +def getQualityMetricsOther() -> set[ int ]: + """Get the set of indexes of other mesh quality metric. + + Returns: + set[int]: Other mesh quality metric indexes + """ + return { metric.getMetricIndex() for metric in list( QualityMetricOtherEnum ) } + + +#: dictionary of cell quality metrics set from cell type +cellQualityMetricsFromCellType: dict[ int, set[ int ] ] = { + VTK_TRIANGLE: getTriangleQualityMeasure(), + VTK_QUAD: getQuadQualityMeasure(), + VTK_TETRA: getTetQualityMeasure(), + VTK_PYRAMID: getPyramidQualityMeasure(), + VTK_WEDGE: getWedgeQualityMeasure(), + VTK_HEXAHEDRON: getHexQualityMeasure(), +} \ No newline at end of file diff --git a/geos-mesh/src/geos/mesh/stats/CellTypeCounterEnhanced.py b/geos-mesh/src/geos/mesh/stats/CellTypeCounterEnhanced.py new file mode 100644 index 00000000..533d759d --- /dev/null +++ b/geos-mesh/src/geos/mesh/stats/CellTypeCounterEnhanced.py @@ -0,0 +1,112 @@ +# 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 ) + +from geos.mesh.model.CellTypeCounts import CellTypeCounts +from geos.mesh.processing.meshQualityMetricHelpers import getAllCellTypes + +__doc__ = """ +CellTypeCounterEnhanced 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.CellTypeCounterEnhanced import CellTypeCounterEnhanced + + # filter inputs + input :vtkUnstructuredGrid + + # instantiate the filter + filter :CellTypeCounterEnhanced = CellTypeCounterEnhanced() + # set input data object + filter.SetInputDataObject(input) + # do calculations + filter.Update() + # get counts + counts :CellTypeCounts = filter.GetCellTypeCountsObject() +""" + + +class CellTypeCounterEnhanced( VTKPythonAlgorithmBase ): + + def __init__( self ) -> None: + """CellTypeCounterEnhanced filter computes mesh stats.""" + super().__init__( nInputPorts=1, nOutputPorts=1, inputType="vtkUnstructuredGrid", outputType="vtkTable" ) + self._counts: CellTypeCounts = 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.reset() + 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 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 GetCellTypeCountsObject( self: Self ) -> CellTypeCounts: + """Get CellTypeCounts object. + + Returns: + CellTypeCounts: CellTypeCounts object. + """ + return self._counts \ No newline at end of file diff --git a/geos-mesh/src/geos/mesh/stats/MeshQualityEnhanced.py b/geos-mesh/src/geos/mesh/stats/MeshQualityEnhanced.py new file mode 100644 index 00000000..420787fe --- /dev/null +++ b/geos-mesh/src/geos/mesh/stats/MeshQualityEnhanced.py @@ -0,0 +1,803 @@ +# 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 import Optional, cast +from typing_extensions import Self +from vtkmodules.vtkFiltersCore import vtkExtractEdges +from vtkmodules.vtkFiltersVerdict import vtkMeshQuality +from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase +from vtkmodules.vtkCommonCore import ( + vtkInformation, + vtkInformationVector, + vtkIdList, + vtkPoints, + vtkDataArray, + vtkIntArray, + vtkDoubleArray, + vtkIdTypeArray, + vtkMath, +) +from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkPolyData, vtkCellData, vtkPointData, vtkFieldData, + vtkCell, vtkCell3D, vtkTetra, vtkCellTypes, vtkPolygon, VTK_TRIANGLE, + VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_HEXAHEDRON, VTK_WEDGE, VTK_POLYGON, + VTK_POLYHEDRON ) +from vtkmodules.util.numpy_support import vtk_to_numpy, numpy_to_vtk + +from geos.mesh.stats.CellTypeCounterEnhanced import CellTypeCounterEnhanced +from geos.mesh.model.CellTypeCounts import CellTypeCounts +from geos.mesh.model.QualityMetricSummary import QualityMetricSummary, StatTypes +from geos.mesh.utils.arrayHelpers import getAttributesFromDataSet +from geos.mesh.processing.meshQualityMetricHelpers import ( + getQualityMeasureNameFromIndex, + getQualityMetricFromIndex, + cellQualityMetricsFromCellType, + VtkCellQualityMetricEnum, + CellQualityMetricAdditionalEnum, + QualityMetricOtherEnum, + MeshQualityMetricEnum, + getAllCellTypesExtended, + getAllCellTypes, + getPolygonCellTypes, + getPolyhedronCellTypes, + getCellQualityMeasureFromCellType, + getChildrenCellTypes, +) + +import geos.utils.geometryFunctions as geom + +__doc__ = """ +MeshQualityEnhanced module is a vtk filter that computes mesh quality stats. + +Mesh quality stats include those from vtkMeshQuality as well as . + +Filter input is a vtkUnstructuredGrid. + +To use the filter: + +.. code-block:: python + + from geos.mesh.stats.MeshQualityEnhanced import MeshQualityEnhanced + + # filter inputs + input :vtkUnstructuredGrid + + # instanciate the filter + filter :MeshQualityEnhanced = MeshQualityEnhanced() + # set input data object + filter.SetInputDataObject(input) + # set metrics to use + filter.SetTriangleMetrics(triangleQualityMetrics) + filter.SetQuadMetrics(quadQualityMetrics) + filter.SetTetraMetrics(tetraQualityMetrics) + filter.SetPyramidMetrics(pyramidQualityMetrics) + filter.SetWedgeMetrics(wedgeQualityMetrics) + filter.SetHexaMetrics(hexaQualityMetrics) + filter.SetOtherMeshQualityMetrics(otherQualityMetrics) + # do calculations + filter.Update() + # get output mesh quality report + outputMesh: vtkUnstructuredGrid = filter.GetOutputDataObject(0) + outputStats: QualityMetricSummary = filter.GetQualityMetricSummary() +""" + +#: name of output quality array from vtkMeshQuality filter +QUALITY_ARRAY_NAME: str = "Quality" + + +def getQualityMetricArrayName( metric: int ) -> str: + """Get the name of the array from quality metric index. + + Args: + metric (int): Metric index + + Returns: + str: name of output array + """ + return QUALITY_ARRAY_NAME + "_" + "".join( getQualityMeasureNameFromIndex( metric ).split( " " ) ) + + +class MeshQualityEnhanced( VTKPythonAlgorithmBase ): + + def __init__( self: Self ) -> None: + """Enhanced vtkMeshQuality filter.""" + super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkUnstructuredGrid" ) + self._outputMesh: vtkUnstructuredGrid + self._cellCounts: CellTypeCounts + self._qualityMetricSummary: QualityMetricSummary = QualityMetricSummary() + + self._MetricsAll: dict[ int, Optional[ set[ int ] ] ] = { + VTK_TRIANGLE: None, + VTK_QUAD: None, + VTK_TETRA: None, + VTK_PYRAMID: None, + VTK_WEDGE: None, + VTK_HEXAHEDRON: None, + } + self._otherMetrics: Optional[ set[ QualityMetricOtherEnum ] ] = None + # for each cell, save cell type for later use + self._cellTypeMask: dict[ int, npt.NDArray[ np.bool_ ] ] = {} + + # static members that can be loaded once to save computational times + self._allCellTypesExtended: tuple[ int, ...] = getAllCellTypesExtended() + self._allCellTypes: tuple[ int, ...] = getAllCellTypes() + + 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, + inInfoVec: list[ vtkInformationVector ], + outInfoVec: vtkInformationVector, + ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestDataObject. + + Args: + request (vtkInformation): Request + inInfoVec (list[vtkInformationVector]): Input objects + outInfoVec (vtkInformationVector): Output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + inData = self.GetInputData( inInfoVec, 0, 0 ) + outData = self.GetOutputData( outInfoVec, 0 ) + assert inData is not None + if outData is None or ( not outData.IsA( inData.GetClassName() ) ): + outData = inData.NewInstance() + outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) + return super().RequestDataObject( request, inInfoVec, outInfoVec ) # type: ignore[no-any-return] + + def GetQualityMetricSummary( self: Self ) -> QualityMetricSummary: + """Get QualityMetricSummary object. + + Returns: + QualityMetricSummary: QualityMetricSummary object + """ + return self._qualityMetricSummary + + def SetTriangleMetrics( self: Self, metrics: Optional[ set[ int ] ] ) -> None: + """Set triangle quality metrics to compute. + + Args: + metrics (Iterable[int]): Metrics to compute + """ + self._MetricsAll[ VTK_TRIANGLE ] = metrics + + def SetQuadMetrics( self: Self, metrics: Optional[ set[ int ] ] = None ) -> None: + """Set triangle quality metrics to compute. + + Args: + metrics (Iterable[int]): Metrics to compute + """ + self._MetricsAll[ VTK_QUAD ] = metrics + + def SetTetraMetrics( self: Self, metrics: Optional[ set[ int ] ] = None ) -> None: + """Set triangle quality metrics to compute. + + Args: + metrics (Iterable[int]): Metrics to compute + """ + self._MetricsAll[ VTK_TETRA ] = metrics + + def SetPyramidMetrics( self: Self, metrics: Optional[ set[ int ] ] = None ) -> None: + """Set triangle quality metrics to compute. + + Args: + metrics (Iterable[int]): Metrics to compute + """ + self._MetricsAll[ VTK_PYRAMID ] = metrics + + def SetWedgeMetrics( self: Self, metrics: Optional[ set[ int ] ] = None ) -> None: + """Set triangle quality metrics to compute. + + Args: + metrics (Iterable[int]): Metrics to compute + """ + self._MetricsAll[ VTK_WEDGE ] = metrics + + def SetHexaMetrics( self: Self, metrics: Optional[ set[ int ] ] = None ) -> None: + """Set triangle quality metrics to compute. + + Args: + metrics (Iterable[int]): Metrics to compute + """ + self._MetricsAll[ VTK_HEXAHEDRON ] = metrics + + def SetCellQualityMetrics( + self: Self, + triangleMetrics: Optional[ set[ int ] ] = None, + quadMetrics: Optional[ set[ int ] ] = None, + tetraMetrics: Optional[ set[ int ] ] = None, + pyramidMetrics: Optional[ set[ int ] ] = None, + wedgeMetrics: Optional[ set[ int ] ] = None, + hexaMetrics: Optional[ set[ int ] ] = None, + ) -> None: + """Set all quality metrics to compute. + + Args: + triangleMetrics (Iterable[int]): Triangle metrics to compute. + + Defaults to [vtkMeshQuality.QualityMeasureTypes.NONE,]. + quadMetrics (Iterable[int]): Quad metrics to compute. + + Defaults to [vtkMeshQuality.QualityMeasureTypes.NONE,]. + tetraMetrics (Iterable[int]): Tetrahedron metrics to compute. + + Defaults to [vtkMeshQuality.QualityMeasureTypes.NONE,]. + pyramidMetrics (Iterable[int]): Pyramid metrics to compute. + + Defaults to [vtkMeshQuality.QualityMeasureTypes.NONE,]. + wedgeMetrics (Iterable[int]): Wedge metrics to compute. + + Defaults to [vtkMeshQuality.QualityMeasureTypes.NONE,]. + hexaMetrics (Iterable[int]): Hexahedron metrics to compute. + + Defaults to [vtkMeshQuality.QualityMeasureTypes.NONE,]. + """ + self.SetTriangleMetrics( triangleMetrics ) + self.SetQuadMetrics( quadMetrics ) + self.SetTetraMetrics( tetraMetrics ) + self.SetPyramidMetrics( pyramidMetrics ) + self.SetWedgeMetrics( wedgeMetrics ) + self.SetHexaMetrics( hexaMetrics ) + + def SetOtherMeshQualityMetrics( self: Self, metrics: set[ QualityMetricOtherEnum ] ) -> None: + """Set additional metrics unrelated to cell types. + + Args: + metrics (set[QualityMetricOtherEnum]): Set of QualityMetricOtherEnum + """ + if len( metrics ) > 0: + self._otherMetrics = metrics + else: + self._otherMetrics = None + + def getComputedMetricsFromCellType( self: Self, cellType: int ) -> Optional[ set[ int ] ]: + """Get the set of metrics computed for input cell type. + + Args: + cellType (int): Cell type index + + Returns: + Optional[set[int]]: Set of computed quality metrics + """ + # child cell type + if cellType in self._allCellTypes: + return self._MetricsAll[ cellType ] + # for parent cell types, gather children metrics + metrics: set[ int ] | None = getCellQualityMeasureFromCellType( cellType ) + if metrics is None: + return None + commonComputedMetricsExists: bool = False + for cellTypeChild in getChildrenCellTypes( cellType ): + computedMetrics: set[ int ] | None = self._MetricsAll[ cellTypeChild ] + if computedMetrics is None: + continue + commonComputedMetricsExists = True + metrics = metrics.intersection( computedMetrics ) + return metrics if commonComputedMetricsExists else None + + 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 ) + self._outputMesh = self.GetOutputData( outInfoVec, 0 ) + assert inData is not None, "Input mesh is undefined." + assert self._outputMesh is not None, "Ouput pipeline is undefined." + self._outputMesh.ShallowCopy( inData ) + + # compute cell type counts + self._computeCellTypeCounts() + + # compute metrics and associated attributes + self._evaluateMeshQualityAll() + + # compute stats summary + self._updateStatsSummary() + + # create field data + self._createFieldDataStatsSummary() + + self._outputMesh.Modified() + return 1 + + def _computeCellTypeCounts( self: Self ) -> None: + """Compute cell type counts.""" + filter: CellTypeCounterEnhanced = CellTypeCounterEnhanced() + filter.SetInputDataObject( self._outputMesh ) + filter.Update() + counts: CellTypeCounts = filter.GetCellTypeCountsObject() + assert counts is not None, "CellTypeCounts is undefined" + self._qualityMetricSummary.setCellTypeCounts( counts ) + + def _evaluateMeshQualityAll( self: Self ) -> None: + """Compute all mesh quality metrics.""" + for cellType, metrics in self._MetricsAll.items(): + if metrics is None: + continue + for metricIndex in metrics: + self._evaluateCellQuality( metricIndex, cellType ) + + if self._otherMetrics is not None: + if QualityMetricOtherEnum.INCIDENT_VERTEX_COUNT.getMetricIndex() in self._otherMetrics: + self._countVertexIncidentEdges() + else: + # TODO: add other metrics + print( "" ) + + def _evaluateCellQuality( self: Self, metricIndex: int, cellType: int ) -> None: + """Compute mesh input quality metric for input cell type. + + Args: + metricIndex (int): Quality metric index + cellType (int): Cell type index + """ + arrayName: str = getQualityMetricArrayName( metricIndex ) + if arrayName in getAttributesFromDataSet( self._outputMesh, False ): + # metric is already computed (by default computed for all cell types if applicable) + return + + # get the list of cell types the metric applies to and check if these cell types are present + metric: MeshQualityMetricEnum | None = getQualityMetricFromIndex( metricIndex ) + if metric is None: + return + cellTypes: Optional[ set[ int ] ] = metric.getApplicableCellTypes() + if cellTypes is None: + return + nbCells: int = 0 + for cellType in cellTypes: + nbCells += self._qualityMetricSummary.getCellTypeCountsOfCellType( cellType ) + if nbCells == 0: + return + + # compute quality metric + output: vtkUnstructuredGrid | None = None + if ( metricIndex == VtkCellQualityMetricEnum.SQUISH_INDEX.metricIndex ): + # redefined Squish index calculation to be computed for any type of polyhedron + self._computeSquishIndex() + elif ( metricIndex in ( CellQualityMetricAdditionalEnum.MAXIMUM_ASPECT_RATIO.metricIndex, ) ): + # extended metric for any type of cells (other than tetra) from tetra metrics + self._computeAdditionalMetrics( metricIndex ) + else: + output = self._applyMeshQualityFilter( metricIndex, cellType ) + assert output is not None, "Output mesh from mesh quality calculation is undefined." + # transfer output cell array to input mesh + # TODO: to test if Shallow copy of vtkMeshQualityFilter result and rename "Quality" array is more efficient than what is done here + self._transferCellAttribute( output, QUALITY_ARRAY_NAME, arrayName, metricIndex ) + + def _applyMeshQualityFilter( self: Self, metric: int, cellType: int ) -> vtkUnstructuredGrid: + """Apply vtkMeshQuality filter. + + Args: + metric (int): Quality metric index + cellType (int): Cell type + + Returns: + vtkUnstructuredGrid: Filtered mesh + """ + meshQualityFilter = vtkMeshQuality() + meshQualityFilter.SetInputDataObject( self._outputMesh ) + if cellType == VTK_TRIANGLE: + meshQualityFilter.SetTriangleQualityMeasure( metric ) + elif cellType == VTK_QUAD: + meshQualityFilter.SetQuadQualityMeasure( metric ) + elif cellType == VTK_TETRA: + meshQualityFilter.SetTetQualityMeasure( metric ) + elif cellType == VTK_PYRAMID: + meshQualityFilter.SetPyramidQualityMeasure( metric ) + elif cellType == VTK_WEDGE: + meshQualityFilter.SetWedgeQualityMeasure( metric ) + elif cellType == VTK_HEXAHEDRON: + meshQualityFilter.SetHexQualityMeasure( metric ) + else: + print( "Cell type is not supported." ) + meshQualityFilter.SaveCellQualityOn() + meshQualityFilter.Update() + return meshQualityFilter.GetOutputDataObject( 0 ) + + def _computeAdditionalMetrics( self: Self, metricIndex: int ) -> None: + """Compute additional metrics from metrics defined for tetrahedron. + + Output is an cell array in output mesh. + + Args: + metricIndex (int): Metric index + """ + metric = getQualityMetricFromIndex( metricIndex ) + assert metric is not None, f"Additional cell quality metric index {metricIndex} is undefined." + # output array + name: str = getQualityMetricArrayName( metric.getMetricIndex() ) + newArray: vtkDoubleArray = vtkDoubleArray() + newArray.SetName( name ) + newArray.SetNumberOfValues( self._outputMesh.GetNumberOfCells() ) + newArray.SetNumberOfComponents( 1 ) + for i in range( self._outputMesh.GetNumberOfCells() ): + cell: vtkCell = self._outputMesh.GetCell( i ) + val: float = self._computeAdditionalMetricsCell( metricIndex, cell ) + newArray.InsertNextValue( val ) + # add array + cellArrays: vtkCellData = self._outputMesh.GetCellData() + assert cellArrays is not None, "Cell data from output mesh is undefined." + cellArrays.AddArray( newArray ) + cellArrays.Modified() + self._outputMesh.Modified() + return + + def _transferCellAttribute( + self: Self, + serverMesh: vtkUnstructuredGrid, + serverAttributeName: str, + clientAttributeName: str, + qualityMetric: int, + ) -> bool: + """Transfer quality attribute to the client mesh. + + The attribute is renamed with quality metric name. Because a default quality + metric is computed if an element does not support the desired metric, this + default metric is replaced by nan values. + + Args: + serverMesh (vtkUnstructuredGrid): Server mesh where Quality metric is + serverAttributeName (str): Name of the attribute in the server mesh + clientAttributeName (str): Name of the attribute in the client mesh + qualityMetric (int): Index of quality metric. + + Returns: + bool: True if the attribute was successfully transfered, False otherwise + """ + cellArrays: vtkCellData = serverMesh.GetCellData() + assert cellArrays is not None, "Cell data from vtkMeshQuality output mesh is undefined." + array: vtkDataArray = cellArrays.GetArray( serverAttributeName ) + assert array is not None, f"{serverAttributeName} attribute is undefined." + # rename array + array.SetName( clientAttributeName ) + # replace irrelevant values + self._replaceIrrelevantValues( array, serverMesh, qualityMetric ) + + # add array to input mesh + inputCellArrays: vtkCellData = self._outputMesh.GetCellData() + assert inputCellArrays is not None, "Cell data from input mesh is undefined." + inputCellArrays.AddArray( array ) + inputCellArrays.Modified() + return True + + def _replaceIrrelevantValues( self: Self, array: vtkDataArray, mesh: vtkUnstructuredGrid, + qualityMetric: int ) -> None: + """Replace irrelevant values. + + Values are irrelevant when a quality metric is computed + whereas input metric does not applies to the cell type. + + Args: + array (vtkDataArray): Array to update + mesh (vtkUnstructuredGrid): Mesh + qualityMetric (int): Quality metric index + """ + for cellId in range( mesh.GetNumberOfCells() ): + cell: vtkCell = mesh.GetCell( cellId ) + cellType: int = cell.GetCellType() + cellTypeQualityMetrics: set[ int ] = cellQualityMetricsFromCellType[ cellType ] + if ( qualityMetric > -1 ) and ( qualityMetric not in cellTypeQualityMetrics ): + array.SetTuple1( cellId, np.nan ) + + def _updateStatsSummary( self: Self ) -> None: + """Compute quality metric statistics.""" + # init cell type masks + self._initCellTypeMasks() + # stats for each cell types individually + count: int = 0 + metrics: set[ int ] | None + for cellType, metrics in self._MetricsAll.items(): + count = self._qualityMetricSummary.getCellTypeCountsOfCellType( cellType ) + if ( count == 0 ) or ( metrics is None ): + continue + for metricIndex in metrics: + self._updateStatsSummaryByCellType( metricIndex, cellType ) + + # stats for polygons and polyhedra + for cellType in ( VTK_POLYGON, VTK_POLYHEDRON ): + count = self._qualityMetricSummary.getCellTypeCountsOfCellType( cellType ) + # get common computed metrics + metrics = self.getComputedMetricsFromCellType( cellType ) + if ( count == 0 ) or ( metrics is None ): + continue + for metricIndex in metrics: + self._updateStatsSummaryByCellType( metricIndex, cellType ) + + def _updateStatsSummaryByCellType( self: Self, metricIndex: int, cellType: int ) -> None: + """Compute input quality metric statistics for input cell types. + + Args: + metricIndex (int): Quality metric index + cellType (int): Cell type index + """ + cellArrays: vtkCellData = self._outputMesh.GetCellData() + assert cellArrays is not None, "Cell data from input mesh is undefined." + arrayName: str = getQualityMetricArrayName( metricIndex ) + array: vtkDataArray | None = cellArrays.GetArray( arrayName ) + if array is None: + return + npArray: npt.NDArray[ np.float64 ] = vtk_to_numpy( array ) + cellTypeMask: npt.NDArray[ np.bool_ ] = self._cellTypeMask[ cellType ] + + self._qualityMetricSummary.setCellStatValueFromStatMetricAndCellType( + metricIndex, cellType, StatTypes.MEAN, StatTypes.MEAN.compute( npArray[ cellTypeMask ] ) ) + self._qualityMetricSummary.setCellStatValueFromStatMetricAndCellType( + metricIndex, cellType, StatTypes.STD_DEV, StatTypes.STD_DEV.compute( npArray[ cellTypeMask ] ) ) + self._qualityMetricSummary.setCellStatValueFromStatMetricAndCellType( + metricIndex, cellType, StatTypes.MIN, StatTypes.MIN.compute( npArray[ cellTypeMask ] ) ) + self._qualityMetricSummary.setCellStatValueFromStatMetricAndCellType( + metricIndex, cellType, StatTypes.MAX, StatTypes.MAX.compute( npArray[ cellTypeMask ] ) ) + self._qualityMetricSummary.setCellStatValueFromStatMetricAndCellType( + metricIndex, cellType, StatTypes.COUNT, StatTypes.COUNT.compute( npArray[ cellTypeMask ] ) ) + + def _initCellTypeMasks( self: Self ) -> None: + """Init _cellTypeMask variable.""" + # compute cell type masks + self._cellTypeMask = { + cellType: np.zeros( self._outputMesh.GetNumberOfCells(), dtype=bool ) + for cellType in self._allCellTypesExtended + } + polyhedronCellTypes: tuple[ int, ...] = getPolyhedronCellTypes() + polygonCellTypes: tuple[ int, ...] = getPolygonCellTypes() + for cellId in range( self._outputMesh.GetNumberOfCells() ): + for cellType in self._allCellTypesExtended: + cellTypes: tuple[ int, ...] = ( cellType, ) + if cellType == VTK_POLYGON: + cellTypes = polygonCellTypes + elif cellType == VTK_POLYHEDRON: + cellTypes = polyhedronCellTypes + self._cellTypeMask[ cellType ][ cellId ] = self._outputMesh.GetCellType( cellId ) in cellTypes + + def _createFieldDataStatsSummary( self: Self ) -> None: + """Create field data arrays with quality statistics.""" + fieldData: vtkFieldData = self._outputMesh.GetFieldData() + assert fieldData is not None, "Field data is undefined." + for cellType in self._allCellTypesExtended: + count: int = self._qualityMetricSummary.getCellTypeCountsOfCellType( cellType ) + metrics: Optional[ set[ int ] ] = self.getComputedMetricsFromCellType( cellType ) + # create count array + name = "_".join( ( vtkCellTypes.GetClassNameFromTypeId( cellType ), StatTypes.COUNT.getString() ) ) + countArray: vtkIntArray = vtkIntArray() + countArray.SetName( name ) + countArray.SetNumberOfValues( 1 ) + countArray.SetNumberOfComponents( 1 ) + countArray.SetValue( 0, count ) + fieldData.AddArray( countArray ) + if ( count == 0 ) or ( metrics is None ): + continue + + # create metric arrays + for metricIndex in metrics: + # one array per statistic number except Count (last one) + for statType in list( StatTypes )[ :-1 ]: + value: int = self._qualityMetricSummary.getCellStatValueFromStatMetricAndCellType( + metricIndex, cellType, statType ) + name = self._createArrayName( cellType, metricIndex, statType ) + metricArray: vtkDoubleArray = vtkDoubleArray() + metricArray.SetName( name ) + metricArray.SetNumberOfValues( 1 ) + metricArray.SetNumberOfComponents( 1 ) + metricArray.SetValue( 0, value ) + fieldData.AddArray( metricArray ) + fieldData.Modified() + + def _createArrayName( self: Self, cellType: int, metricIndex: int, statType: StatTypes ) -> str: + """Concatenate cell type, metric name, and statistic name in array name. + + Args: + cellType (int): Cell type index + metricIndex (int): Quality metric index + statType (StatTypes): Statistic type + + Returns: + str: array name + """ + return "_".join( ( vtkCellTypes.GetClassNameFromTypeId( cellType ), + getQualityMeasureNameFromIndex( metricIndex ).replace( " ", "" ), statType.getString() ) ) + + def _computeAdditionalMetricsCell( self: Self, metricIndex: int, cell: vtkCell ) -> float: + """Compute additional metrics from metrics defined for tetrahedron for a cell. + + Args: + metricIndex (int): Metric index + cell (vtkCell): Cell + + Returns: + float: outout value + """ + meshQualityFilter: vtkMeshQuality = vtkMeshQuality() + # triangulate cell faces + listSimplexPts = vtkPoints() + idList = vtkIdList() + cell.Triangulate( 1, idList, listSimplexPts ) + + simplexAspectRatio: list[ float ] = [] + index: int = 0 + while index != listSimplexPts.GetNumberOfPoints(): + # create tetra + tetra: vtkTetra = vtkTetra() + tetraPts: vtkPoints = tetra.GetPoints() + for i in range( 4 ): + tetraPts.SetPoint( i, listSimplexPts.GetPoint( index ) ) + tetraPts.Modified() + index += 1 + # compute aspect ratio of tetra + if metricIndex == CellQualityMetricAdditionalEnum.MAXIMUM_ASPECT_RATIO.getMetricIndex(): + simplexAspectRatio.append( meshQualityFilter.TetAspectRatio( tetra ) ) + else: + # metric is not supported + simplexAspectRatio.append( np.nan ) + if any( np.isfinite( simplexAspectRatio ) ): + return np.nanmax( simplexAspectRatio ) + return np.nan + + def _countVertexIncidentEdges( self: Self ) -> None: + """Compute edge length and vertex incident edge number.""" + metric: QualityMetricOtherEnum = QualityMetricOtherEnum.INCIDENT_VERTEX_COUNT + # edge are extracted as "cell" of dimension 1 + extractEdges = vtkExtractEdges() + extractEdges.SetInputData( self._outputMesh ) + extractEdges.Update() + polyData: vtkPolyData = extractEdges.GetOutput() + incidentCounts: npt.NDArray[ np.int64 ] = np.zeros( self._outputMesh.GetNumberOfPoints(), dtype=int ) + for edg in range( polyData.GetNumberOfCells() ): + if polyData.GetCell( edg ).GetCellDimension() != 1: + # not an edge + continue + + edgesPointIds: vtkIdList = polyData.GetCell( edg ).GetPointIds() + for i in range( edgesPointIds.GetNumberOfIds() ): + incidentCounts[ edgesPointIds.GetId( i ) ] += 1 + + # create point attribute + pointData: vtkPointData = self._outputMesh.GetPointData() + assert pointData is not None, "Point data is undefined." + countArray: vtkIntArray = numpy_to_vtk( incidentCounts, deep=1 ) + metricName: str = metric.getMetricName().replace( " ", "" ) + name: str = QUALITY_ARRAY_NAME + "_" + metricName + countArray.SetName( name ) + pointData.AddArray( countArray ) + pointData.Modified() + + fieldData: vtkFieldData = self._outputMesh.GetFieldData() + assert fieldData is not None, "Field data is undefined." + for statType in list( StatTypes ): + name = metricName + "_" + statType.getString() + val: float | int = statType.compute( incidentCounts ) + # add values to quality summary stats + self._qualityMetricSummary.setOtherStatValueFromMetric( metric.getMetricIndex(), statType, val ) + metricArray: vtkDoubleArray = vtkDoubleArray() + metricArray.SetName( name ) + metricArray.SetNumberOfValues( 1 ) + metricArray.SetNumberOfComponents( 1 ) + metricArray.SetValue( 0, val ) # type: ignore[arg-type] + fieldData.AddArray( metricArray ) + fieldData.Modified() + self._outputMesh.Modified() + + def _computeSquishIndex( self: Self ) -> None: + """Compute Squish index for all element type. + + Squish index is the maximum value of the cosine of the deviation angle between + cell center to face center vector and face normal vector. + + Output is a new cell array. + """ + # output array + name: str = getQualityMetricArrayName( VtkCellQualityMetricEnum.SQUISH_INDEX.getMetricIndex() ) + newArray: vtkDoubleArray = vtkDoubleArray() + newArray.SetName( name ) + newArray.SetNumberOfValues( self._outputMesh.GetNumberOfCells() ) + newArray.SetNumberOfComponents( 1 ) + # copy input data to prevent modifications from GetCellNeighbors method + copyData: vtkUnstructuredGrid = vtkUnstructuredGrid() + copyData.ShallowCopy( self._outputMesh ) + points: vtkPoints = copyData.GetPoints() + for c in range( copyData.GetNumberOfCells() ): + cell: vtkCell = copyData.GetCell( c ) + # applies only to polyhedra + if cell.GetCellDimension() != 3: + continue + # get cell center + cellCenter: npt.NDArray[ np.float64 ] = self._getCellCenter( cell ) + # compute deviation cosine for each face + squishIndex: npt.NDArray[ np.float64 ] = np.full( cell.GetNumberOfFaces(), np.nan ) + for f in range( cell.GetNumberOfFaces() ): + face: vtkCell = cell.GetFace( f ) + # get face center + ptsIds: vtkIdTypeArray = vtkIdTypeArray() + ptsIds.Allocate( face.GetNumberOfPoints() ) + ptsIdsList: vtkIdList = face.GetPointIds() + for i in range( ptsIdsList.GetNumberOfIds() ): + ptsIds.InsertNextValue( ptsIdsList.GetId( i ) ) + faceCenter: npt.NDArray[ np.float64 ] = self._getCellCenter( face, ptsIds, points ) + faceNormal: npt.NDArray[ np.float64 ] = self._getNormalVector( points, face ) + + vec: npt.NDArray[ np.float64 ] = cellCenter - faceCenter + angle: float = vtkMath.AngleBetweenVectors( vec, faceNormal ) # type: ignore[arg-type] + squishIndex[ f ] = np.sin( angle ) + newArray.InsertValue( c, np.nanmax( squishIndex ) ) + + # add array + cellArrays: vtkCellData = self._outputMesh.GetCellData() + assert cellArrays is not None, "Cell data from output mesh is undefined." + cellArrays.AddArray( newArray ) + cellArrays.Modified() + self._outputMesh.Modified() + + def _getCellCenter( self: Self, + cell: vtkCell, + ptsIds: Optional[ vtkIdTypeArray ] = None, + points: Optional[ vtkPoints ] = None ) -> npt.NDArray[ np.float64 ]: + """Compute cell center. + + Args: + cell (vtkCell): Input cell + ptsIds (vtkIdTypeArray | None): Cell point ids. Defaults to None. + points (vtkPoints | None): Mesh point coordinates. Defaults to None. + + Returns: + npt.NDArray[np.float64]: output cell center + """ + cellCenter: npt.NDArray[ np.float64 ] = np.zeros( 3 ) + if cell.GetCellDimension() == 2: + # polygonal cell + assert ptsIds is not None, "Point ids are required for computing polygonal cell center." + assert points is not None, "Points are required for computing polygonal cell center." + cell.GetPointIds() + vtkPolygon.ComputeCentroid( ptsIds, points, cellCenter ) # type: ignore[call-overload] + elif cell.GetCellDimension() == 3: + # volume cell + cell3D: vtkCell3D = cast( vtkCell3D, cell ) + cell3D.GetCentroid( cellCenter ) # type: ignore[arg-type] + else: + raise TypeError( "Cell must be polygonal or volumic." ) + return cellCenter + + def _getNormalVector( self: Self, points: vtkPoints, face: vtkCell ) -> npt.NDArray[ np.float64 ]: + """Get the normal to the input face. + + .. NOTE:: This method consider the faces as planes. + + Args: + points (vtkPoints): Mesh points + face (vtkCell): Face + + Returns: + npt.NDArray[np.float64]: Coordinates of the normal vector + """ + assert face.GetCellDimension() == 2, "Cell must be a planar polygon." + facePtsIds: vtkIdList = face.GetPointIds() + # need only 3 points among all to get the normal of the face since we suppose face is a plane + ptsCoords: npt.NDArray[ np.float64 ] = np.zeros( ( 3, 3 ), dtype=float ) + for i in range( 3 ): + points.GetPoint( facePtsIds.GetId( i ), ptsCoords[ i ] ) + return geom.computeNormalFromPoints( ptsCoords[ 0 ], ptsCoords[ 1 ], ptsCoords[ 2 ] ) \ No newline at end of file diff --git a/geos-mesh/src/geos/mesh/stats/__init__.py b/geos-mesh/src/geos/mesh/stats/__init__.py new file mode 100644 index 00000000..b7db2541 --- /dev/null +++ b/geos-mesh/src/geos/mesh/stats/__init__.py @@ -0,0 +1 @@ +# Empty diff --git a/geos-mesh/src/geos/mesh/utils/__init__.py b/geos-mesh/src/geos/mesh/utils/__init__.py index e69de29b..b7db2541 100644 --- a/geos-mesh/src/geos/mesh/utils/__init__.py +++ b/geos-mesh/src/geos/mesh/utils/__init__.py @@ -0,0 +1 @@ +# Empty diff --git a/geos-mesh/src/geos/mesh/utils/arrayHelpers.py b/geos-mesh/src/geos/mesh/utils/arrayHelpers.py index 3139d67f..32824aa8 100644 --- a/geos-mesh/src/geos/mesh/utils/arrayHelpers.py +++ b/geos-mesh/src/geos/mesh/utils/arrayHelpers.py @@ -571,66 +571,6 @@ def AsDF( surface: vtkPolyData, attributeNames: tuple[ str, ...] ) -> pd.DataFra return data -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 computeCellCenterCoordinates( mesh: vtkDataSet ) -> vtkDataArray: """Get the coordinates of Cell center. diff --git a/geos-mesh/tests/conftest.py b/geos-mesh/tests/conftest.py index 56a1de08..2e3231fc 100644 --- a/geos-mesh/tests/conftest.py +++ b/geos-mesh/tests/conftest.py @@ -44,7 +44,10 @@ def _get_dataset( datasetType: str ): vtkFilename = "data/domain_res5_id.vtu" elif datasetType == "polydata": reader: vtkXMLUnstructuredGridReader = vtkXMLUnstructuredGridReader() - vtkFilename = "data/surface.vtu" + vtkFilename = "data/triangulatedSurface.vtu" + elif datasetType == "tetraVolume": + reader: vtkXMLUnstructuredGridReader = vtkXMLUnstructuredGridReader() + vtkFilename = "data/tetraVolume.vtu" datapath: str = os.path.join( os.path.dirname( os.path.realpath( __file__ ) ), vtkFilename ) reader.SetFileName( datapath ) reader.Update() diff --git a/geos-mesh/tests/data/tetraVolume.vtu b/geos-mesh/tests/data/tetraVolume.vtu new file mode 100644 index 00000000..4b619c92 Binary files /dev/null and b/geos-mesh/tests/data/tetraVolume.vtu differ diff --git a/geos-mesh/tests/data/surface.vtu b/geos-mesh/tests/data/triangulatedSurface.vtu similarity index 100% rename from geos-mesh/tests/data/surface.vtu rename to geos-mesh/tests/data/triangulatedSurface.vtu diff --git a/geos-mesh/tests/test_CellTypeCounterEnhanced.py b/geos-mesh/tests/test_CellTypeCounterEnhanced.py new file mode 100644 index 00000000..d948d430 --- /dev/null +++ b/geos-mesh/tests/test_CellTypeCounterEnhanced.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.CellTypeCounterEnhanced import CellTypeCounterEnhanced +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_CellTypeCounterEnhanced_single( test_case: TestCase ) -> None: + """Test of CellTypeCounterEnhanced filter. + + Args: + test_case (TestCase): Test case + """ + filter: CellTypeCounterEnhanced = CellTypeCounterEnhanced() + filter.SetInputDataObject( test_case.mesh ) + filter.Update() + countsObs: CellTypeCounts = filter.GetCellTypeCountsObject() + 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_CellTypeCounterEnhanced_multi( test_case: TestCase ) -> None: + """Test of CellTypeCounterEnhanced filter. + + Args: + test_case (TestCase): Test case + """ + filter: CellTypeCounterEnhanced = CellTypeCounterEnhanced() + filter.SetInputDataObject( test_case.mesh ) + filter.Update() + countsObs: CellTypeCounts = filter.GetCellTypeCountsObject() + 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}." \ No newline at end of file diff --git a/geos-mesh/tests/test_CellTypeCounts.py b/geos-mesh/tests/test_CellTypeCounts.py index 32d3c7ea..b5278ce0 100644 --- a/geos-mesh/tests/test_CellTypeCounts.py +++ b/geos-mesh/tests/test_CellTypeCounts.py @@ -102,11 +102,7 @@ def test_CellTypeCounts_init() -> None: @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 - """ + """Test of CellTypeCounts init method.""" counts: CellTypeCounts = CellTypeCounts() for _ in range( test_case.nbVertex ): counts.addType( VTK_VERTEX ) diff --git a/geos-mesh/tests/test_MeshQualityEnhanced.py b/geos-mesh/tests/test_MeshQualityEnhanced.py new file mode 100644 index 00000000..9b69d08e --- /dev/null +++ b/geos-mesh/tests/test_MeshQualityEnhanced.py @@ -0,0 +1,162 @@ +# 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 matplotlib.figure import Figure +from dataclasses import dataclass +import numpy as np +import pandas as pd +import pytest +from typing import ( + Iterator, + Optional, +) + +from geos.mesh.processing.meshQualityMetricHelpers import ( + getAllCellTypesExtended, ) +from geos.mesh.stats.MeshQualityEnhanced import MeshQualityEnhanced +from geos.mesh.model.QualityMetricSummary import QualityMetricSummary + +from vtkmodules.vtkFiltersVerdict import vtkMeshQuality +from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkCellData, vtkFieldData, vtkCellTypes, VTK_TRIANGLE, + VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON ) +from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridReader + +# input data +meshName_all: tuple[ str, ...] = ( + "polydata", + "tetraVolume", +) +cellTypes_all: set[ int ] = ( VTK_TRIANGLE, VTK_TETRA ) +qualityMetrics_all: tuple[ set[ int ], ...] = ( + ( int( vtkMeshQuality.QualityMeasureTypes.ASPECT_RATIO ), int( vtkMeshQuality.QualityMeasureTypes.SCALED_JACOBIAN ), + int( vtkMeshQuality.QualityMeasureTypes.MAX_ANGLE ) ), + ( int( vtkMeshQuality.QualityMeasureTypes.SCALED_JACOBIAN ), + int( vtkMeshQuality.QualityMeasureTypes.EQUIANGLE_SKEW ), + int( vtkMeshQuality.QualityMeasureTypes.SQUISH_INDEX ) ), +) +cellTypeCounts_all: tuple[ tuple[ int, ...], ...] = ( ( + 26324, + 0, + 0, + 0, + 0, + 0, + 26324, + 0, +), ( + 0, + 0, + 368606, + 0, + 0, + 0, + 0, + 368606, +) ) +metricsSummary_all = ( + ( ( 1.07, 0.11, 1.00, 1.94, 26324.0 ), ( 1.07, 0.11, 1.00, 1.94, 26324.0 ), ( 64.59, 6.73, 60.00, 110.67, + 26324.0 ) ), + ( ( 1.71, 0.32, 1.02, 3.3, 368606.0 ), ( 1.71, 0.32, 1.02, 3.3, 368606.0 ), ( 0.65, 0.15, 0.05, 0.94, 368606.0 ) ), +) + + +@dataclass( frozen=True ) +class TestCase: + """Test case.""" + __test__ = False + #: mesh + meshName: str + cellType: vtkCellTypes + qualityMetrics: set[ int ] + cellTypeCounts: tuple[ int ] + metricsSummary: tuple[ float ] + + +def __generate_test_data() -> Iterator[ TestCase ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: Iterator on test cases + """ + for meshName, cellType, qualityMetrics, cellTypeCounts, metricsSummary in zip( meshName_all, + cellTypes_all, + qualityMetrics_all, + cellTypeCounts_all, + metricsSummary_all, + strict=True ): + yield TestCase( meshName, cellType, qualityMetrics, cellTypeCounts, metricsSummary ) + + +ids: list[ str ] = [ os.path.splitext( name )[ 0 ] for name in meshName_all ] + + +@pytest.mark.parametrize( "test_case", __generate_test_data(), ids=ids ) +def test_MeshQualityEnhanced( test_case: TestCase, dataSetTest: vtkUnstructuredGrid ) -> None: + """Test of CellTypeCounterEnhanced filter. + + Args: + test_case (TestCase): Test case + dataSetTest: vtkUnstructuredGrid + """ + mesh: vtkUnstructuredGrid = dataSetTest( test_case.meshName ) + filter: MeshQualityEnhanced = MeshQualityEnhanced() + filter.SetInputDataObject( mesh ) + if test_case.cellType == VTK_TRIANGLE: + filter.SetTriangleMetrics( test_case.qualityMetrics ) + elif test_case.cellType == VTK_QUAD: + filter.SetQuadMetrics( test_case.qualityMetrics ) + elif test_case.cellType == VTK_TETRA: + filter.SetTetraMetrics( test_case.qualityMetrics ) + elif test_case.cellType == VTK_PYRAMID: + filter.SetPyramidMetrics( test_case.qualityMetrics ) + elif test_case.cellType == VTK_WEDGE: + filter.SetWedgeMetrics( test_case.qualityMetrics ) + elif test_case.cellType == VTK_HEXAHEDRON: + filter.SetHexaMetrics( test_case.qualityMetrics ) + filter.Update() + + # test method getComputedMetricsFromCellType + for i, cellType in enumerate( getAllCellTypesExtended() ): + metrics: Optional[ set[ int ] ] = filter.getComputedMetricsFromCellType( cellType ) + if test_case.cellTypeCounts[ i ] > 0: + assert metrics is not None, f"Metrics from {vtkCellTypes.GetClassNameFromTypeId(cellType)} cells is undefined." + + # test attributes + outputMesh: vtkUnstructuredGrid = filter.GetOutputDataObject( 0 ) + cellData: vtkCellData = outputMesh.GetCellData() + assert cellData is not None, "Cell data is undefined." + + nbMetrics: int = len( test_case.qualityMetrics ) + nbCellArrayExp: int = mesh.GetCellData().GetNumberOfArrays() + nbMetrics + assert cellData.GetNumberOfArrays() == nbCellArrayExp, f"Number of cell arrays is expected to be {nbCellArrayExp}." + + # test field data + fieldData: vtkFieldData = outputMesh.GetFieldData() + assert fieldData is not None, "Field data is undefined." + tmp = np.array( test_case.cellTypeCounts ) > 0 + nbPolygon: int = np.sum( tmp[ :2 ].astype( int ) ) + nbPolygon = 0 if nbPolygon == 0 else nbPolygon + 1 + nbPolyhedra: int = np.sum( tmp[ 2:6 ].astype( int ) ) + nbPolyhedra = 0 if nbPolyhedra == 0 else nbPolyhedra + 1 + nbFieldArrayExp: int = mesh.GetFieldData().GetNumberOfArrays() + tmp.size + 4 * nbMetrics * ( nbPolygon + + nbPolyhedra ) + assert fieldData.GetNumberOfArrays( + ) == nbFieldArrayExp, f"Number of field data arrays is expected to be {nbFieldArrayExp}." + + stats: QualityMetricSummary = filter.GetQualityMetricSummary() + for i, cellType in enumerate( getAllCellTypesExtended() ): + # test Counts + assert stats.getCellTypeCountsOfCellType( cellType ) == test_case.cellTypeCounts[ + i ], f"Number of {vtkCellTypes.GetClassNameFromTypeId(cellType)} cells is expected to be {test_case.cellTypeCounts[i]}" + if stats.getCellTypeCountsOfCellType( cellType ) == 0: + continue + + # test metric summary + for j, metricIndex in enumerate( test_case.qualityMetrics ): + subStats: pd.Series = stats.getStatsFromMetricAndCellType( metricIndex, cellType ) + assert np.round( subStats, 2 ).tolist() == list( + test_case.metricsSummary[ j ] ), f"Stats at metric index {j} are wrong." + + fig: Figure = stats.plotSummaryFigure() + assert len( fig.get_axes() ) == 4, "Number of Axes is expected to be 4." diff --git a/geos-mesh/tests/test_QualityMetricSummary.py b/geos-mesh/tests/test_QualityMetricSummary.py new file mode 100644 index 00000000..8e6d7716 --- /dev/null +++ b/geos-mesh/tests/test_QualityMetricSummary.py @@ -0,0 +1,228 @@ +# 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 numpy as np +import numpy.typing as npt +import pandas as pd +import random as rd +import pytest +from typing import Iterator +from matplotlib.figure import Figure + +from geos.mesh.model.QualityMetricSummary import QualityMetricSummary, StatTypes + +from vtkmodules.vtkFiltersVerdict import vtkMeshQuality +from vtkmodules.vtkCommonDataModel import ( VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_HEXAHEDRON, VTK_WEDGE, + VTK_POLYGON, VTK_POLYHEDRON ) + +# inputs +statTypes: tuple[ StatTypes, ...] = tuple( StatTypes ) +size: int = 100 +nbNan: int = 5 +expectedStatTypes_all: tuple[ float | int, ...] = ( 0.5957, 0.6508, -2.2055, 1.6259, size - nbNan ) + +cellType_all: tuple[ int, ...] = ( VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON, + VTK_POLYGON, VTK_POLYHEDRON ) +cellMetricIndexes_all: tuple[ int, ...] = ( + ( vtkMeshQuality.QualityMeasureTypes.EDGE_RATIO, vtkMeshQuality.QualityMeasureTypes.AREA, + vtkMeshQuality.QualityMeasureTypes.CONDITION ), + ( vtkMeshQuality.QualityMeasureTypes.MED_ASPECT_FROBENIUS, vtkMeshQuality.QualityMeasureTypes.AREA, + vtkMeshQuality.QualityMeasureTypes.CONDITION ), + ( vtkMeshQuality.QualityMeasureTypes.SHAPE, vtkMeshQuality.QualityMeasureTypes.VOLUME, + vtkMeshQuality.QualityMeasureTypes.ASPECT_GAMMA ), + ( vtkMeshQuality.QualityMeasureTypes.SHAPE, vtkMeshQuality.QualityMeasureTypes.VOLUME, + vtkMeshQuality.QualityMeasureTypes.JACOBIAN, 100 ), + ( vtkMeshQuality.QualityMeasureTypes.SHAPE, vtkMeshQuality.QualityMeasureTypes.VOLUME, + vtkMeshQuality.QualityMeasureTypes.MEAN_ASPECT_FROBENIUS, 100 ), + ( vtkMeshQuality.QualityMeasureTypes.SHAPE, vtkMeshQuality.QualityMeasureTypes.VOLUME, + vtkMeshQuality.QualityMeasureTypes.ODDY, 100 ), + ( vtkMeshQuality.QualityMeasureTypes.SHAPE, vtkMeshQuality.QualityMeasureTypes.AREA ), + ( vtkMeshQuality.QualityMeasureTypes.SHAPE, vtkMeshQuality.QualityMeasureTypes.VOLUME ), +) + +metricIndexesFail_all: tuple[ int, ...] = ( + ( vtkMeshQuality.QualityMeasureTypes.MED_ASPECT_FROBENIUS, vtkMeshQuality.QualityMeasureTypes.VOLUME, 200 ), + ( vtkMeshQuality.QualityMeasureTypes.NORMALIZED_INRADIUS, vtkMeshQuality.QualityMeasureTypes.VOLUME, 200 ), + ( vtkMeshQuality.QualityMeasureTypes.ODDY, vtkMeshQuality.QualityMeasureTypes.AREA, 200 ), + ( vtkMeshQuality.QualityMeasureTypes.ASPECT_RATIO, vtkMeshQuality.QualityMeasureTypes.AREA, 200 ), + ( vtkMeshQuality.QualityMeasureTypes.ASPECT_RATIO, vtkMeshQuality.QualityMeasureTypes.AREA, 200 ), + ( vtkMeshQuality.QualityMeasureTypes.ASPECT_RATIO, vtkMeshQuality.QualityMeasureTypes.AREA, 200 ), + ( vtkMeshQuality.QualityMeasureTypes.SQUISH_INDEX, vtkMeshQuality.QualityMeasureTypes.VOLUME, 200 ), + ( vtkMeshQuality.QualityMeasureTypes.ASPECT_RATIO, vtkMeshQuality.QualityMeasureTypes.AREA, 100, 200 ), +) + + +@dataclass( frozen=True ) +class TestCase_StatsType: + """Test case.""" + __test__ = False + statType: StatTypes + values: tuple[ float, ...] + expected: int | float + + +def __generate_test_data_for_StatsType() -> Iterator[ TestCase_StatsType ]: + """Generate test cases for StatsType. + + Yields: + Iterator[ TestCase ]: Iterator on test cases for StatsType + """ + rd.seed( 10 ) + np.random.seed( 10 ) + for statType, expected in zip( statTypes, expectedStatTypes_all, strict=True ): + loc: float = rd.random() + scale: float = rd.random() + values: npt.NDArray[ np.float64 ] = np.random.normal( loc, scale, size ) + # insert nan values + for _ in range( nbNan ): + index = rd.randint( 0, size ) + values[ index ] = np.nan + yield TestCase_StatsType( statType, values, expected ) + + +ids = [ statType.getString() for statType in statTypes ] + + +@pytest.mark.parametrize( "test_case", __generate_test_data_for_StatsType(), ids=ids ) +def test_StatsType_compute( test_case: TestCase_StatsType ) -> None: + """Test of StatsType compute method. + + Args: + test_case (TestCase_StatsType): Test case + """ + obs: int | float = test_case.statType.compute( test_case.values ) + assert abs( + obs - test_case.expected + ) < 1e-4, f"Observed value is {round(obs, 4)} whereas expected value is {test_case.expected} for {test_case.statType.getString()}." + + +@dataclass( frozen=True ) +class TestCase: + """Test case.""" + __test__ = False + cellMetricIndexes: tuple[ int, ...] + otherMetricIndexes: tuple[ int, ...] + cellType: int + statTypes: StatTypes + values: tuple[ float, float, float, float, int ] + + +def __generate_test_data() -> Iterator[ TestCase ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: Iterator on test cases + """ + rd.seed( 10 ) + for cellMetricIndexes, cellType in zip( cellMetricIndexes_all, cellType_all, strict=True ): + values: tuple[ float, float, float, float, int ] = ( + rd.uniform( 0.0, 5.0 ), + rd.uniform( 0.0, 5.0 ), + rd.uniform( 0.0, 5.0 ), + rd.uniform( 0.0, 5.0 ), + rd.randint( 1, 100000 ), + ) + otherMetricIndexes: tuple[ int, ...] = ( 200, ) + yield TestCase( cellMetricIndexes, otherMetricIndexes, cellType, statTypes, values ) + + +def __generate_failed_test_data() -> Iterator[ TestCase ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: Iterator on test cases + """ + for metricIndexes, cellType in zip( metricIndexesFail_all, cellType_all, strict=True ): + otherMetricIndexes: tuple = () + values: tuple[ float, float, float, float, int ] = ( 0., 0., 0., 0., 0., 0 ) + yield TestCase( metricIndexes, otherMetricIndexes, cellType, statTypes, values ) + + +def test_QualityMetricSummary_init() -> None: + """Test of QualityMetricSummary init method.""" + stats: QualityMetricSummary = QualityMetricSummary() + assert stats.getAllCellStats() is not None, "Stats member is undefined." + assert ( stats.getAllCellStats().shape[ 0 ] + == 5 ) and stats.getAllCellStats().shape[ 1 ] == 122, "Stats shape is wrong." + + +@pytest.mark.parametrize( "test_case", __generate_failed_test_data() ) +def test_QualityMetricSummary_setter( test_case: TestCase ) -> None: + """Test of setStatsToMetricAndCellType method for IndexError. + + Args: + test_case (TestCase): Test case + """ + stats: QualityMetricSummary = QualityMetricSummary() + statType: StatTypes = StatTypes.COUNT + val: float = 1.0 + for cellMetricIndex in test_case.cellMetricIndexes: + with pytest.raises( IndexError ) as pytest_wrapped_e: + stats.setCellStatValueFromStatMetricAndCellType( cellMetricIndex, test_case.cellType, statType, val ) + assert pytest_wrapped_e.type is IndexError + + +@pytest.mark.parametrize( "test_case", __generate_test_data() ) +def test_QualityMetricSummary_setterGetter( test_case: TestCase ) -> None: + """Test of setter and getter methods. + + Args: + test_case (TestCase): Test case + """ + stats: QualityMetricSummary = QualityMetricSummary() + for cellMetricIndex in test_case.cellMetricIndexes: + for statType, value in zip( test_case.statTypes, test_case.values, strict=True ): + stats.setCellStatValueFromStatMetricAndCellType( cellMetricIndex, test_case.cellType, statType, value ) + + for otherMetricIndex in test_case.otherMetricIndexes: + for statType, value in zip( test_case.statTypes, test_case.values, strict=True ): + stats.setOtherStatValueFromMetric( otherMetricIndex, statType, value ) + + assert np.any( stats.getAllCellStats().to_numpy() > 0 ), "Stats values were not correctly set." + for cellMetricIndex in test_case.cellMetricIndexes: + for statType, val in zip( test_case.statTypes, test_case.values, strict=True ): + subSet: pd.DataFrame = stats.getCellStatsFromCellType( test_case.cellType ) + assert subSet[ cellMetricIndex ][ statType.getIndex( + ) ] == val, f"Stats at ({cellMetricIndex}, {test_case.cellType}, {statType}) from getCellStatsFromCellType is expected to be equal to {val}." + + subSet2: pd.DataFrame = stats.getStatsFromMetric( cellMetricIndex ) + assert subSet2[ test_case.cellType ][ statType.getIndex( + ) ] == val, f"Stats at ({cellMetricIndex}, {test_case.cellType}, {statType}) from getStatsFromMetric is expected to be equal to {val}." + + subSet3: pd.Series = stats.getStatsFromMetricAndCellType( cellMetricIndex, test_case.cellType ) + assert subSet3[ statType.getIndex( + ) ] == val, f"Stats at ({cellMetricIndex}, {test_case.cellType}, {statType}) from getStatsFromMetricAndCellType is expected to be equal to {val}." + + valObs: float = stats.getCellStatValueFromStatMetricAndCellType( cellMetricIndex, test_case.cellType, + statType ) + assert valObs == val, f"Stats at ({cellMetricIndex}, {test_case.cellType}, {statType}) from getStatValueFromMetricAndCellType is expected to be equal to {val}." + + for cellMetricIndex in test_case.otherMetricIndexes: + for statType, val in zip( test_case.statTypes, test_case.values, strict=True ): + subSet4: pd.DataFrame = stats.getStatsFromMetric( cellMetricIndex ) + assert subSet4[ statType.getIndex( + ) ] == val, f"Stats at ({cellMetricIndex}, {statType}) from getStatsFromMetric is expected to be equal to {val}." + + +@pytest.mark.parametrize( "test_case", __generate_test_data() ) +def test_QualityMetricSummary_plotSummaryFigure( test_case: TestCase ) -> None: + """Test of plotSummaryFigure method. + + Args: + test_case (TestCase): Test case + """ + stats: QualityMetricSummary = QualityMetricSummary() + for cellMetricIndex in test_case.cellMetricIndexes: + for statType, value in zip( test_case.statTypes, test_case.values, strict=True ): + stats.setCellStatValueFromStatMetricAndCellType( cellMetricIndex, test_case.cellType, statType, value ) + + for otherMetricIndex in test_case.otherMetricIndexes: + for statType, value in zip( test_case.statTypes, test_case.values, strict=True ): + stats.setOtherStatValueFromMetric( otherMetricIndex, statType, value ) + + fig: Figure = stats.plotSummaryFigure() + assert fig is not None, "Figure is undefined" + # metrics + counts + nbAxesExp: int = len( test_case.cellMetricIndexes ) + len( test_case.otherMetricIndexes ) + 1 + assert len( fig.get_axes() ) == nbAxesExp, f"Number of Axes is expected to be {nbAxesExp}." \ No newline at end of file diff --git a/geos-mesh/tests/test_meshQualityHelpers.py b/geos-mesh/tests/test_meshQualityHelpers.py new file mode 100644 index 00000000..f961b32b --- /dev/null +++ b/geos-mesh/tests/test_meshQualityHelpers.py @@ -0,0 +1,364 @@ +# 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 vtkmodules.vtkFiltersVerdict import vtkMeshQuality +from vtkmodules.vtkCommonDataModel import ( VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON, + VTK_POLYGON, VTK_POLYHEDRON ) + +from geos.mesh.processing.meshQualityMetricHelpers import ( + VtkCellQualityMetricEnum, + CellQualityMetricAdditionalEnum, + QualityRange, + getCellQualityMeasureFromCellType, + getTriangleQualityMeasure, + getQuadQualityMeasure, + getTetQualityMeasure, + getPyramidQualityMeasure, + getWedgeQualityMeasure, + getHexQualityMeasure, + getCommonPolygonQualityMeasure, + getCommonPolyhedraQualityMeasure, + getQualityMeasureNameFromIndex, + getQualityMeasureIndexFromName, + getAllCellTypesExtended, + getAllCellTypes, +) + +triangleQualityMeasureExp: set[ int ] = { + vtkMeshQuality.QualityMeasureTypes.AREA, vtkMeshQuality.QualityMeasureTypes.ASPECT_RATIO, + vtkMeshQuality.QualityMeasureTypes.ASPECT_FROBENIUS, vtkMeshQuality.QualityMeasureTypes.CONDITION, + vtkMeshQuality.QualityMeasureTypes.DISTORTION, vtkMeshQuality.QualityMeasureTypes.EDGE_RATIO, + vtkMeshQuality.QualityMeasureTypes.EQUIANGLE_SKEW, vtkMeshQuality.QualityMeasureTypes.MAX_ANGLE, + vtkMeshQuality.QualityMeasureTypes.MIN_ANGLE, vtkMeshQuality.QualityMeasureTypes.NORMALIZED_INRADIUS, + vtkMeshQuality.QualityMeasureTypes.RADIUS_RATIO, vtkMeshQuality.QualityMeasureTypes.RELATIVE_SIZE_SQUARED, + vtkMeshQuality.QualityMeasureTypes.SCALED_JACOBIAN, vtkMeshQuality.QualityMeasureTypes.SHAPE, + vtkMeshQuality.QualityMeasureTypes.SHAPE_AND_SIZE +} + +quadQualityMeasureExp: set[ int ] = { + vtkMeshQuality.QualityMeasureTypes.AREA, + vtkMeshQuality.QualityMeasureTypes.ASPECT_RATIO, + vtkMeshQuality.QualityMeasureTypes.CONDITION, + vtkMeshQuality.QualityMeasureTypes.DISTORTION, + vtkMeshQuality.QualityMeasureTypes.EDGE_RATIO, + vtkMeshQuality.QualityMeasureTypes.EQUIANGLE_SKEW, + vtkMeshQuality.QualityMeasureTypes.JACOBIAN, + vtkMeshQuality.QualityMeasureTypes.MED_ASPECT_FROBENIUS, + vtkMeshQuality.QualityMeasureTypes.MAX_ASPECT_FROBENIUS, + vtkMeshQuality.QualityMeasureTypes.MAX_EDGE_RATIO, + vtkMeshQuality.QualityMeasureTypes.MAX_ANGLE, + vtkMeshQuality.QualityMeasureTypes.MIN_ANGLE, + vtkMeshQuality.QualityMeasureTypes.ODDY, + vtkMeshQuality.QualityMeasureTypes.RADIUS_RATIO, + vtkMeshQuality.QualityMeasureTypes.RELATIVE_SIZE_SQUARED, + vtkMeshQuality.QualityMeasureTypes.SCALED_JACOBIAN, + vtkMeshQuality.QualityMeasureTypes.SHAPE, + vtkMeshQuality.QualityMeasureTypes.SHAPE_AND_SIZE, + vtkMeshQuality.QualityMeasureTypes.SHEAR, + vtkMeshQuality.QualityMeasureTypes.SHEAR_AND_SIZE, + vtkMeshQuality.QualityMeasureTypes.SKEW, + vtkMeshQuality.QualityMeasureTypes.STRETCH, + vtkMeshQuality.QualityMeasureTypes.TAPER, + vtkMeshQuality.QualityMeasureTypes.WARPAGE, +} + +tetQualityMeasureExp: set[ int ] = { + vtkMeshQuality.QualityMeasureTypes.ASPECT_FROBENIUS, + vtkMeshQuality.QualityMeasureTypes.ASPECT_GAMMA, + vtkMeshQuality.QualityMeasureTypes.ASPECT_RATIO, + vtkMeshQuality.QualityMeasureTypes.COLLAPSE_RATIO, + vtkMeshQuality.QualityMeasureTypes.CONDITION, + vtkMeshQuality.QualityMeasureTypes.DISTORTION, + vtkMeshQuality.QualityMeasureTypes.EDGE_RATIO, + vtkMeshQuality.QualityMeasureTypes.EQUIANGLE_SKEW, + vtkMeshQuality.QualityMeasureTypes.EQUIVOLUME_SKEW, + vtkMeshQuality.QualityMeasureTypes.JACOBIAN, + vtkMeshQuality.QualityMeasureTypes.MEAN_RATIO, + vtkMeshQuality.QualityMeasureTypes.MIN_ANGLE, + vtkMeshQuality.QualityMeasureTypes.NORMALIZED_INRADIUS, + vtkMeshQuality.QualityMeasureTypes.RADIUS_RATIO, + vtkMeshQuality.QualityMeasureTypes.RELATIVE_SIZE_SQUARED, + vtkMeshQuality.QualityMeasureTypes.SCALED_JACOBIAN, + vtkMeshQuality.QualityMeasureTypes.SHAPE, + vtkMeshQuality.QualityMeasureTypes.SHAPE_AND_SIZE, + vtkMeshQuality.QualityMeasureTypes.SQUISH_INDEX, + vtkMeshQuality.QualityMeasureTypes.VOLUME, +} + +pyrQualityMeasureExp: set[ int ] = { + vtkMeshQuality.QualityMeasureTypes.EQUIANGLE_SKEW, + vtkMeshQuality.QualityMeasureTypes.JACOBIAN, + vtkMeshQuality.QualityMeasureTypes.SCALED_JACOBIAN, + vtkMeshQuality.QualityMeasureTypes.SHAPE, + vtkMeshQuality.QualityMeasureTypes.VOLUME, + vtkMeshQuality.QualityMeasureTypes.SQUISH_INDEX, + CellQualityMetricAdditionalEnum.MAXIMUM_ASPECT_RATIO.metricIndex, +} + +wedgeQualityMeasureExp: set[ int ] = { + vtkMeshQuality.QualityMeasureTypes.CONDITION, + vtkMeshQuality.QualityMeasureTypes.DISTORTION, + vtkMeshQuality.QualityMeasureTypes.EDGE_RATIO, + vtkMeshQuality.QualityMeasureTypes.EQUIANGLE_SKEW, + vtkMeshQuality.QualityMeasureTypes.JACOBIAN, + vtkMeshQuality.QualityMeasureTypes.MAX_ASPECT_FROBENIUS, + vtkMeshQuality.QualityMeasureTypes.MAX_STRETCH, + vtkMeshQuality.QualityMeasureTypes.MEAN_ASPECT_FROBENIUS, + vtkMeshQuality.QualityMeasureTypes.SCALED_JACOBIAN, + vtkMeshQuality.QualityMeasureTypes.SHAPE, + vtkMeshQuality.QualityMeasureTypes.VOLUME, + vtkMeshQuality.QualityMeasureTypes.SQUISH_INDEX, + CellQualityMetricAdditionalEnum.MAXIMUM_ASPECT_RATIO.metricIndex, +} + +hexQualityMeasureExp: set[ int ] = { + vtkMeshQuality.QualityMeasureTypes.CONDITION, + vtkMeshQuality.QualityMeasureTypes.DIAGONAL, + vtkMeshQuality.QualityMeasureTypes.DIMENSION, + vtkMeshQuality.QualityMeasureTypes.DISTORTION, + vtkMeshQuality.QualityMeasureTypes.EDGE_RATIO, + vtkMeshQuality.QualityMeasureTypes.EQUIANGLE_SKEW, + vtkMeshQuality.QualityMeasureTypes.JACOBIAN, + vtkMeshQuality.QualityMeasureTypes.MAX_EDGE_RATIO, + vtkMeshQuality.QualityMeasureTypes.MAX_ASPECT_FROBENIUS, + vtkMeshQuality.QualityMeasureTypes.MED_ASPECT_FROBENIUS, + vtkMeshQuality.QualityMeasureTypes.NODAL_JACOBIAN_RATIO, + vtkMeshQuality.QualityMeasureTypes.ODDY, + vtkMeshQuality.QualityMeasureTypes.RELATIVE_SIZE_SQUARED, + vtkMeshQuality.QualityMeasureTypes.SCALED_JACOBIAN, + vtkMeshQuality.QualityMeasureTypes.SHAPE, + vtkMeshQuality.QualityMeasureTypes.SHAPE_AND_SIZE, + vtkMeshQuality.QualityMeasureTypes.SHEAR, + vtkMeshQuality.QualityMeasureTypes.SHEAR_AND_SIZE, + vtkMeshQuality.QualityMeasureTypes.SKEW, + vtkMeshQuality.QualityMeasureTypes.STRETCH, + vtkMeshQuality.QualityMeasureTypes.TAPER, + vtkMeshQuality.QualityMeasureTypes.VOLUME, + vtkMeshQuality.QualityMeasureTypes.SQUISH_INDEX, + CellQualityMetricAdditionalEnum.MAXIMUM_ASPECT_RATIO.metricIndex, +} + + +@dataclass( frozen=True ) +class TestCase: + """Test case.""" + __test__ = False + #: input mesh + qualityMetricIndex: int + #: expected points + qualityMetricName: str + + +def __generate_test_data() -> Iterator[ TestCase ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: Iterator on test cases + """ + for metric in list( VtkCellQualityMetricEnum ): + yield TestCase( metric.metricIndex, metric.metricName ) + + +def test_CellQualityMetricEnum_Order() -> None: + """Test VtkCellQualityMetricEnum ordering is correct.""" + for i, metric in enumerate( list( VtkCellQualityMetricEnum ) ): + assert metric.metricIndex == i + + +def test_CellQualityMetricEnum_QualityRange() -> None: + """Test VtkCellQualityMetricEnum.getQualityRange returns the right number of values.""" + for metric in list( VtkCellQualityMetricEnum ): + for cellType in getAllCellTypes(): + qualityRange: QualityRange = metric.getQualityRange( cellType ) + if qualityRange is not None: + assert ( len( qualityRange.fullRange ) == 2 ), "Full range length is expected to be 2" + assert ( len( qualityRange.normalRange ) == 2 ), "Normal range length is expected to be 2" + assert ( len( qualityRange.acceptableRange ) == 2 ), "Acceptable range length is expected to be 2" + + for cellType in ( VTK_POLYGON, VTK_POLYHEDRON ): + assert metric.getQualityRange( cellType ) is None, "QualityRange should be undefined." + + +ids: list[ str ] = [ metric.metricName for metric in list( VtkCellQualityMetricEnum ) ] + + +@pytest.mark.parametrize( "test_case", __generate_test_data(), ids=ids ) +def test_getQualityMeasureNameFromIndex( test_case: TestCase ) -> None: + """Test of getQualityMeasureNameFromIndex method.""" + name: str = getQualityMeasureNameFromIndex( test_case.qualityMetricIndex ) + assert name == test_case.qualityMetricName + + +@pytest.mark.parametrize( "test_case", __generate_test_data(), ids=ids ) +def test_getQualityMeasureIndexFromName( test_case: TestCase ) -> None: + """Test of getQualityMeasureIndexFromName method.""" + index: int = getQualityMeasureIndexFromName( test_case.qualityMetricName ) + assert index == test_case.qualityMetricIndex + + +def test_getQualityMeasureFromCellType_exception() -> None: + """Test of supported cell type from getCellQualityMeasureFromCellType method.""" + for i in range( 20 ): + if i in ( VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON, VTK_POLYGON, + VTK_POLYHEDRON ): + assert len( getCellQualityMeasureFromCellType( i ) ) > 0 + else: + with pytest.raises( ValueError ) as pytest_wrapped_e: + getCellQualityMeasureFromCellType( i ) + assert pytest_wrapped_e.type is ValueError + + +def test_getTriangleQualityMeasure() -> None: + """Test of getTriangleQualityMeasure method.""" + obs: set[ int ] = getTriangleQualityMeasure() + diffAdditional: set[ int ] = obs.difference( triangleQualityMeasureExp ) + diffMissing: set[ int ] = triangleQualityMeasureExp.difference( obs ) + messAdditional: str = f"{len(diffAdditional)} additional elements" + if len( diffAdditional ) > 0: + messAdditional += f" including {[getQualityMeasureNameFromIndex(index) for index in diffAdditional]}" + messMissing: str = f" and {len(diffMissing)} missing elements" + if len( diffMissing ): + messMissing += f" including {[getQualityMeasureNameFromIndex(index) for index in diffMissing]}" + expNames: list[ str ] = [ getQualityMeasureNameFromIndex( index ) for index in sorted( triangleQualityMeasureExp ) ] + assert sorted( + obs ) == sorted( triangleQualityMeasureExp + ), f"Expected metrics are {expNames}. Observed metrics contains " + messAdditional + messMissing + + +def test_getQuadQualityMeasure() -> None: + """Test of getQuadQualityMeasure method.""" + obs: set[ int ] = getQuadQualityMeasure() + diffAdditional: set[ int ] = obs.difference( quadQualityMeasureExp ) + diffMissing: set[ int ] = quadQualityMeasureExp.difference( obs ) + messAdditional: str = f"{len(diffAdditional)} additional elements" + if len( diffAdditional ) > 0: + messAdditional += f" including {[getQualityMeasureNameFromIndex(index) for index in diffAdditional]}" + messMissing: str = f" and {len(diffMissing)} missing elements" + if len( diffMissing ): + messMissing += f" including {[getQualityMeasureNameFromIndex(index) for index in diffMissing]}" + expNames: list[ str ] = [ getQualityMeasureNameFromIndex( index ) for index in sorted( quadQualityMeasureExp ) ] + assert sorted( + obs ) == sorted( quadQualityMeasureExp + ), f"Expected metrics are {expNames}. Observed metrics contains " + messAdditional + messMissing + + +def test_getCommonPolygonQualityMeasure() -> None: + """Test of getCommonPolygonQualityMeasure method.""" + obs: set[ int ] = getCommonPolygonQualityMeasure() + exp: set[ int ] = quadQualityMeasureExp.intersection( triangleQualityMeasureExp ) + diffAdditional: set[ int ] = obs.difference( exp ) + diffMissing: set[ int ] = exp.difference( obs ) + messAdditional: str = f"{len(diffAdditional)} additional elements" + if len( diffAdditional ) > 0: + messAdditional += f" including {[getQualityMeasureNameFromIndex(index) for index in diffAdditional]}" + messMissing: str = f" and {len(diffMissing)} missing elements" + if len( diffMissing ): + messMissing += f" including {[getQualityMeasureNameFromIndex(index) for index in diffMissing]}" + expNames: list[ str ] = [ getQualityMeasureNameFromIndex( index ) for index in sorted( exp ) ] + assert sorted( obs ) == sorted( + exp ), f"Expected metrics are {expNames}. Observed metrics contains " + messAdditional + messMissing + + +def test_getTetQualityMeasure() -> None: + """Test of getTetQualityMeasure method.""" + obs: set[ int ] = getTetQualityMeasure() + diffAdditional: set[ int ] = obs.difference( tetQualityMeasureExp ) + diffMissing: set[ int ] = tetQualityMeasureExp.difference( obs ) + messAdditional: str = f"{len(diffAdditional)} additional elements" + if len( diffAdditional ) > 0: + messAdditional += f" including {[getQualityMeasureNameFromIndex(index) for index in diffAdditional]}" + messMissing: str = f" and {len(diffMissing)} missing elements" + if len( diffMissing ): + messMissing += f" including {[getQualityMeasureNameFromIndex(index) for index in diffMissing]}" + expNames: list[ str ] = [ getQualityMeasureNameFromIndex( index ) for index in sorted( tetQualityMeasureExp ) ] + assert sorted( + obs ) == sorted( tetQualityMeasureExp + ), f"Expected metrics are {expNames}. Observed metrics contains " + messAdditional + messMissing + + +def test_getPyramidQualityMeasure() -> None: + """Test of getPyramidQualityMeasure method.""" + obs: set[ int ] = getPyramidQualityMeasure() + diffAdditional: set[ int ] = obs.difference( pyrQualityMeasureExp ) + diffMissing: set[ int ] = pyrQualityMeasureExp.difference( obs ) + messAdditional: str = f"{len(diffAdditional)} additional elements" + if len( diffAdditional ) > 0: + messAdditional += f" including {[getQualityMeasureNameFromIndex(index) for index in diffAdditional]}" + messMissing: str = f" and {len(diffMissing)} missing elements" + if len( diffMissing ): + messMissing += f" including {[getQualityMeasureNameFromIndex(index) for index in diffMissing]}" + expNames: list[ str ] = [ getQualityMeasureNameFromIndex( index ) for index in sorted( pyrQualityMeasureExp ) ] + assert sorted( + obs ) == sorted( pyrQualityMeasureExp + ), f"Expected metrics are {expNames}. Observed metrics contains " + messAdditional + messMissing + + +def test_getWedgeQualityMeasure() -> None: + """Test of getWedgeQualityMeasure method.""" + obs: set[ int ] = getWedgeQualityMeasure() + diffAdditional: set[ int ] = obs.difference( wedgeQualityMeasureExp ) + diffMissing: set[ int ] = wedgeQualityMeasureExp.difference( obs ) + messAdditional: str = f"{len(diffAdditional)} additional elements" + if len( diffAdditional ) > 0: + messAdditional += f" including {[getQualityMeasureNameFromIndex(index) for index in diffAdditional]}" + messMissing: str = f" and {len(diffMissing)} missing elements" + if len( diffMissing ): + messMissing += f" including {[getQualityMeasureNameFromIndex(index) for index in diffMissing]}" + expNames: list[ str ] = [ getQualityMeasureNameFromIndex( index ) for index in sorted( wedgeQualityMeasureExp ) ] + assert sorted( + obs ) == sorted( wedgeQualityMeasureExp + ), f"Expected metrics are {expNames}. Observed metrics contains " + messAdditional + messMissing + + +def test_getHexQualityMeasure() -> None: + """Test of getHexQualityMeasure method.""" + obs: set[ int ] = getHexQualityMeasure() + diffAdditional: set[ int ] = obs.difference( hexQualityMeasureExp ) + diffMissing: set[ int ] = hexQualityMeasureExp.difference( obs ) + messAdditional: str = f"{len(diffAdditional)} additional elements" + if len( diffAdditional ) > 0: + messAdditional += f" including {[getQualityMeasureNameFromIndex(index) for index in diffAdditional]}" + messMissing: str = f" and {len(diffMissing)} missing elements" + if len( diffMissing ): + messMissing += f" including {[getQualityMeasureNameFromIndex(index) for index in diffMissing]}" + expNames: list[ str ] = [ getQualityMeasureNameFromIndex( index ) for index in sorted( hexQualityMeasureExp ) ] + assert sorted( + obs ) == sorted( hexQualityMeasureExp + ), f"Expected metrics are {expNames}. Observed metrics contains " + messAdditional + messMissing + + +def test_getCommonPolyhedraQualityMeasure() -> None: + """Test of getCommonPolyhedraQualityMeasure method.""" + obs: set[ int ] = getCommonPolyhedraQualityMeasure() + exp: set[ int ] = tetQualityMeasureExp.intersection( pyrQualityMeasureExp ).intersection( + wedgeQualityMeasureExp ).intersection( hexQualityMeasureExp ) + diffAdditional: set[ int ] = obs.difference( exp ) + diffMissing: set[ int ] = exp.difference( obs ) + messAdditional: str = f"{len(diffAdditional)} additional elements" + if len( diffAdditional ) > 0: + messAdditional += f" including {[getQualityMeasureNameFromIndex(index) for index in diffAdditional]}" + messMissing: str = f" and {len(diffMissing)} missing elements" + if len( diffMissing ): + messMissing += f" including {[getQualityMeasureNameFromIndex(index) for index in diffMissing]}" + expNames: list[ str ] = [ getQualityMeasureNameFromIndex( index ) for index in sorted( exp ) ] + assert sorted( obs ) == sorted( + exp ), f"Expected metrics are {expNames}. Observed metrics contains " + messAdditional + messMissing + + +def test_getAllCellTypesExtended() -> None: + """Test of getAllCellTypesExtended method.""" + cellTypesExp: list[ int ] = [ + VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON, VTK_POLYGON, VTK_POLYHEDRON + ] + assert cellTypesExp == getAllCellTypesExtended(), "Cell types differ." + + +def test_getAllCellTypes() -> None: + """Test of getAllCellTypes method.""" + cellTypesExp: list[ int ] = [ VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON ] + assert cellTypesExp == getAllCellTypes(), "Cell types differ." \ No newline at end of file diff --git a/geos-pv/src/PVplugins/PVCellTypeCounterEnhanced.py b/geos-pv/src/PVplugins/PVCellTypeCounterEnhanced.py new file mode 100644 index 00000000..25453fc9 --- /dev/null +++ b/geos-pv/src/PVplugins/PVCellTypeCounterEnhanced.py @@ -0,0 +1,158 @@ +# 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.CellTypeCounterEnhanced import CellTypeCounterEnhanced +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>PVCellTypeCounterEnhanced. +* Select the input mesh. +* Apply the filter. + +""" + + +@smproxy.filter( name="PVCellTypeCounterEnhanced", label="Cell Type Counter" ) +@smhint.xml( '' ) +@smproperty.input( name="Input", port_index=0 ) +@smdomain.datatype( + dataTypes=[ "vtkUnstructuredGrid" ], + composite_data_supported=True, +) +class PVCellTypeCounterEnhanced( 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: CellTypeCounterEnhanced = CellTypeCounterEnhanced() + filter.SetInputDataObject( inputMesh ) + filter.Update() + outputTable.ShallowCopy( filter.GetOutputDataObject( 0 ) ) + + # print counts in Output Messages view + counts: CellTypeCounts = filter.GetCellTypeCountsObject() + + 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 \ No newline at end of file diff --git a/geos-pv/src/PVplugins/PVSurfaceMeshQualityEnhanced.py b/geos-pv/src/PVplugins/PVSurfaceMeshQualityEnhanced.py new file mode 100644 index 00000000..7d5b2d48 --- /dev/null +++ b/geos-pv/src/PVplugins/PVSurfaceMeshQualityEnhanced.py @@ -0,0 +1,277 @@ +# 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, Optional + +from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] + VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, +) + +from vtkmodules.vtkCommonCore import ( + vtkInformation, + vtkInformationVector, + vtkDataArraySelection, +) +from vtkmodules.vtkCommonDataModel import ( + vtkUnstructuredGrid, ) + +# 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.model.QualityMetricSummary import QualityMetricSummary +from geos.mesh.stats.MeshQualityEnhanced import MeshQualityEnhanced + +from geos.mesh.processing.meshQualityMetricHelpers import ( + getQualityMetricsOther, + getQualityMeasureNameFromIndex, + getQualityMeasureIndexFromName, + getQuadQualityMeasure, + getTriangleQualityMeasure, + getCommonPolygonQualityMeasure, +) +from geos.pv.utils.checkboxFunction import ( # type: ignore[attr-defined] + createModifiedCallback, ) +from geos.pv.utils.paraviewTreatments import getArrayChoices + +__doc__ = """ +Compute mesh quality metrics on surface meshes. + +To use it: + +* Load the module in Paraview: Tools>Manage Plugins...>Load new>PVSurfaceMeshQualityEnhanced. +* Select the input mesh. +* Select the metric to compute +* Apply the filter. + +""" + + +@smproxy.filter( name="PVSurfaceMeshQualityEnhanced", label="Surface Mesh Quality Enhanced" ) +@smhint.xml( '' ) +@smproperty.input( name="Input", port_index=0 ) +@smdomain.datatype( + dataTypes=[ "vtkUnstructuredGrid" ], + composite_data_supported=True, +) +class PVSurfaceMeshQualityEnhanced( VTKPythonAlgorithmBase ): + + def __init__( self: Self ) -> None: + """Merge collocated points.""" + super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkUnstructuredGrid" ) + + self._filename: Optional[ str ] = None + self._saveToFile: bool = True + self._blockIndex: int = 0 + # used to concatenate results if vtkMultiBlockDataSet + self._metricsAll: list[ float ] = [] + self._commonMeshQualityMetric: vtkDataArraySelection = vtkDataArraySelection() + self._commonCellQualityMetric: vtkDataArraySelection = vtkDataArraySelection() + self._triangleQualityMetric: vtkDataArraySelection = vtkDataArraySelection() + self._quadsQualityMetric: vtkDataArraySelection = vtkDataArraySelection() + self._initQualityMetricSelection() + + def _initQualityMetricSelection( self: Self ) -> None: + self._commonCellQualityMetric.RemoveAllArrays() + self._commonCellQualityMetric.AddObserver( + "ModifiedEvent", # type: ignore[arg-type] + createModifiedCallback( self ) ) + commonCellMetrics: set[ int ] = getCommonPolygonQualityMeasure() + for measure in commonCellMetrics: + self._commonCellQualityMetric.AddArray( getQualityMeasureNameFromIndex( measure ) ) + + self._triangleQualityMetric.RemoveAllArrays() + self._triangleQualityMetric.AddObserver( + "ModifiedEvent", # type: ignore[arg-type] + createModifiedCallback( self ) ) + for measure in getTriangleQualityMeasure().difference( commonCellMetrics ): + self._triangleQualityMetric.AddArray( getQualityMeasureNameFromIndex( measure ) ) + + self._quadsQualityMetric.RemoveAllArrays() + self._quadsQualityMetric.AddObserver( + "ModifiedEvent", # type: ignore[arg-type] + createModifiedCallback( self ) ) + for measure in getQuadQualityMeasure().difference( commonCellMetrics ): + self._quadsQualityMetric.AddArray( getQualityMeasureNameFromIndex( measure ) ) + + otherMetrics: set[ int ] = getQualityMetricsOther() + for measure in otherMetrics: + self._commonMeshQualityMetric.AddArray( getQualityMeasureNameFromIndex( measure ) ) + + @smproperty.dataarrayselection( name="CommonCellQualityMetric" ) + def a01SetCommonMetrics( self: Self ) -> vtkDataArraySelection: + """Set polygon quality metrics selection.""" + return self._commonCellQualityMetric + + @smproperty.dataarrayselection( name="TriangleSpecificQualityMetric" ) + def a02SetTriangleMetrics( self: Self ) -> vtkDataArraySelection: + """Set triangle quality metrics selection.""" + return self._triangleQualityMetric + + @smproperty.dataarrayselection( name="QuadSpecificQualityMetric" ) + def a03sSetQuadMetrics( self: Self ) -> vtkDataArraySelection: + """Set quad quality metrics selection.""" + return self._quadsQualityMetric + + @smproperty.dataarrayselection( name="OtherMeshQualityMetric" ) + def a04SetOtherMeshMetrics( self: Self ) -> vtkDataArraySelection: + """Set other mesh quality metrics selection.""" + return self._commonMeshQualityMetric + + @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 b01SetSaveToFile( 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 b02SetFileName( 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 b03GroupAdvancedOutputParameters( self: Self ) -> None: + """Organize groups.""" + self.Modified() + + def Modified( self: Self ) -> None: + """Overload Modified method to reset _blockIndex.""" + self._blockIndex = 0 + super().Modified() + + def RequestDataObject( + self: Self, + request: vtkInformation, + inInfoVec: list[ vtkInformationVector ], + outInfoVec: vtkInformationVector, + ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestDataObject. + + Args: + request (vtkInformation): Request + inInfoVec (list[vtkInformationVector]): Input objects + outInfoVec (vtkInformationVector): Output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + inData = self.GetInputData( inInfoVec, 0, 0 ) + outData = self.GetOutputData( outInfoVec, 0 ) + assert inData is not None + if outData is None or ( not outData.IsA( inData.GetClassName() ) ): + outData = inData.NewInstance() + outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) + return super().RequestDataObject( request, inInfoVec, outInfoVec ) + + def _getQualityMetricsToUse( self: Self, selection: vtkDataArraySelection ) -> set[ int ]: + """Get mesh quality metric indexes from user selection. + + Returns: + list[int]: List of quality metric indexes + """ + metricsNames: set[ str ] = getArrayChoices( selection ) + return { getQualityMeasureIndexFromName( name ) for name in metricsNames } + + 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: vtkUnstructuredGrid = self.GetInputData( inInfoVec, 0, 0 ) + outputMesh: vtkUnstructuredGrid = vtkUnstructuredGrid.GetData( outInfoVec, 0 ) + assert inputMesh is not None, "Input server mesh is null." + assert outputMesh is not None, "Output pipeline is null." + + triangleMetrics: set[ int ] = self._getQualityMetricsToUse( self._commonCellQualityMetric ).union( + self._getQualityMetricsToUse( self._triangleQualityMetric ) ) + quadMetrics: set[ int ] = self._getQualityMetricsToUse( self._commonCellQualityMetric ).union( + self._getQualityMetricsToUse( self._quadsQualityMetric ) ) + otherMetrics: set[ int ] = self._getQualityMetricsToUse( self._commonMeshQualityMetric ) + filter: MeshQualityEnhanced = MeshQualityEnhanced() + filter.SetInputDataObject( inputMesh ) + filter.SetCellQualityMetrics( triangleMetrics=triangleMetrics, quadMetrics=quadMetrics ) + filter.SetOtherMeshQualityMetrics( otherMetrics ) + filter.Update() + + outputMesh.ShallowCopy( filter.GetOutputDataObject( 0 ) ) + + # save to file if asked + if self._saveToFile: + stats: QualityMetricSummary = filter.GetQualityMetricSummary() + self.saveFile( stats ) + self._blockIndex += 1 + return 1 + + def saveFile( self: Self, stats: QualityMetricSummary ) -> None: + """Export mesh quality metric summary file.""" + try: + if self._filename is None: + print( "Mesh quality summary report file path is undefined." ) + return + # add index for multiblock meshes + index: int = self._filename.rfind( '.' ) + filename: str = self._filename[ :index ] + f"_{self._blockIndex}" + self._filename[ index: ] + fig = stats.plotSummaryFigure() + fig.savefig( filename, dpi=150 ) + print( f"File {filename} was successfully written." ) + except Exception as e: + print( "Error while exporting the file due to:" ) + print( str( e ) ) \ No newline at end of file diff --git a/geos-pv/src/PVplugins/PVVolumeMeshQualityEnhanced.py b/geos-pv/src/PVplugins/PVVolumeMeshQualityEnhanced.py new file mode 100644 index 00000000..43660c56 --- /dev/null +++ b/geos-pv/src/PVplugins/PVVolumeMeshQualityEnhanced.py @@ -0,0 +1,306 @@ +# 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, Optional + +from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] + VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, +) + +from vtkmodules.vtkCommonCore import ( + vtkInformation, + vtkInformationVector, + vtkDataArraySelection, +) +from vtkmodules.vtkCommonDataModel import ( + vtkUnstructuredGrid, ) + +# 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.model.QualityMetricSummary import QualityMetricSummary +from geos.mesh.stats.MeshQualityEnhanced import MeshQualityEnhanced + +from geos.mesh.processing.meshQualityMetricHelpers import ( + getQualityMetricsOther, + getQualityMeasureNameFromIndex, + getQualityMeasureIndexFromName, + getTetQualityMeasure, + getPyramidQualityMeasure, + getWedgeQualityMeasure, + getHexQualityMeasure, + getCommonPolyhedraQualityMeasure, +) +from geos.pv.utils.checkboxFunction import ( # type: ignore[attr-defined] + createModifiedCallback, ) +from geos.pv.utils.paraviewTreatments import getArrayChoices + +__doc__ = """ +Compute mesh quality metrics on volume meshes. + +To use it: + +* Load the module in Paraview: Tools>Manage Plugins...>Load new>PVVolumeMeshQualityEnhanced. +* Select the input mesh. +* Select the metric to compute +* Apply the filter. + +""" + + +@smproxy.filter( name="PVVolumeMeshQualityEnhanced", label="Volume Mesh Quality Enhanced" ) +@smhint.xml( '' ) +@smproperty.input( name="Input", port_index=0 ) +@smdomain.datatype( + dataTypes=[ "vtkUnstructuredGrid" ], + composite_data_supported=True, +) +class PVVolumeMeshQualityEnhanced( VTKPythonAlgorithmBase ): + + def __init__( self: Self ) -> None: + """Merge collocated points.""" + super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkUnstructuredGrid" ) + + self._filename: Optional[ str ] = None + self._saveToFile: bool = True + self._blockIndex: int = 0 + # used to concatenate results if vtkMultiBlockDataSet + self._metricsAll: list[ float ] = [] + self._commonMeshQualityMetric: vtkDataArraySelection = vtkDataArraySelection() + self._commonCellQualityMetric: vtkDataArraySelection = vtkDataArraySelection() + self._tetQualityMetric: vtkDataArraySelection = vtkDataArraySelection() + self._PyrQualityMetric: vtkDataArraySelection = vtkDataArraySelection() + self._WedgeQualityMetric: vtkDataArraySelection = vtkDataArraySelection() + self._HexQualityMetric: vtkDataArraySelection = vtkDataArraySelection() + self._initQualityMetricSelection() + + def _initQualityMetricSelection( self: Self ) -> None: + self._commonCellQualityMetric.RemoveAllArrays() + self._commonCellQualityMetric.AddObserver( + "ModifiedEvent", # type: ignore[arg-type] + createModifiedCallback( self ) ) + commonCellMetrics: set[ int ] = getCommonPolyhedraQualityMeasure() + for measure in commonCellMetrics: + self._commonCellQualityMetric.AddArray( getQualityMeasureNameFromIndex( measure ) ) + + self._tetQualityMetric.RemoveAllArrays() + self._tetQualityMetric.AddObserver( "ModifiedEvent", createModifiedCallback( self ) ) # type: ignore[arg-type] + for measure in getTetQualityMeasure().difference( commonCellMetrics ): + self._tetQualityMetric.AddArray( getQualityMeasureNameFromIndex( measure ) ) + + self._PyrQualityMetric.RemoveAllArrays() + self._PyrQualityMetric.AddObserver( "ModifiedEvent", createModifiedCallback( self ) ) # type: ignore[arg-type] + for measure in getPyramidQualityMeasure().difference( commonCellMetrics ): + self._PyrQualityMetric.AddArray( getQualityMeasureNameFromIndex( measure ) ) + + self._WedgeQualityMetric.RemoveAllArrays() + self._WedgeQualityMetric.AddObserver( + "ModifiedEvent", # type: ignore[arg-type] + createModifiedCallback( self ) ) + for measure in getWedgeQualityMeasure().difference( commonCellMetrics ): + self._WedgeQualityMetric.AddArray( getQualityMeasureNameFromIndex( measure ) ) + + self._HexQualityMetric.RemoveAllArrays() + self._HexQualityMetric.AddObserver( "ModifiedEvent", createModifiedCallback( self ) ) # type: ignore[arg-type] + for measure in getHexQualityMeasure().difference( commonCellMetrics ): + self._HexQualityMetric.AddArray( getQualityMeasureNameFromIndex( measure ) ) + + otherMetrics: set[ int ] = getQualityMetricsOther() + for measure in otherMetrics: + self._commonMeshQualityMetric.AddArray( getQualityMeasureNameFromIndex( measure ) ) + + @smproperty.dataarrayselection( name="CommonCellQualityMetric" ) + def a01SetCommonMetrics( self: Self ) -> vtkDataArraySelection: + """Set polygon quality metrics selection.""" + return self._commonCellQualityMetric + + @smproperty.dataarrayselection( name="TetrahedronSpecificQualityMetric" ) + def a02SetTetMetrics( self: Self ) -> vtkDataArraySelection: + """Set tetra quality metrics selection.""" + return self._tetQualityMetric + + @smproperty.dataarrayselection( name="PyramidSpecificQualityMetric" ) + def a03sSetPyrMetrics( self: Self ) -> vtkDataArraySelection: + """Set Pyramid quality metrics selection.""" + return self._PyrQualityMetric + + @smproperty.dataarrayselection( name="WedgeSpecificQualityMetric" ) + def a03sSetWedgeMetrics( self: Self ) -> vtkDataArraySelection: + """Set Wedge quality metrics selection.""" + return self._WedgeQualityMetric + + @smproperty.dataarrayselection( name="HexahedronSpecificQualityMetric" ) + def a03sSetHexMetrics( self: Self ) -> vtkDataArraySelection: + """Set Hexahdron quality metrics selection.""" + return self._HexQualityMetric + + @smproperty.dataarrayselection( name="OtherMeshQualityMetric" ) + def a04SetOtherMeshMetrics( self: Self ) -> vtkDataArraySelection: + """Set other mesh quality metrics selection.""" + return self._commonMeshQualityMetric + + @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 b01SetSaveToFile( 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 b02SetFileName( 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 b03GroupAdvancedOutputParameters( self: Self ) -> None: + """Organize groups.""" + self.Modified() + + def Modified( self: Self ) -> None: + """Overload Modified method to reset _blockIndex.""" + self._blockIndex = 0 + super().Modified() + + def RequestDataObject( + self: Self, + request: vtkInformation, + inInfoVec: list[ vtkInformationVector ], + outInfoVec: vtkInformationVector, + ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestDataObject. + + Args: + request (vtkInformation): Request + inInfoVec (list[vtkInformationVector]): Input objects + outInfoVec (vtkInformationVector): Output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + inData = self.GetInputData( inInfoVec, 0, 0 ) + outData = self.GetOutputData( outInfoVec, 0 ) + assert inData is not None + if outData is None or ( not outData.IsA( inData.GetClassName() ) ): + outData = inData.NewInstance() + outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) + return super().RequestDataObject( request, inInfoVec, outInfoVec ) + + 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: vtkUnstructuredGrid = self.GetInputData( inInfoVec, 0, 0 ) + outputMesh: vtkUnstructuredGrid = vtkUnstructuredGrid.GetData( outInfoVec, 0 ) + assert inputMesh is not None, "Input server mesh is null." + assert outputMesh is not None, "Output pipeline is null." + + tetraMetrics: set[ int ] = self._getQualityMetricsToUse( self._commonCellQualityMetric ).union( + self._getQualityMetricsToUse( self._tetQualityMetric ) ) + pyrMetrics: set[ int ] = self._getQualityMetricsToUse( self._commonCellQualityMetric ).union( + self._getQualityMetricsToUse( self._PyrQualityMetric ) ) + wedgeMetrics: set[ int ] = self._getQualityMetricsToUse( self._commonCellQualityMetric ).union( + self._getQualityMetricsToUse( self._PyrQualityMetric ) ) + hexaMetrics: set[ int ] = self._getQualityMetricsToUse( self._commonCellQualityMetric ).union( + self._getQualityMetricsToUse( self._PyrQualityMetric ) ) + otherMetrics: set[ int ] = self._getQualityMetricsToUse( self._commonMeshQualityMetric ) + filter: MeshQualityEnhanced = MeshQualityEnhanced() + filter.SetInputDataObject( inputMesh ) + filter.SetCellQualityMetrics( tetraMetrics=tetraMetrics, + pyramidMetrics=pyrMetrics, + wedgeMetrics=wedgeMetrics, + hexaMetrics=hexaMetrics ) + filter.SetOtherMeshQualityMetrics( otherMetrics ) + filter.Update() + + outputMesh.ShallowCopy( filter.GetOutputDataObject( 0 ) ) + + # save to file if asked + if self._saveToFile: + stats: QualityMetricSummary = filter.GetQualityMetricSummary() + self._saveFile( stats ) + self._blockIndex += 1 + return 1 + + def _getQualityMetricsToUse( self: Self, selection: vtkDataArraySelection ) -> set[ int ]: + """Get mesh quality metric indexes from user selection. + + Returns: + list[int]: List of quality metric indexes + """ + metricsNames: set[ str ] = getArrayChoices( selection ) + return { getQualityMeasureIndexFromName( name ) for name in metricsNames } + + def _saveFile( self: Self, stats: QualityMetricSummary ) -> None: + """Export mesh quality metric summary file.""" + try: + if self._filename is None: + print( "Mesh quality summary report file path is undefined." ) + return + # add index for multiblock meshes + index: int = self._filename.rfind( '.' ) + filename: str = self._filename[ :index ] + f"_{self._blockIndex}" + self._filename[ index: ] + fig = stats.plotSummaryFigure() + fig.savefig( filename, dpi=150 ) + print( f"File {filename} was successfully written." ) + except Exception as e: + print( "Error while exporting the file due to:" ) + print( str( e ) ) \ No newline at end of file diff --git a/geos-utils/src/geos/utils/geometryFunctions.py b/geos-utils/src/geos/utils/geometryFunctions.py index 467bac40..138b310f 100644 --- a/geos-utils/src/geos/utils/geometryFunctions.py +++ b/geos-utils/src/geos/utils/geometryFunctions.py @@ -26,11 +26,11 @@ def getChangeOfBasisMatrix( C is then Vc = P.Vb Args: - basisFrom (tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]): origin basis - basisTo (tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]): destination basis + basisFrom (tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]): Origin basis + basisTo (tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]): Destination basis Returns: - npt.NDArray[np.floating[Any]]: change of basis matrix. + npt.NDArray[np.floating[Any]]: Change of basis matrix. """ assert ( basisFrom[ 0 ].size == basisFrom[ 1 ].size ) and ( basisFrom[ 0 ].size == basisFrom[ 2 ].size ), ( "Origin space vectors must have the same size." ) @@ -52,11 +52,11 @@ def computeCoordinatesInNewBasis( vec: npt.NDArray[ np.floating[ Any ] ], """Computes the coordinates of a matrix from a basis B in the new basis B'. Args: - vec (npt.NDArray[np.floating[Any]]): vector to compute the new coordinates + vec (npt.NDArray[np.floating[Any]]): Vector to compute the new coordinates changeOfBasisMatrix (npt.NDArray[np.floating[Any]]): Change of basis matrix Returns: - npt.NDArray[np.floating[Any]]: the new coordinates of the matrix in the basis + npt.NDArray[np.floating[Any]]: The new coordinates of the matrix in the basis B'. """ assert ( vec.size == changeOfBasisMatrix.shape[ 1 ] ), """The size of the vector @@ -80,7 +80,7 @@ def computePlaneFrom3Points( pt3 (npt.NDArray[np.floating[Any]]): 3rd point of the plane. Returns: - tuple[float, float, float, float]: tuple of the 4 coefficients. + tuple[float, float, float, float]: Tuple of the 4 coefficients. """ # plane vectors v1: npt.NDArray[ np.floating[ Any ] ] = pt2 - pt1 @@ -145,3 +145,112 @@ def getPointSideAgainstPlane( dot: float = np.dot( planeNormal, vec ) assert abs( dot ) > EPSILON, "The point is on the plane." return dot > 0 + + +def computeAngleFromPoints( pt1: npt.NDArray[ np.float64 ], pt2: npt.NDArray[ np.float64 ], + pt3: npt.NDArray[ np.float64 ] ) -> float: + """Compute angle from 3 points. + + Args: + pt1 (npt.NDArray[np.float64]): First point + pt2 (npt.NDArray[np.float64]): Second point + pt3 (npt.NDArray[np.float64]): Third point + + Returns: + float: Angle + """ + # compute vectors + vec1: npt.NDArray[ np.float64 ] = pt1 - pt2 + vec2: npt.NDArray[ np.float64 ] = pt3 - pt2 + return computeAngleFromVectors( vec1, vec2 ) + + +def computeAngleFromVectors( + vec1: npt.NDArray[ np.float64 ], + vec2: npt.NDArray[ np.float64 ], +) -> float: + """Compute angle from 2 vectors. + + Args: + vec1 (npt.NDArray[np.float64]): First vector + vec2 (npt.NDArray[np.float64]): Second vector + + Returns: + float: angle + """ + assert abs( np.linalg.norm( vec1 ) ) > 0., "First vector cannot be null" + assert abs( np.linalg.norm( vec2 ) ) > 0., "Second vector cannot be null" + # normalization + vec1_norm: npt.NDArray[ np.float64 ] = vec1 / np.linalg.norm( vec1 ) + vec2_norm: npt.NDArray[ np.float64 ] = vec2 / np.linalg.norm( vec2 ) + dot: float = np.dot( vec1, vec2 ) + det: float = np.linalg.det( ( vec1_norm, vec2_norm ) ) + # truncation due to numerical approximations + if dot > 1.: + dot = 1. + if dot < -1.: + dot = -1. + teta: float = np.arccos( dot ) + if det < 0: + teta = 2.0 * np.pi - teta + return teta + + +def computeCosineFromVectors( + vec1: npt.NDArray[ np.float64 ], + vec2: npt.NDArray[ np.float64 ], +) -> float: + """Compute cosine from 2 vectors. + + Args: + vec1 (npt.NDArray[np.float64]): First vector + vec2 (npt.NDArray[np.float64]): Second vector + + Returns: + float: Cosine + """ + assert abs( np.linalg.norm( vec1 ) ) > 0., "First vector cannot be null" + assert abs( np.linalg.norm( vec2 ) ) > 0., "Second vector cannot be null" + # normalization + vec1_norm: npt.NDArray[ np.float64 ] = vec1 / np.linalg.norm( vec1 ) + vec2_norm: npt.NDArray[ np.float64 ] = vec2 / np.linalg.norm( vec2 ) + return np.dot( vec1_norm, vec2_norm ) + + +def computeNormalFromPoints( pt1: npt.NDArray[ np.float64 ], pt2: npt.NDArray[ np.float64 ], + pt3: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]: + """Compute the normal of a plane defined from 3 points. + + Args: + pt1 (npt.NDArray[np.float64]): First point + pt2 (npt.NDArray[np.float64]): Second point + pt3 (npt.NDArray[np.float64]): Third point + + Returns: + npt.NDArray[np.float64]: Normal vector coordinates + """ + # compute vectors + vec1: npt.NDArray[ np.float64 ] = pt1 - pt2 + vec2: npt.NDArray[ np.float64 ] = pt3 - pt2 + return computeNormalFromVectors( vec1, vec2 ) + + +def computeNormalFromVectors( + vec1: npt.NDArray[ np.float64 ], + vec2: npt.NDArray[ np.float64 ], +) -> npt.NDArray[ np.float64 ]: + """Compute the normal of a plane defined from 2 vectors. + + Args: + vec1 (npt.NDArray[np.float64]): First vector + vec2 (npt.NDArray[np.float64]): Second vector + + Returns: + npt.NDArray[np.float64]: Normal vector coordinates + """ + assert abs( np.linalg.norm( vec1 ) ) > 0., "First and second points must be different" + assert abs( np.linalg.norm( vec2 ) ) > 0., "Second and third points must be different" + # normalization + vec1_norm = vec1 / np.linalg.norm( vec1 ) + vec2_norm = vec2 / np.linalg.norm( vec2 ) + return np.cross( vec1_norm, vec2_norm ) diff --git a/geos-utils/tests/testsGeometryFunctions.py b/geos-utils/tests/testsGeometryFunctions.py index 674123d9..d21442ca 100644 --- a/geos-utils/tests/testsGeometryFunctions.py +++ b/geos-utils/tests/testsGeometryFunctions.py @@ -2,12 +2,13 @@ # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. # SPDX-FileContributor: Martin Lemay # ruff: noqa: E402 # disable Module level import not at top of file -import unittest - -import geos.utils.geometryFunctions as fcts +import pytest import numpy as np import numpy.typing as npt -from typing_extensions import Self +from typing import Iterator +from dataclasses import dataclass +from itertools import combinations +import geos.utils.geometryFunctions as fcts basisCanon: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ] = ( np.array( [ 1.0, 0.0, 0.0 ] ), np.array( [ 0.0, 1.0, 0.0 ] ), @@ -24,112 +25,239 @@ np.array( [ 0.0, 0.0, 2.0 ] ) ) -class TestsGeometryFunctions( unittest.TestCase ): - - def test_getChangeOfBasisMatrixToCanonic( self: Self ) -> None: - """Test change of basis matrix using canonic basis.""" - obtained: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisCanon ) - # matrix where the columns are the vectors - expected: npt.NDArray[ np.float64 ] = np.transpose( np.array( basisTo0 ) ) - self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) ) - - def test_getChangeOfBasisMatrix( self: Self ) -> None: - """Test change of basis matrix format from basis vectors.""" - obtained: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisTo1 ) - expected: npt.NDArray[ np.float64 ] = np.array( [ [ 2.5, 4.5, -0.5 ], [ 0.5, 0.5, 0.5 ], [ -0.5, -0.5, 0.5 ] ] ) - self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) ) - - def test_computeCoordinatesInNewBasisIdentity( self: Self ) -> None: - """Test calculation of coordinates of a vector in the same basis.""" - vec: npt.NDArray[ np.float64 ] = np.array( [ 2.0, 3.0, 4.0 ] ) - - # get change of basis matrix - changeOfBasisMatrix: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisCanon, basisCanon ) - obtained: npt.NDArray[ np.float64 ] = fcts.computeCoordinatesInNewBasis( vec, changeOfBasisMatrix ) - expected: npt.NDArray[ np.float64 ] = vec - self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) ) - - def test_computeCoordinatesInNewBasis( self: Self ) -> None: - """Test calculation of coordinates of a vector in another basis.""" - vec: npt.NDArray[ np.float64 ] = np.array( [ 2.0, 3.0, 4.0 ] ) - - # get change of basis matrix and the inverse - changeOfBasisMatrix: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisTo1 ) - obtained = fcts.computeCoordinatesInNewBasis( vec, changeOfBasisMatrix ) - expected: npt.NDArray[ np.float64 ] = np.array( [ 16.5, 4.5, -0.5 ] ) - self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) ) - - def test_computePlaneFrom3Points( self: Self ) -> None: - """Test calculation of plane coefficients from 3 points.""" - pt1: npt.NDArray[ np.float64 ] = np.array( [ 1.0, 2.0, 1.0 ] ) - pt2: npt.NDArray[ np.float64 ] = np.array( [ 1.0, 1.0, 2.0 ] ) - pt3: npt.NDArray[ np.float64 ] = np.array( [ 3.0, 2.0, 2.0 ] ) - obtained: tuple[ float, float, float, float ] = fcts.computePlaneFrom3Points( pt1, pt2, pt3 ) - expected: tuple[ float, float, float, float ] = ( -1.0, 2.0, 2.0, -5.0 ) - self.assertSequenceEqual( obtained, expected ) - - def test_getPointSideAgainstPlaneAssertion( self: Self ) -> None: - """Test get point side against a plane.""" - planePt: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 0.0, 0.0 ] ) - - # assertion error - Point on the plane - planeNormal: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 0.0, 1.0 ] ) - self.assertRaises( AssertionError, fcts.getPointSideAgainstPlane, planePt, planePt, planeNormal ) - - def test_getPointSideAgainstPlaneHorizontal( self: Self ) -> None: - """Test get point side against a horizontal plane.""" - # horizontal plane - planePt: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 0.0, 0.0 ] ) - listPtsCoords: list[ npt.NDArray[ np.float64 ] ] = [ - np.array( [ 0.5, 0.5, 0.5 ] ), - np.array( [ 0.5, 0.5, -0.5 ] ), - ] - planeNormal: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 0.0, 1.0 ] ) - obtained: list[ bool ] = [] - for ptCoords in listPtsCoords: - side: bool = fcts.getPointSideAgainstPlane( ptCoords, planePt, planeNormal ) - obtained += [ side ] - expected: tuple[ bool, bool ] = ( True, False ) - self.assertSequenceEqual( obtained, expected ) - - def test_getPointSideAgainstPlaneVertical( self: Self ) -> None: - """Test get point side against a vertical plane.""" - # vertical plane - planePt: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 0.0, 0.0 ] ) - listPtsCoords: list[ npt.NDArray[ np.float64 ] ] = [ - np.array( [ 0.5, 0.5, 0.5 ] ), - np.array( [ -0.5, 0.5, 0.5 ] ), - ] - planeNormal: npt.NDArray[ np.float64 ] = np.array( [ 1.0, 0.0, 0.0 ] ) - obtained: list[ bool ] = [] - for ptCoords in listPtsCoords: - side: bool = fcts.getPointSideAgainstPlane( ptCoords, planePt, planeNormal ) - obtained += [ side ] - expected: tuple[ bool, bool ] = ( True, False ) - self.assertSequenceEqual( obtained, expected ) - - def test_getCellSideAgainstPlaneRandom( self: Self ) -> None: - """Test get cell side against a plane.""" - # random plane - planePt = np.array( [ 125.58337, 1386.0465, -2782.502 ] ) - listCellPtsCoords = [ - np.array( [ - [ 135.49551, 1374.7644, -2786.884 ], - [ 125.58337, 1376.7441, -2782.502 ], - [ 132.19525, 1382.2516, -2771.0508 ], - [ 125.58337, 1386.0465, -2782.502 ], - ] ), - np.array( [ - [ 111.9148, 1377.0265, -2764.875 ], - [ 132.19525, 1382.2516, -2771.0508 ], - [ 125.58337, 1376.7441, -2782.502 ], - [ 125.58337, 1386.0465, -2782.502 ], - ] ), - ] - planeNormal = np.array( [ 0.8660075, 0.0, -0.5000311 ] ) - obtained = [] - for cellPtsCoords in listCellPtsCoords: - side: bool = fcts.getCellSideAgainstPlane( cellPtsCoords, planePt, planeNormal ) - obtained += [ side ] - expected = ( True, False ) - self.assertSequenceEqual( obtained, expected ) +def test_getChangeOfBasisMatrixToCanonic() -> None: + """Test change of basis matrix using canonic basis.""" + obtained: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisCanon ) + # matrix where the columns are the vectors + expected: npt.NDArray[ np.float64 ] = np.transpose( np.array( basisTo0 ) ) + assert np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), + equal_nan=True ), f"Expected array is {np.round( expected, 2 ).tolist()}" + + +def test_getChangeOfBasisMatrix() -> None: + """Test change of basis matrix format from basis vectors.""" + obtained: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisTo1 ) + expected: npt.NDArray[ np.float64 ] = np.array( [ [ 2.5, 4.5, -0.5 ], [ 0.5, 0.5, 0.5 ], [ -0.5, -0.5, 0.5 ] ] ) + assert np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), + equal_nan=True ), f"Expected array is {np.round( expected, 2 ).tolist()}" + + +def test_computeCoordinatesInNewBasisIdentity() -> None: + """Test calculation of coordinates of a vector in the same basis.""" + vec: npt.NDArray[ np.float64 ] = np.array( [ 2.0, 3.0, 4.0 ] ) + + # get change of basis matrix + changeOfBasisMatrix: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisCanon, basisCanon ) + obtained: npt.NDArray[ np.float64 ] = fcts.computeCoordinatesInNewBasis( vec, changeOfBasisMatrix ) + expected: npt.NDArray[ np.float64 ] = vec + assert np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), + equal_nan=True ), f"Expected array is {np.round( expected, 2 ).tolist()}" + + +def test_computeCoordinatesInNewBasis() -> None: + """Test calculation of coordinates of a vector in another basis.""" + vec: npt.NDArray[ np.float64 ] = np.array( [ 2.0, 3.0, 4.0 ] ) + + # get change of basis matrix and the inverse + changeOfBasisMatrix: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisTo1 ) + obtained = fcts.computeCoordinatesInNewBasis( vec, changeOfBasisMatrix ) + expected: npt.NDArray[ np.float64 ] = np.array( [ 16.5, 4.5, -0.5 ] ) + assert np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), + equal_nan=True ), f"Expected array is {np.round( expected, 2 ).tolist()}" + + +def test_computePlaneFrom3Points() -> None: + """Test calculation of plane coefficients from 3 points.""" + pt1: npt.NDArray[ np.float64 ] = np.array( [ 1.0, 2.0, 1.0 ] ) + pt2: npt.NDArray[ np.float64 ] = np.array( [ 1.0, 1.0, 2.0 ] ) + pt3: npt.NDArray[ np.float64 ] = np.array( [ 3.0, 2.0, 2.0 ] ) + obtained: tuple[ float, float, float, float ] = fcts.computePlaneFrom3Points( pt1, pt2, pt3 ) + expected: tuple[ float, float, float, float ] = ( -1.0, 2.0, 2.0, -5.0 ) + assert obtained == expected, f"Expected tuple is {expected}" + + +def test_getPointSideAgainstPlaneAssertion() -> None: + """Test get point side against a plane.""" + planePt: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 0.0, 0.0 ] ) + + # assertion error - Point on the plane + planeNormal: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 0.0, 1.0 ] ) + with pytest.raises( AssertionError ): + fcts.getPointSideAgainstPlane( planePt, planePt, planeNormal ) + + +listPtsCoords_all = ( + [ + np.array( [ 0.5, 0.5, 0.5 ] ), + np.array( [ 0.5, 0.5, -0.5 ] ), + ], + [ + np.array( [ 0.5, 0.5, 0.5 ] ), + np.array( [ -0.5, 0.5, 0.5 ] ), + ], +) +planePt_all = ( + np.array( [ 0.0, 0.0, 0.0 ] ), + np.array( [ 0.0, 0.0, 0.0 ] ), +) +planeNormal_all = ( + np.array( [ 0.0, 0.0, 1.0 ] ), + np.array( [ 1.0, 0.0, 0.0 ] ), +) +expected_all = ( + ( True, False ), + ( True, False ), +) + + +@dataclass( frozen=True ) +class TestCasePointSideAgainstPlane: + """Test case.""" + __test__ = False + #: list of points + listPtsCoords: list[ npt.NDArray[ np.float64 ] ] + #: plane point + planePt: npt.NDArray[ np.float64 ] + #: plane normal + planeNormal: npt.NDArray[ np.float64 ] + #: expected result + expected: list[ bool ] + + +def __generate_PointSideAgainstPlane_test_data() -> Iterator[ TestCasePointSideAgainstPlane ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: Iterator on test cases + """ + for listPtsCoords, planePt, planeNormal, expected in zip( listPtsCoords_all, + planePt_all, + planeNormal_all, + expected_all, + strict=True ): + yield TestCasePointSideAgainstPlane( listPtsCoords, planePt, planeNormal, list( expected ) ) + + +ids: tuple[ str, str ] = ( "Horizontal plane", "Vertical plane" ) + + +@pytest.mark.parametrize( "test_case", __generate_PointSideAgainstPlane_test_data(), ids=ids ) +def test_getPointSideAgainstPlane( test_case: TestCasePointSideAgainstPlane ) -> None: + """Test get point side against a plane.""" + obtained: list[ bool ] = [] + for ptCoords in test_case.listPtsCoords: + side: bool = fcts.getPointSideAgainstPlane( ptCoords, test_case.planePt, test_case.planeNormal ) + obtained += [ side ] + assert obtained == test_case.expected, f"Expected tuple is {test_case.expected}" + + +def test_getCellSideAgainstPlaneRandom() -> None: + """Test get cell side against a plane.""" + # random plane + planePt: npt.NDArray[ np.float64 ] = np.array( [ 125.58337, 1386.0465, -2782.502 ] ) + listCellPtsCoords: list[ npt.NDArray[ np.float64 ] ] = [ + np.array( [ + [ 135.49551, 1374.7644, -2786.884 ], + [ 125.58337, 1376.7441, -2782.502 ], + [ 132.19525, 1382.2516, -2771.0508 ], + [ 125.58337, 1386.0465, -2782.502 ], + ] ), + np.array( [ + [ 111.9148, 1377.0265, -2764.875 ], + [ 132.19525, 1382.2516, -2771.0508 ], + [ 125.58337, 1376.7441, -2782.502 ], + [ 125.58337, 1386.0465, -2782.502 ], + ] ), + ] + planeNormal: npt.NDArray[ np.float64 ] = np.array( [ 0.8660075, 0.0, -0.5000311 ] ) + obtained: list[ bool ] = [] + for cellPtsCoords in listCellPtsCoords: + side: bool = fcts.getCellSideAgainstPlane( cellPtsCoords, planePt, planeNormal ) + obtained += [ side ] + expected: list[ bool ] = [ True, False ] + assert obtained == expected, f"Expected tuple is {expected}" + + +pts_all: tuple[ npt.NDArray[ np.float64 ], ...] = ( + np.array( [ 0.0, 0.0, 0.0 ] ), + np.array( [ 1.0, 0.0, 0.0 ] ), + np.array( [ 0.5, 0.5, 0.0 ] ), + np.array( [ -0.5, 0.5, 0.0 ] ), + np.array( [ -0.5, -0.5, 0.0 ] ), + np.array( [ 0.5, -0.5, -1.0 ] ), +) +angleExp_all: tuple[ float, ...] = ( + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., + 0., +) + + +@dataclass( frozen=True ) +class TestCaseAngle: + """Test case.""" + __test__ = False + pt1: npt.NDArray[ np.float64 ] + pt2: npt.NDArray[ np.float64 ] + pt3: npt.NDArray[ np.float64 ] + angleExp: float + + +def __generate_Angle_test_data() -> Iterator[ TestCaseAngle ]: + """Generate test cases. + + Yields: + Iterator[ TestCase ]: Iterator on test cases + """ + print( len( list( combinations( pts_all, 3 ) ) ) ) + for pts, angle in zip( list( combinations( pts_all, 3 ) ), angleExp_all, strict=True ): + yield TestCaseAngle( pts[ 0 ], pts[ 1 ], pts[ 2 ], angle ) + + +@pytest.mark.skip( "Test to fix" ) +@pytest.mark.parametrize( "test_case", __generate_Angle_test_data() ) +def test_computeAngleFromPoints( test_case: TestCaseAngle ) -> None: + """Test computeAngleFromPoints method.""" + obs: float = fcts.computeAngleFromPoints( test_case.pt1, test_case.pt2, test_case.pt3 ) + assert obs == test_case.angleExp + pass + + +@pytest.mark.skip( "Test to fix" ) +@pytest.mark.parametrize( "test_case", __generate_Angle_test_data() ) +def test_computeAngleFromVectors( test_case: TestCaseAngle ) -> None: + """Test computeAngleFromVectors method.""" + vec1: npt.NDArray[ np.float64 ] = test_case.pt1 - test_case.pt2 + vec2: npt.NDArray[ np.float64 ] = test_case.pt3 - test_case.pt2 + obs: float = fcts.computeAngleFromVectors( vec1, vec2 ) + print( f"{test_case.__str__}: {obs}" ) + assert obs == test_case.angleExp + + +@pytest.mark.skip( "Test to fix" ) +def test_computeNormalFromPoints() -> None: + """Test computeNormalFromPoints method.""" + pass + + +@pytest.mark.skip( "Test to fix" ) +def test_computeNormalFromVectors() -> None: + """Test computeNormalFromVectors method.""" + pass