Add a vtk_mesh_is_valid helper to check validity of VTK meshes

src/ansys/dpf/core/vtk_helper.py

@@ -22,7 +22,9 @@
 
 """Provides for vtk helper functions."""
 
+from dataclasses import dataclass
 from typing import Union
+import warnings
 
 import numpy as np
 import pyvista as pv
@@ -363,6 +365,7 @@
     mesh: dpf.MeshedRegion,
     nodes: Union[dpf.Field, None] = None,
     as_linear: bool = True,
+    check_validity: bool = False,
 ) -> pv.UnstructuredGrid:
     """Return a pyvista UnstructuredGrid given a pydpf MeshedRegion.
 
@@ -377,18 +380,159 @@
     as_linear:
         Export quadratic surface elements as linear.
 
-    export_faces:
-        Whether to export face elements along with volume elements for fluid meshes.
+    check_validity:
+        Whether to run the VTK cell validity check on the generated mesh and warn if not valid.
 
     Returns
     -------
     grid:
         UnstructuredGrid corresponding to the DPF mesh.
     """
     try:
-        return dpf_mesh_to_vtk_op(mesh, nodes, as_linear)
+        grid = dpf_mesh_to_vtk_op(mesh, nodes, as_linear)
     except (AttributeError, KeyError, errors.DPFServerException):
-        return dpf_mesh_to_vtk_py(mesh, nodes, as_linear)
+        grid = dpf_mesh_to_vtk_py(mesh, nodes, as_linear)
+    if check_validity:
+        validity = vtk_mesh_is_valid(grid)
+        if not validity.valid:
+            warnings.warn(f"\nVTK mesh validity check\n{validity.msg}")
+    return grid
+
+
+@dataclass
+class VTKMeshValidity:
+    """Dataclass containing the results of a call to vtk_mesh_is_valid.
+
+    valid:
+        Whether the vtk mesh is valid according to the vtkCellValidator.
+    message:
+        Output message.
+    validity_grid:
+        A copy of the original grid, with validity fields.
+    wrong_number_of_points:
+        List of indexes of elements with the wrong number of points.
+    intersecting_edges:
+        List of indexes of elements with intersecting edges.
+    intersecting_faces:
+        List of indexes of elements with intersecting faces.
+    non_contiguous_edges:
+        List of indexes of elements with non-contiguous edges.
+    non_convex:
+        List of indexes of elements with non-convex shape.
+    inverted_faces:
+        List of indexes of elements with inverted faces.
+    """
+
+    valid: bool
+    msg: str
+    grid: pv.UnstructuredGrid
+    wrong_number_of_points: np.ndarray
+    intersecting_edges: np.ndarray
+    intersecting_faces: np.ndarray
+    non_contiguous_edges: np.ndarray
+    non_convex: np.ndarray
+    inverted_faces: np.ndarray
+
+
+def vtk_mesh_is_valid(grid: pv.UnstructuredGrid, verbose: bool = False) -> VTKMeshValidity:
+    """Run a vtk.CellValidator filter on the input grid.
+
+    Parameters
+    ----------
+    grid:
+        A vtk mesh to validate.
+    verbose:
+        Whether to print the complete validation.
+
+    Returns
+    -------
+    validity:
+        A dataclass containing the results of the validator.
+    """
+    from enum import Enum
+
+    from vtkmodules.util.numpy_support import vtk_to_numpy
+    from vtkmodules.vtkFiltersGeneral import vtkCellValidator
+
+    # Prepare the Enum of possible validity states
+    class State(Enum):
+        Valid = 0
+        WrongNumberOfPoints = (1,)
+        IntersectingEdges = (2,)
+        IntersectingFaces = (4,)
+        NoncontiguousEdges = (8,)
+        Nonconvex = (16,)
+        FacesAreOrientedIncorrectly = (32,)
+
+    # Run the cell validator
+    cell_validator = vtkCellValidator()
+    cell_validator.SetInputData(grid)
+    cell_validator.Update()
+    # Get the states for all cells as a numpy array
+    validity_grid = cell_validator.GetUnstructuredGridOutput()
+    cell_states = vtk_to_numpy(validity_grid.GetCellData().GetArray("ValidityState"))
+    # Check for invalid states
+    elem_with_wrong_number_of_nodes = np.where(cell_states & State.WrongNumberOfPoints.value)[0]
+    elem_with_intersecting_edges = np.where(cell_states & State.IntersectingEdges.value)[0]
+    elem_with_intersecting_faces = np.where(cell_states & State.IntersectingFaces.value)[0]
+    elem_with_non_contiguous_edges = np.where(cell_states & State.NoncontiguousEdges.value)[0]
+    elem_with_non_convex_shape = np.where(cell_states & State.Nonconvex.value)[0]
+    elem_with_badly_oriented_faces = np.where(
+        cell_states & State.FacesAreOrientedIncorrectly.value
+    )[0]
+
+    # Build list of number of elements failing each test
+    failing_elements_number = [
+        len(elem_with_wrong_number_of_nodes),
+        len(elem_with_intersecting_edges),
+        len(elem_with_intersecting_faces),
+        len(elem_with_non_contiguous_edges),
+        len(elem_with_non_convex_shape),
+        len(elem_with_badly_oriented_faces),
+    ]
+    # Define whether mesh is valid
+    mesh_is_valid = np.sum(failing_elements_number) == 0
+    # Build output message
+    out_msg = ""
+    if mesh_is_valid:
+        out_msg += "Mesh is valid."
+    else:
+        out_msg += "Mesh is invalid because of (by index):\n"
+        if failing_elements_number[0] > 0:
+            out_msg += (
+                f"  - {failing_elements_number[0]} elements with the wrong number of points:\n"
+            )
+            out_msg += f"      {elem_with_wrong_number_of_nodes}\n"
+        if failing_elements_number[1] > 0:
+            out_msg += f"  - {failing_elements_number[1]} elements with intersecting edges:\n"
+            out_msg += f"      {elem_with_intersecting_edges}\n"
+        if failing_elements_number[2] > 0:
+            out_msg += f"  - {failing_elements_number[2]} elements with intersecting faces:\n"
+            out_msg += f"      {elem_with_intersecting_faces}\n"
+        if failing_elements_number[3] > 0:
+            out_msg += f"  - {failing_elements_number[3]} elements with non contiguous edges:\n"
+            out_msg += f"      {elem_with_non_contiguous_edges}\n"
+        if failing_elements_number[4] > 0:
+            out_msg += f"  - {failing_elements_number[4]} elements with non convex shape:\n"
+            out_msg += f"      {elem_with_non_convex_shape}\n"
+        if failing_elements_number[5] > 0:
+            out_msg += f"  - {failing_elements_number[5]} elements with bad face orientations:\n"
+            out_msg += f"      {elem_with_badly_oriented_faces}\n"
+    if verbose:
+        print(out_msg)
+    validity_grid = pv.UnstructuredGrid(validity_grid)
+    validity_grid.set_active_scalars("ValidityState")
+    return VTKMeshValidity(
+        valid=mesh_is_valid,
+        msg=out_msg,
+        grid=validity_grid,
+        wrong_number_of_points=elem_with_wrong_number_of_nodes,
+        intersecting_edges=elem_with_intersecting_edges,
+        intersecting_faces=elem_with_intersecting_faces,
+        non_contiguous_edges=elem_with_non_contiguous_edges,
+        non_convex=elem_with_non_convex_shape,
+        inverted_faces=elem_with_badly_oriented_faces,
+    )
 
 
 def vtk_update_coordinates(vtk_grid, coordinates_array):

