Improvements in mesh_doctor (#2274)

- Check and fixes for wrong support nodes ordering.
- Checks and fixes polyhedron with normals not directed inward. (Not exposed through the CLI though).
- Checks if the mesh contains "non-conformal" elements (elements are one in front of the other with non matching nodes).
- Checks that all elements are supported by geos, including that polyhedra can be converted into supported elements.
- Generates a cube in vtk.
- The registration of checks has also been improved. And many other code improvements.

Author: TotoGaz

geosx_mesh_doctor/checks/element_volumes.py

@@ -1,12 +1,19 @@
-from dataclasses import dataclass
 import logging
+from dataclasses import dataclass
 from typing import List, Tuple
 import uuid
 
+from vtkmodules.vtkCommonDataModel import (
+    VTK_HEXAHEDRON,
+    VTK_PYRAMID,
+    VTK_TETRA,
+    VTK_WEDGE,
+)
 from vtkmodules.vtkFiltersVerdict import (
-    vtkCellSizeFilter
+    vtkCellSizeFilter,
+    vtkMeshQuality,
 )
-from vtk.util.numpy_support import (
+from vtkmodules.util.numpy_support import (
     vtk_to_numpy,
 )
 
@@ -25,28 +32,48 @@ class Result:
 
 
 def __check(mesh, options: Options) -> Result:
-    f = vtkCellSizeFilter()
+    cs = vtkCellSizeFilter()
 
-    f.ComputeAreaOff()
-    f.ComputeLengthOff()
-    f.ComputeSumOff()
-    f.ComputeVertexCountOff()
-    f.ComputeVolumeOn()
+    cs.ComputeAreaOff()
+    cs.ComputeLengthOff()
+    cs.ComputeSumOff()
+    cs.ComputeVertexCountOff()
+    cs.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))
+    cs.SetVolumeArrayName(volume_array_name)
+
+    cs.SetInputData(mesh)
+    cs.Update()
+
+    mq = vtkMeshQuality()
+    SUPPORTED_TYPES = [VTK_HEXAHEDRON, VTK_TETRA]
+
+    mq.SetTetQualityMeasureToVolume()
+    mq.SetHexQualityMeasureToVolume()
+    if hasattr(mq, "SetPyramidQualityMeasureToVolume"):  # This feature is quite recent
+        mq.SetPyramidQualityMeasureToVolume()
+        SUPPORTED_TYPES.append(VTK_PYRAMID)
+        mq.SetWedgeQualityMeasureToVolume()
+        SUPPORTED_TYPES.append(VTK_WEDGE)
+    else:
+        logging.debug("Your \"pyvtk\" version does not bring pyramid not wedge support with vtkMeshQuality. Using the fallback solution.")
+
+    mq.SetInputData(mesh)
+    mq.Update()
+
+    volume = cs.GetOutput().GetCellData().GetArray(volume_array_name)
+    quality = mq.GetOutput().GetCellData().GetArray("Quality")  # Name is imposed by vtk.
+
     assert volume is not None
+    assert quality is not None
+    volume = vtk_to_numpy(volume)
+    quality = vtk_to_numpy(quality)
     small_volumes: List[Tuple[int, float]] = []
-    for i, v in enumerate(volume):
-        if v < options.min_volume:
-            small_volumes.append((i, v))
+    for i, pack in enumerate(zip(volume, quality)):
+        v, q = pack
+        vol = q if mesh.GetCellType(i) in SUPPORTED_TYPES else v
+        if vol < options.min_volume:
+            small_volumes.append((i, vol))
     return Result(element_volumes=small_volumes)
 
 

geosx_mesh_doctor/checks/fix_elements_orderings.py

@@ -0,0 +1,57 @@
+from dataclasses import dataclass
+import logging
+from typing import List, Dict, Set
+
+from vtkmodules.vtkCommonCore import (
+    vtkIdList,
+)
+
+from . import vtk_utils
+from .vtk_utils import (
+    to_vtk_id_list,
+    VtkOutput,
+)
+
+
+@dataclass(frozen=True)
+class Options:
+    vtk_output: VtkOutput
+    cell_type_to_ordering: Dict[int, List[int]]
+
+
+@dataclass(frozen=True)
+class Result:
+    output: str
+    unchanged_cell_types: Set[int]
+
+
+def __check(mesh, options: Options):
+    cell_type_to_ordering: Dict[int, List[int]] = options.cell_type_to_ordering
+    unchanged_cell_types = []
+    output_mesh = mesh.NewInstance()  # keeping the same instance type.
+    output_mesh.CopyStructure(mesh)
+    output_mesh.CopyAttributes(mesh)
+
+    # `output_mesh` now contains a full copy of the input mesh.
+    # We'll now modify the support nodes orderings in place if needed.
+    cells = output_mesh.GetCells()
+    for cell_idx in range(output_mesh.GetNumberOfCells()):
+        cell_type = output_mesh.GetCell(cell_idx).GetCellType()
+        new_ordering = cell_type_to_ordering.get(cell_type)
+        if new_ordering:
+            support_point_ids = vtkIdList()
+            cells.GetCellAtId(cell_idx, support_point_ids)
+            new_support_point_ids = []
+            for i, v in enumerate(new_ordering):
+                new_support_point_ids.append(support_point_ids.GetId(new_ordering[i]))
+            cells.ReplaceCellAtId(cell_idx, to_vtk_id_list(new_support_point_ids))
+        else:
+            unchanged_cell_types.insert(cell_type)
+    is_written_error = vtk_utils.write_mesh(output_mesh, options.vtk_output)
+    return Result(output=options.vtk_output.output if not is_written_error else "",
+                  unchanged_cell_types=set(unchanged_cell_types))
+
+
+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/generate_cube.py

