Added volume, self intersecting and wrong support nodes checks to mesh_doctor (#2158)

Author: TotoGaz

geosx_mesh_doctor/checks/__init__.py

@@ -1,5 +1 @@
-COLLOCATES_NODES = "collocated_nodes"
-GENERATE_FRACTURES = "generate_fractures"
-GENERATE_GLOBAL_IDS = "generate_global_ids"
-
-all_checks = dict()
+# Empty

geosx_mesh_doctor/checks/collocated_nodes.py

@@ -1,7 +1,7 @@
 from collections import defaultdict
 from dataclasses import dataclass
 import logging
-from typing import Tuple
+from typing import List
 
 import numpy
 
@@ -22,7 +22,8 @@ class Options:
 
 @dataclass(frozen=True)
 class Result:
-    nodes_buckets: Tuple[Tuple[int]]
+    nodes_buckets: List[List[int]]  # Each bucket contains the duplicated node indices.
+    wrong_support_elements: List[int]  # Element indices with support node indices appearing more than once.
 
 
 def __check(mesh, options: Options) -> Result:
@@ -34,7 +35,6 @@ def __check(mesh, options: Options) -> Result:
     locator.InitPointInsertion(output, points.GetBounds())
 
     # original ids to/from filtered ids.
-    # original_to_filtered = numpy.ones(points.GetNumberOfPoints(), dtype=int) * -1
     filtered_to_original = numpy.ones(points.GetNumberOfPoints(), dtype=int) * -1
 
     rejected_points = defaultdict(list)
@@ -59,7 +59,16 @@ def __check(mesh, options: Options) -> Result:
     for n, ns in rejected_points.items():
         tmp.append((n, *ns))
 
-    return Result(nodes_buckets=tmp)
+    # Checking that the support node indices appear only once per element.
+    wrong_support_elements = []
+    for c in range(mesh.GetNumberOfCells()):
+        cell = mesh.GetCell(c)
+        num_points_per_cell = cell.GetNumberOfPoints()
+        if len({cell.GetPointId(i) for i in range(num_points_per_cell)}) != num_points_per_cell:
+            wrong_support_elements.append(c)
+
+    return Result(nodes_buckets=tmp,
+                  wrong_support_elements=wrong_support_elements)
 
 
 def check(vtk_input_file: str, options: Options) -> Result:

geosx_mesh_doctor/checks/element_volumes.py

@@ -0,0 +1,55 @@
+from dataclasses import dataclass
+import logging
+from typing import List, Tuple
+import uuid
+
+from vtkmodules.vtkFiltersVerdict import (
+    vtkCellSizeFilter
+)
+from vtk.util.numpy_support import (
+    vtk_to_numpy,
+)
+
+
+from . import vtk_utils
+
+
+@dataclass(frozen=True)
+class Options:
+    min_volume: float
+
+
+@dataclass(frozen=True)
+class Result:
+    element_volumes: List[Tuple[int, float]]
+
+
+def __check(mesh, options: Options) -> Result:
+    f = vtkCellSizeFilter()
+
+    f.ComputeAreaOff()
+    f.ComputeLengthOff()
+    f.ComputeSumOff()
+    f.ComputeVertexCountOff()
+    f.ComputeVolumeOn()
+    volume_array_name = "__MESH_DOCTOR_VOLUME-" + str(uuid.uuid4())  # Making the name unique
+    f.SetVolumeArrayName(volume_array_name)
+
+    f.SetInputData(mesh)
+    f.Update()
+    output = f.GetOutput()
+    cd = output.GetCellData()
+    for i in range(cd.GetNumberOfArrays()):
+        if cd.GetArrayName(i) == volume_array_name:
+            volume = vtk_to_numpy(cd.GetArray(i))
+    assert volume is not None
+    small_volumes: List[Tuple[int, float]] = []
+    for i, v in enumerate(volume):
+        if v < options.min_volume:
+            small_volumes.append((i, v))
+    return Result(element_volumes=small_volumes)
+
+
+def check(vtk_input_file: str, options: Options) -> Result:
+    mesh = vtk_utils.read_mesh(vtk_input_file)
+    return __check(mesh, options)

geosx_mesh_doctor/checks/register.py

@@ -0,0 +1,91 @@
+from dataclasses import dataclass
+import logging
+from typing import List
+
+from vtkmodules.vtkFiltersGeneral import (
+    vtkCellValidator
+)
+from vtkmodules.vtkCommonCore import (
+    vtkOutputWindow,
+    vtkFileOutputWindow
+)
+from vtk.util.numpy_support import (
+    vtk_to_numpy,
+)
+
+from . import vtk_utils
+
+
+@dataclass(frozen=True)
+class Options:
+    tolerance: float
+
+
+@dataclass(frozen=True)
+class Result:
+    wrong_number_of_points_elements: List[int]
+    intersecting_edges_elements: List[int]
+    intersecting_faces_elements: List[int]
+    non_contiguous_edges_elements: List[int]
+    non_convex_elements: List[int]
+    faces_are_oriented_incorrectly_elements: List[int]
+
+
+def __check(mesh, options: Options) -> Result:
+    err_out = vtkFileOutputWindow()
+    err_out.SetFileName("/dev/null")  # vtkCellValidator outputs loads for each cell...
+    vtk_std_err_out = vtkOutputWindow()
+    vtk_std_err_out.SetInstance(err_out)
+
+    valid = 0x0
+    wrong_number_of_points = 0x01
+    intersecting_edges = 0x02
+    intersecting_faces = 0x04
+    non_contiguous_edges = 0x08
+    non_convex = 0x10
+    faces_are_oriented_incorrectly = 0x20
+
+    wrong_number_of_points_elements: List[int] = []
+    intersecting_edges_elements: List[int] = []
+    intersecting_faces_elements: List[int] = []
+    non_contiguous_edges_elements: List[int] = []
+    non_convex_elements: List[int] = []
+    faces_are_oriented_incorrectly_elements: List[int] = []
+
+    f = vtkCellValidator()
+    f.SetTolerance(options.tolerance)
+
+    f.SetInputData(mesh)
+    f.Update()
+    output = f.GetOutput()
+
+    cd = output.GetCellData()
+    for i in range(cd.GetNumberOfArrays()):
+        if cd.GetArrayName(i) == "ValidityState":  # Could not change name using the vtk interface.
+            validity = vtk_to_numpy(cd.GetArray(i))
+    assert validity is not None
+    for i, v in enumerate(validity):
+        if not v & valid:
+            if v & wrong_number_of_points:
+                wrong_number_of_points_elements.append(i)
+            if v & intersecting_edges:
+                intersecting_edges_elements.append(i)
+            if v & intersecting_faces:
+                intersecting_faces_elements.append(i)
+            if v & non_contiguous_edges:
+                non_contiguous_edges_elements.append(i)
+            if v & non_convex:
+                non_convex_elements.append(i)
+            if v & faces_are_oriented_incorrectly:
+                faces_are_oriented_incorrectly_elements.append(i)
+    return Result(wrong_number_of_points_elements=wrong_number_of_points_elements,
+                  intersecting_edges_elements=intersecting_edges_elements,
+                  intersecting_faces_elements=intersecting_faces_elements,
+                  non_contiguous_edges_elements=non_contiguous_edges_elements,
+                  non_convex_elements=non_convex_elements,
+                  faces_are_oriented_incorrectly_elements=faces_are_oriented_incorrectly_elements)
+
+
+def check(vtk_input_file: str, options: Options) -> Result:
+    mesh = vtk_utils.read_mesh(vtk_input_file)
+    return __check(mesh, options)

geosx_mesh_doctor/checks/self_intersecting_elements.py

@@ -1,21 +1,36 @@
 import os.path
 import logging
+import sys
 
 from vtkmodules.vtkCommonDataModel import (
-    vtkUnstructuredGrid, )
+    vtkUnstructuredGrid,
+)
 from vtkmodules.vtkIOLegacy import (
     vtkUnstructuredGridWriter,
     vtkUnstructuredGridReader,
 )
 from vtkmodules.vtkIOXML import (
-    vtkXMLUnstructuredGridReader, )
+    vtkXMLUnstructuredGridReader,
+)
 
 
 def read_mesh(vtk_input_file: str) -> vtkUnstructuredGrid:
