Merge pull request #6303 from nilsdeppe/tetrahedral_connectivity

Add support for tetrahedral connectivity and use scipy to generate it

src/Visualization/Python/CMakeLists.txt

@@ -7,6 +7,7 @@ spectre_python_add_module(
   Visualization
   PYTHON_FILES
   __init__.py
+  GenerateTetrahedralConnectivity.py
   GenerateXdmf.py
   InterpolateToMesh.py
   OpenVolfiles.py

src/Visualization/Python/GenerateTetrahedralConnectivity.py

@@ -0,0 +1,215 @@
+#!/usr/bin/env python
+
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+import logging
+from typing import Optional
+
+import click
+import h5py
+import numpy as np
+import rich
+import scipy.spatial as spatial
+
+from spectre.support.CliExceptions import RequiredChoiceError
+from spectre.support.Logging import configure_logging
+from spectre.Visualization.GenerateXdmf import available_subfiles
+
+logger = logging.getLogger(__name__)
+
+
+def generate_tetrahedral_connectivity(
+    h5file,
+    subfile_name: str,
+    relative_paths: bool = True,
+    start_time: Optional[float] = None,
+    stop_time: Optional[float] = None,
+    stride: int = 1,
+    coordinates: str = "InertialCoordinates",
+    force: bool = False,
+):
+    """Generate tetrahedral connectivity using scipy.spatial.Delaunay
+
+    Given the coordinates, this generates a tetrahedral connectivity that can
+    be read in by ParaView or other visualization software (e.g. VisIt or yt).
+    It uses the scipy.spatial.Delaunay class to generate the connectivity, which
+    uses the Qhull library (github.com/qhull/qhull/), which uses the quickhull
+    algorithm and scales as O(n^2) in the worst case. Thus, generating the
+    connectivity can take several minutes per temporal ID. Unfortunately,
+    depending on the particular grid point distribution, generating the
+    connectivity may even fail in qhull, which is unfortunately difficult
+    to fix.
+
+    You should combine the volume HDF5 files into one before running this in
+    order to get a fully connected domain.
+
+    Note that while ParaView has a Delaunay3D filter, it is much slower than
+    qhull, needs to be rerun every time ParaView is opened while the qhull
+    output is stored, and the Delaunay3D filter sometimes produces nonsense
+    connectivity that is very difficult to debug and fix.
+
+    After the tetrahedral connectivity has been written you can run
+    'generate-xdmf' with the flag '--use-tetrahedral-connectivity'. ParaView
+    volume renderings are sometimes nicer with tetrahedral connectivity and this
+    can be used to fill gaps between finite-difference or Gauss elements.
+
+    If this algorithm is too slow, one possible improvement is to apply
+    qhull to each block of the domain and then connect the blocks to each
+    other separately. This keeps the number of grid points lower for each
+    invocation of qhull, which likely reduces total runtime and may also
+    reduce or eliminate failure cases. In the ideal case we would apply
+    qhull to each element and then connect elements that are using FD or
+    Gauss points to their neighbors.
+
+    \f
+    Arguments:
+      h5file: The HDF5 file on which to run.
+      subfile_name: Volume data subfile in the H5 files.
+      start_time: Optional. The earliest time at which to start visualizing. The
+        start-time value is included.
+      stop_time: Optional. The time at which to stop visualizing. The stop-time
+        value is not included.
+      stride: Optional. View only every stride'th time step.
+      coordinates: Optional. Name of coordinates dataset. Default:
+        "InertialCoordinates".
+      force: Optional. Overwrite the existing tetrahedral connectivity.
+        Default: False
+    """
+    filename = h5file
+    h5file = h5py.File(filename, "a")
+
+    if not subfile_name:
+        subfiles = available_subfiles(h5file, extension=".vol")
+        if len(subfiles) == 1:
+            subfile_name = subfiles[0]
+            logger.info(
+                f"Selected subfile {subfile_name} (the only available one)."
+            )
+        else:
+            raise RequiredChoiceError(
+                (
+                    "Specify '--subfile-name' / '-d' to select a"
+                    " subfile containing volume data."
+                ),
+                choices=subfiles,
+            )
+
+    if not subfile_name.endswith(".vol"):
+        subfile_name += ".vol"
+
+    # Open subfile
+    try:
+        vol_subfile = h5file[subfile_name]
+    except KeyError as err:
+        raise ValueError(
+            f"Could not open subfile name '{subfile_name}' in"
+            f" '{filename}'. Available subfiles: "
+            + str(available_subfiles(h5file, extension=".vol"))
+        ) from err
+
+    # Sort timesteps by time
+    temporal_ids_and_values = sorted(
+        [
+            (key, vol_subfile[key].attrs["observation_value"])
+            for key in vol_subfile.keys()
+        ],
+        key=lambda key_and_time: key_and_time[1],
+    )
+
+    # Stride through timesteps
+    for temporal_id, time in temporal_ids_and_values[::stride]:
+        # Filter by start and end time
+        if start_time is not None and time < start_time:
+            continue
+        if stop_time is not None and time > stop_time:
+            break
+
+        data_at_id = vol_subfile[temporal_id]
+
+        if force and ("tetrahedral_connectivity" in data_at_id):
+            del data_at_id["tetrahedral_connectivity"]
+
+        x_coords = np.asarray(data_at_id[coordinates + "_x"])
+        y_coords = np.asarray(data_at_id[coordinates + "_y"])
+        if (coordinates + "_z") in data_at_id:
+            z_coords = np.asarray(data_at_id[coordinates + "_z"])
+            coords = np.column_stack((x_coords, y_coords, z_coords))
+        else:
+            coords = np.column_stack((x_coords, y_coords))
+
+        logger.info(
+            "Generating tetrahedral connectivity at"
+            f" {temporal_id}/{time}. This may take a few minutes and"
+            " may even fail depending on the grid structure."
+        )
+        delaunay = spatial.Delaunay(coords)
+        data_at_id.create_dataset(
+            "tetrahedral_connectivity",
+            data=delaunay.simplices.flatten(),
+        )
+
+    h5file.close()
+
+
+@click.command(
+    name="generate-tetrahedral-connectivity",
+    help=generate_tetrahedral_connectivity.__doc__,
+)
+@click.argument(
+    "h5file",
+    type=click.Path(exists=True, file_okay=True, dir_okay=False, readable=True),
+    nargs=1,
+    required=True,
+)
+@click.option(
+    "--subfile-name",
+    "-d",
+    help=(
+        "Name of the volume data subfile in the H5 files. A '.vol' extension is"
+        " added if needed. If unspecified, and the first H5 file contains only"
+        " a single '.vol' subfile, choose that. Otherwise, list all '.vol'"
+        " subfiles and exit."
+    ),
+)
+@click.option(
+    "--stride", default=1, type=int, help="View only every stride'th time step"
+)
+@click.option(
+    "--start-time",
+    type=float,
+    help=(
+        "The earliest time at which to start visualizing. The start-time "
+        "value is included."
+    ),
+)
+@click.option(
+    "--stop-time",
+    type=float,
+    help=(
+        "The time at which to stop visualizing. The stop-time value is "
+        "included."
+    ),
+)
+@click.option(
+    "--coordinates",
+    default="InertialCoordinates",
+    show_default=True,
+    help="The coordinates to use for visualization",
+)
+@click.option(
+    "--force",
+    is_flag=True,
+    default=False,
+    help="Overwrite existing tetrahedral connectivity.",
+)
+def generate_tetrahedral_connectivity_command(**kwargs):
+    _rich_traceback_guard = True  # Hide traceback until here
+    generate_tetrahedral_connectivity(**kwargs)
+
+
+if __name__ == "__main__":
+    configure_logging(log_level=logging.INFO)
+    generate_tetrahedral_connectivity_command(
+        help_option_names=["-h", "--help"]
+    )

src/Visualization/Python/GenerateXdmf.py

@@ -32,6 +32,8 @@ def _list_tensor_components(observation):
     components.remove("connectivity")
     if "pole_connectivity" in components:
         components.remove("pole_connectivity")
+    if "tetrahedral_connectivity" in components:
+        components.remove("tetrahedral_connectivity")
     components.remove("total_extents")
     components.remove("grid_names")
     components.remove("bases")
@@ -57,6 +59,7 @@ def _xmf_topology(
     num_vertices = {
         "Hexahedron": 8,
         "Quadrilateral": 4,
+        "Tetrahedron": 4,
         "Triangle": 3,
     }[topology_type]
     num_cells = len(observation[connectivity_name]) // num_vertices
@@ -162,6 +165,7 @@ def _xmf_grid(
     temporal_id: str,
     coordinates: str,
     filling_poles: bool = False,
+    use_tetrahedral_connectivity: bool = False,
 ) -> ET.Element:
     # Make sure the coordinates are found in the file. We assume there should
     # always be an x-coordinate.
@@ -219,10 +223,17 @@ def _xmf_grid(
             )
     else:
         # Cover volume
+        if use_tetrahedral_connectivity:
+            topology_type = {3: "Tetrahedron", 2: "Triangle"}[topo_dim]
+            connectivity_name = "tetrahedral_connectivity"
+        else:
+            topology_type = {3: "Hexahedron", 2: "Quadrilateral"}[topo_dim]
+            connectivity_name = "connectivity"
+
         xmf_topology = _xmf_topology(
             observation,
-            topology_type={3: "Hexahedron", 2: "Quadrilateral"}[topo_dim],
-            connectivity_name="connectivity",
+            topology_type=topology_type,
+            connectivity_name=connectivity_name,
             grid_path=grid_path,
         )
     xmf_grid.append(xmf_topology)
@@ -279,6 +290,7 @@ def generate_xdmf(
     stop_time: Optional[float] = None,
     stride: int = 1,
     coordinates: str = "InertialCoordinates",
+    use_tetrahedral_connectivity: bool = False,
 ):
     """Generate an XDMF file for ParaView and VisIt
 
@@ -304,6 +316,8 @@ def generate_xdmf(
       stride: Optional. View only every stride'th time step.
       coordinates: Optional. Name of coordinates dataset. Default:
         "InertialCoordinates".
+      use_tetrahedral_connectivity: Optional. Use "tetrahedral_connectivity".
+        Default: False
     """
     h5files = [(h5py.File(filename, "r"), filename) for filename in h5files]
 
@@ -402,6 +416,7 @@ def generate_xdmf(
                     subfile_name=subfile_name,
                     temporal_id=temporal_id,
                     coordinates=coordinates,
+                    use_tetrahedral_connectivity=use_tetrahedral_connectivity,
                 )
             )
             # Connect poles if the data is a 2D surface in 3D
@@ -415,6 +430,9 @@ def generate_xdmf(
                         temporal_id=temporal_id,
                         coordinates=coordinates,
                         filling_poles=True,
+                        use_tetrahedral_connectivity=(
+                            use_tetrahedral_connectivity
+                        ),
                     )
                 )
 