tests/test_vtk_translate.py

@@ -32,6 +32,7 @@
     dpf_mesh_to_vtk,
     dpf_meshes_to_vtk,
     dpf_property_field_to_vtk,
+    vtk_mesh_is_valid,
 )
 import conftest
 
@@ -230,3 +231,90 @@ def test_append_fields_container_to_grid(simple_rst, server_type):
     assert isinstance(ug, pv.UnstructuredGrid)
     assert "disp {'time': 1}" in ug.point_data.keys()
     assert "volume {'time': 1}" in ug.cell_data.keys()
+
+
+@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista")
+def test_vtk_mesh_is_valid_polyhedron():
+    # Element type is polyhedron
+    cell_types = [pv.CellType.POLYHEDRON]
+
+    # Start with a valid element
+    nodes_1 = [
+        [0.0, 0.0, 0.0],
+        [1.0, 0.0, 0.0],
+        [0.0, 1.0, 0.0],
+        [0.0, 0.0, 0.5],
+        [1.0, 0.0, 0.5],
+        [0.0, 1.0, 0.5],
+    ]
+    cells_1 = [5, 4, 4, 1, 2, 5, 4, 3, 0, 1, 4, 3, 2, 1, 0, 3, 3, 4, 5, 4, 5, 2, 0, 3]
+    grid = pv.UnstructuredGrid([len(cells_1), *cells_1], cell_types, nodes_1)
+    validity = vtk_mesh_is_valid(grid)
+    print(validity)
+    assert validity.valid
+    assert "valid" in validity.msg
+    assert validity.grid.active_scalars_name == "ValidityState"
+    assert len(validity.wrong_number_of_points) == 0
+    assert len(validity.intersecting_edges) == 0
+    assert len(validity.intersecting_faces) == 0
+    assert len(validity.non_contiguous_edges) == 0
+    assert len(validity.non_convex) == 0
+    assert len(validity.inverted_faces) == 0
+
+    # Move one node
+    nodes_2 = [
+        [0.0, 0.0, 0.0],
+        [1.0, 0.0, 0.0],
+        [0.0, 1.0, 0.0],
+        [0.0, -0.05, 0.5],  # Moved one node along Y axis
+        [1.0, 0.0, 0.5],
+        [0.0, 1.0, 0.5],
+    ]
+    grid = pv.UnstructuredGrid([len(cells_1), *cells_1], cell_types, nodes_2)
+    validity = vtk_mesh_is_valid(grid)
+    print(validity)
+    assert not validity.valid  # For some reason this element is found to be non-convex
+    assert len(validity.wrong_number_of_points) == 0
+    assert len(validity.intersecting_edges) == 0
+    assert len(validity.intersecting_faces) == 0
+    assert len(validity.non_contiguous_edges) == 0
+    assert len(validity.non_convex) == 1
+    assert len(validity.inverted_faces) == 0
+
+    # Invert one face
+    cells_2 = [
+        5,
+        4,
+        4,
+        1,
+        2,
+        5,
+        4,
+        3,
+        0,
+        1,
+        4,
+        3,
+        2,
+        1,
+        0,
+        3,
+        5,
+        4,
+        3,  # Inverted face
+        4,
+        5,
+        2,
+        0,
+        3,
+    ]
+    grid = pv.UnstructuredGrid([len(cells_2), *cells_2], cell_types, nodes_1)
+    validity = vtk_mesh_is_valid(grid)
+    print(validity)
+    assert not validity.valid  # Non-convex AND bad face orientation
+    assert len(validity.wrong_number_of_points) == 0
+    assert len(validity.intersecting_edges) == 0
+    assert len(validity.intersecting_faces) == 0
+    assert len(validity.non_contiguous_edges) == 0
+    assert len(validity.non_convex) == 1
+    assert len(validity.inverted_faces) == 1