-    reader = vtkXMLUnstructuredGridReader()    # TODO Find a generic way to read the vtk mesh.
+    # Testing legacy file format.
+    reader = vtkUnstructuredGridReader()
     reader.SetFileName(vtk_input_file)
-    reader.Update()
-    return reader.GetOutput()
+    if reader.IsFileUnstructuredGrid():
+        logging.debug(f"Reading file \"{vtk_input_file}\" using legacy format reader.")
+        reader.Update()
+        return reader.GetOutput()
+    # Now it's xml time!
+    reader = vtkXMLUnstructuredGridReader()
+    if reader.CanReadFile(vtk_input_file):
+        logging.debug(f"Reading file \"{vtk_input_file}\" using XML format reader.")
+        reader.SetFileName(vtk_input_file)
+        reader.Update()
+        return reader.GetOutput()
+    logging.critical(f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\". Dying...")
+    sys.exit(1)
 
 
 def write_mesh(mesh: vtkUnstructuredGrid, output: str) -> None:

geosx_mesh_doctor/checks/vtk_utils.py

@@ -1,17 +1,16 @@
 import logging
 import sys
 
-from checks import all_checks
-from parsing import register
-from parsing import all_checks_helpers
+from parsing import CheckHelper
 from parsing.cli_parsing import parse, parse_and_set_verbosity