@@ -492,7 +510,7 @@ def generate_xdmf(
     type=float,
     help=(
         "The time at which to stop visualizing. The stop-time value is "
-        "not included."
+        "included."
     ),
 )
 @click.option(
@@ -501,6 +519,17 @@ def generate_xdmf(
     show_default=True,
     help="The coordinates to use for visualization",
 )
+@click.option(
+    "--use-tetrahedral-connectivity",
+    is_flag=True,
+    default=False,
+    help=(
+        "Use a tetrahedral connectivity called tetrahedral_connectivity in "
+        "the HDF5 file. See the generate-tetrahedral-connectivity CLI for "
+        "information on how to generate tetrahedral connectivity and what it "
+        "can be useful for."
+    ),
+)
 def generate_xdmf_command(**kwargs):
     _rich_traceback_guard = True  # Hide traceback until here
     generate_xdmf(**kwargs)

support/Python/__main__.py

@@ -30,6 +30,7 @@ def list_commands(self, ctx):
             "extract-dat",
             "extract-input",
             "find-radial-surface",
+            "generate-tetrahedral-connectivity",
             "generate-xdmf",
             "interpolate-to-mesh",
             "interpolate-to-points",
@@ -91,6 +92,12 @@ def get_command(self, ctx, name):
             )
 
             return find_radial_surface_command
+        elif name == "generate-tetrahedral-connectivity":
+            from spectre.Visualization.GenerateTetrahedralConnectivity import (
+                generate_tetrahedral_connectivity_command,
+            )
+
+            return generate_tetrahedral_connectivity_command
         elif name == "generate-xdmf":
             from spectre.Visualization.GenerateXdmf import generate_xdmf_command
 

tests/Unit/Visualization/Python/CMakeLists.txt

@@ -1,6 +1,12 @@
 # Distributed under the MIT License.
 # See LICENSE.txt for details.
 
+spectre_add_python_bindings_test(
+  "Unit.Visualization.Python.GenerateTetrahedralConnectivity"
+  Test_GenerateTetrahedralConnectivity.py
+  "unit;visualization;python"
+  None)
+
 spectre_add_python_bindings_test(
   "Unit.Visualization.Python.GenerateXdmf"
   Test_GenerateXdmf.py