@@ -0,0 +1,160 @@
+from dataclasses import dataclass
+import logging
+from typing import Sequence, Iterable
+
+import numpy
+
+from vtkmodules.vtkCommonCore import (
+    vtkPoints,
+)
+from vtkmodules.vtkCommonDataModel import (
+    VTK_HEXAHEDRON,
+    vtkCellArray,
+    vtkHexahedron,
+    vtkRectilinearGrid,
+    vtkUnstructuredGrid,
+)
+from vtkmodules.util.numpy_support import (
+    numpy_to_vtk,
+)
+
+from . import vtk_utils
+from .vtk_utils import (
+    VtkOutput,
+)
+
+from .generate_global_ids import __build_global_ids
+
+
+@dataclass(frozen=True)
+class Result:
+    info: str
+
+
+@dataclass(frozen=True)
+class FieldInfo:
+    name: str
+    dimension: int
+    support: str
+
+
+@dataclass(frozen=True)
+class Options:
+    vtk_output: VtkOutput
+    generate_cells_global_ids: bool
+    generate_points_global_ids: bool
+    xs: Sequence[float]
+    ys: Sequence[float]
+    zs: Sequence[float]
+    nxs: Sequence[int]
+    nys: Sequence[int]
+    nzs: Sequence[int]
+    fields: Iterable[FieldInfo]
+
+
+@dataclass(frozen=True)
+class XYZ:
+    x: numpy.array
+    y: numpy.array
+    z: numpy.array
+
+
+def build_rectilinear_blocks_mesh(xyzs: Iterable[XYZ]) -> vtkUnstructuredGrid:
+    """
+    Builds an unstructured vtk grid from the `xyzs` blocks. Kind of InternalMeshGenerator.
+    :param xyzs: The blocks.
+    :return: The unstructured mesh, even if it's topologically structured.
+    """
+    rgs = []
+    for xyz in xyzs:
+        rg = vtkRectilinearGrid()
+        rg.SetDimensions(len(xyz.x), len(xyz.y), len(xyz.z))
+        rg.SetXCoordinates(numpy_to_vtk(xyz.x))
+        rg.SetYCoordinates(numpy_to_vtk(xyz.y))
+        rg.SetZCoordinates(numpy_to_vtk(xyz.z))
+        rgs.append(rg)
+
+    num_points = sum(map(lambda r: r.GetNumberOfPoints(), rgs))
+    num_cells = sum(map(lambda r: r.GetNumberOfCells(), rgs))
+
+    points = vtkPoints()
+    points.Allocate(num_points)
+    for rg in rgs:
+        for i in range(rg.GetNumberOfPoints()):
+            points.InsertNextPoint(rg.GetPoint(i))
+
+    cell_types = [VTK_HEXAHEDRON] * num_cells
+    cells = vtkCellArray()
+    cells.AllocateExact(num_cells, num_cells * 8)
+
+    m = (0, 1, 3, 2, 4, 5, 7, 6)  # VTK_VOXEL and VTK_HEXAHEDRON do not share the same ordering.
+    offset = 0
+    for rg in rgs:
+        for i in range(rg.GetNumberOfCells()):
+            c = rg.GetCell(i)
+            new_cell = vtkHexahedron()
+            for j in range(8):
+                new_cell.GetPointIds().SetId(j, offset + c.GetPointId(m[j]))
+            cells.InsertNextCell(new_cell)
+        offset += rg.GetNumberOfPoints()
+
+    mesh = vtkUnstructuredGrid()
+    mesh.SetPoints(points)
+    mesh.SetCells(cell_types, cells)
+
+    return mesh
+
+
+def __add_fields(mesh: vtkUnstructuredGrid, fields: Iterable[FieldInfo]) -> vtkUnstructuredGrid:
+    for field_info in fields:
+        if field_info.support == "CELLS":
+            data = mesh.GetCellData()
+            n = mesh.GetNumberOfCells()
+        elif field_info.support == "POINTS":
+            data = mesh.GetPointData()
+            n = mesh.GetNumberOfPoints()
+        array = numpy.ones((n, field_info.dimension), dtype=float)
+        vtk_array = numpy_to_vtk(array)
+        vtk_array.SetName(field_info.name)
+        data.AddArray(vtk_array)
+    return mesh
+
+
+def __build(options: Options):
+    def build_coordinates(positions, num_elements):
+        result = []
+        it = zip(zip(positions, positions[1:]), num_elements)
+        try:
+            coords, n = next(it)
+            while True:
+                start, stop = coords
+                end_point = False
+                tmp = numpy.linspace(start=start, stop=stop, num=n+end_point, endpoint=end_point)
+                coords, n = next(it)
+                result.append(tmp)
+        except StopIteration:
+            end_point = True
+            tmp = numpy.linspace(start=start, stop=stop, num=n+end_point, endpoint=end_point)
+            result.append(tmp)
+        return numpy.concatenate(result)
+    x = build_coordinates(options.xs, options.nxs)
+    y = build_coordinates(options.ys, options.nys)
+    z = build_coordinates(options.zs, options.nzs)
+    cube = build_rectilinear_blocks_mesh((XYZ(x, y, z),))
+    cube = __add_fields(cube, options.fields)
+    __build_global_ids(cube, options.generate_cells_global_ids, options.generate_points_global_ids)
+    return cube
+
+
+def __check(options: Options) -> Result:
+    output_mesh = __build(options)
+    vtk_utils.write_mesh(output_mesh, options.vtk_output)
+    return Result(info=f"Mesh was written to {options.vtk_output.output}")
+
+
+def check(vtk_input_file: str, options: Options) -> Result:
+    try:
+        return __check(options)
+    except BaseException as e:
+        logging.error(e)
+        return Result(info="Something went wrong.")