+import register
 
 
 def main():
     logging.basicConfig(format='[%(asctime)s][%(levelname)s] %(message)s')
     parse_and_set_verbosity(sys.argv)
-    register.register()
-    args = parse(sys.argv)
+    all_checks, all_checks_helpers = register.register()
+    args = parse(all_checks_helpers, sys.argv)
     logging.info(f"Checking mesh \"{args.vtk_input_file}\".")
     for check_name, check_options in args.checks.items():
         # If there is no option, this means that the check was not requested by the user
@@ -22,7 +21,7 @@ def main():
         except KeyError as e:
             logging.critical(f"Check {check_name} is not a valid check.")
             sys.exit(1)
-        helper = all_checks_helpers[check_name]
+        helper: CheckHelper = all_checks_helpers[check_name]
         result = check(args.vtk_input_file, check_options)
         helper.display_results(check_options, result)
 

geosx_mesh_doctor/mesh_doctor.py

@@ -2,23 +2,15 @@
 from typing import Callable, Any
 
 
+COLLOCATES_NODES = "collocated_nodes"
+ELEMENT_VOLUMES = "element_volumes"
+GENERATE_FRACTURES = "generate_fractures"
+GENERATE_GLOBAL_IDS = "generate_global_ids"
+SELF_INTERSECTING_ELEMENTS = "self_intersecting_elements"
+
+
 @dataclass(frozen=True)
 class CheckHelper:
     parse_cli_options: Callable[[str], Any]
     display_results: Callable[[Any, Any], None]
     get_help: Callable[[None], str]
-
-
-def build_check_helper(module) -> CheckHelper:
-    """
-    If a module has functions `parse_cli_options`, `display_results` and `get_help`,
-    the `CheckHelper` built from those functions.
-    :param module: Any module
-    :return: The CheckHelper instance.
-    """
-    return CheckHelper(parse_cli_options=module.parse_cli_options,
-                       display_results=module.display_results,
-                       get_help=module.get_help)
-
-
-all_checks_helpers = dict()