diff --git a/src/Visualization/Python/CMakeLists.txt b/src/Visualization/Python/CMakeLists.txt index eb60bb2ac035..828a17ffdc9d 100644 --- a/src/Visualization/Python/CMakeLists.txt +++ b/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 diff --git a/src/Visualization/Python/GenerateTetrahedralConnectivity.py b/src/Visualization/Python/GenerateTetrahedralConnectivity.py new file mode 100755 index 000000000000..36e53a718453 --- /dev/null +++ b/src/Visualization/Python/GenerateTetrahedralConnectivity.py @@ -0,0 +1,208 @@ +#!/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( + h5files, + subfile_name: str, + relative_paths: bool = True, + start_time: Optional[float] = None, + stop_time: Optional[float] = None, + stride: int = 1, + coordinates: str = "InertialCoordinates", + delete_tetrahedral_connectivity: 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. + + 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: + h5files: List of H5 volume data files. + 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". + delete_tetrahedral_connectivity: Optional. Delete the dataset and return. + Default: False + """ + h5files = [(h5py.File(filename, "a"), filename) for filename in h5files] + + if not subfile_name: + subfiles = available_subfiles( + (h5file for h5file, _ in h5files), extension=".vol" + ) + raise RequiredChoiceError( + ( + "Specify '--subfile-name' / '-d' to select a" + " subfile containing volume data." + ), + choices=subfiles, + ) + + if not subfile_name.endswith(".vol"): + subfile_name += ".vol" + + for h5file, filename in h5files: + # Open subfile + try: + vol_subfile = h5file[subfile_name] + except KeyError as err: + raise ValueError( + f"Could not open subfile name '{subfile_name}' in '{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 delete_tetrahedral_connectivity: + if "tetrahedral_connectivity" in data_at_id: + del data_at_id["tetrahedral_connectivity"] + continue + + 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( + f"Generating tetrahedral connectivity at {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() + ) + + for h5file in h5files: + h5file[0].close() + + +@click.command( + name="generate-tetrahedral-connectivity", + help=generate_tetrahedral_connectivity.__doc__, +) +@click.argument( + "h5files", + 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 " + "not included." + ), +) +@click.option( + "--coordinates", + default="InertialCoordinates", + show_default=True, + help="The coordinates to use for visualization", +) +@click.option( + "--delete-tetrahedral-connectivity", + is_flag=True, + default=False, + help=( + "Delete all the tetrahedral connectivity between start and " + "stop time with the given stride." + ), +) +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"] + ) diff --git a/support/Python/__main__.py b/support/Python/__main__.py index 647e02376fb2..e67320ff261b 100644 --- a/support/Python/__main__.py +++ b/support/Python/__main__.py @@ -29,6 +29,7 @@ def list_commands(self, ctx): "extend-connectivity", "extract-dat", "extract-input", + "generate-tetrahedral-connectivity", "generate-xdmf", "interpolate-to-mesh", "interpolate-to-points", @@ -84,6 +85,12 @@ def get_command(self, ctx, name): ) return extract_input_source_from_h5_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 diff --git a/tests/Unit/Visualization/Python/CMakeLists.txt b/tests/Unit/Visualization/Python/CMakeLists.txt index c3d3c65c4de1..b5fe498818c9 100644 --- a/tests/Unit/Visualization/Python/CMakeLists.txt +++ b/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 diff --git a/tests/Unit/Visualization/Python/Test_GenerateTetrahedralConnectivity.py b/tests/Unit/Visualization/Python/Test_GenerateTetrahedralConnectivity.py new file mode 100644 index 000000000000..7a5675f5c6fe --- /dev/null +++ b/tests/Unit/Visualization/Python/Test_GenerateTetrahedralConnectivity.py @@ -0,0 +1,159 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +import glob +import logging +import os +import shutil +import sys +import unittest + +import h5py +import numpy as np +from click.testing import CliRunner + +import spectre.Informer as spectre_informer +from spectre.support.Logging import configure_logging +from spectre.Visualization.GenerateTetrahedralConnectivity import ( + generate_tetrahedral_connectivity, + generate_tetrahedral_connectivity_command, +) + + +class TestGenerateTetrahedralConnectivity(unittest.TestCase): + def setUp(self): + self.data_dir = os.path.join( + spectre_informer.unit_test_src_path(), "Visualization/Python" + ) + self.test_dir = os.path.join( + spectre_informer.unit_test_build_path(), + "Visualization/TetrahedralConnectivity", + ) + os.makedirs(self.test_dir, exist_ok=True) + self.maxDiff = None + + def tearDown(self): + shutil.rmtree(self.test_dir) + + def test_generate_tetrahedral_connectivity(self): + original_file = os.path.join(self.data_dir, "VolTestData0.h5") + data_file = os.path.join(self.test_dir, "VolTestData0.h5") + shutil.copy2(original_file, data_file) + + self.assertTrue( + "tetrahedral_connectivity" + in h5py.File(data_file, "r")["element_data.vol"][ + "ObservationId1090594013349131584" + ] + ) + + generate_tetrahedral_connectivity( + h5files=[data_file], + subfile_name="element_data", + start_time=0.0, + stop_time=1.0, + stride=1, + coordinates="InertialCoordinates", + delete_tetrahedral_connectivity=True, + ) + + self.assertFalse( + "tetrahedral_connectivity" + in h5py.File(data_file, "r")["element_data.vol"][ + "ObservationId1090594013349131584" + ] + ) + + generate_tetrahedral_connectivity( + h5files=[data_file], + subfile_name="element_data", + start_time=0.0, + stop_time=1.0, + stride=1, + coordinates="InertialCoordinates", + ) + + self.assertTrue( + "tetrahedral_connectivity" + in h5py.File(data_file, "r")["element_data.vol"][ + "ObservationId1090594013349131584" + ] + ) + + original_data = np.asarray( + h5py.File(original_file, "r")["element_data.vol"][ + "ObservationId1090594013349131584" + ]["tetrahedral_connectivity"] + ) + new_data = np.asarray( + h5py.File(data_file, "r")["element_data.vol"][ + "ObservationId1090594013349131584" + ]["tetrahedral_connectivity"] + ) + + def test_cli(self): + original_file = os.path.join(self.data_dir, "VolTestData0.h5") + data_file = os.path.join(self.test_dir, "VolTestData0.h5") + shutil.copy2(original_file, data_file) + data_files = glob.glob(os.path.join(self.test_dir, "VolTestData0.h5")) + + self.assertTrue( + "tetrahedral_connectivity" + in h5py.File(data_file, "r")["element_data.vol"][ + "ObservationId1090594013349131584" + ] + ) + + runner = CliRunner() + result = runner.invoke( + generate_tetrahedral_connectivity_command, + [ + *data_files, + "-d", + "element_data", + "--delete-tetrahedral-connectivity", + ], + catch_exceptions=False, + ) + self.assertEqual(result.exit_code, 0) + + self.assertFalse( + "tetrahedral_connectivity" + in h5py.File(data_file, "r")["element_data.vol"][ + "ObservationId1090594013349131584" + ] + ) + + result = runner.invoke( + generate_tetrahedral_connectivity_command, + [ + *data_files, + "-d", + "element_data", + ], + catch_exceptions=False, + ) + self.assertEqual(result.exit_code, 0) + + self.assertTrue( + "tetrahedral_connectivity" + in h5py.File(data_file, "r")["element_data.vol"][ + "ObservationId1090594013349131584" + ] + ) + + # List available subfiles + result = runner.invoke( + generate_tetrahedral_connectivity_command, + [ + *data_files, + ], + catch_exceptions=False, + ) + self.assertNotEqual(result.exit_code, 0) + self.assertIn("element_data", result.output) + + +if __name__ == "__main__": + configure_logging(log_level=logging.DEBUG) + unittest.main(verbosity=2)