Skip to content

Commit

Permalink
fixup! Add GenerateTetrahedralConnectivity
Browse files Browse the repository at this point in the history
  • Loading branch information
@nilsdeppe
nilsdeppe committed Sep 28, 2024
1 parent 13782d3 commit fc9da37
Show file tree
Hide file tree
Showing 2 changed files with 99 additions and 102 deletions.
159 changes: 83 additions & 76 deletions src/Visualization/Python/GenerateTetrahedralConnectivity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,14 +20,14 @@


def generate_tetrahedral_connectivity(
h5files,
h5file,
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,
force: bool = False,
):
"""Generate tetrahedral connectivity using scipy.spatial.Delaunay
Expand All @@ -41,6 +41,14 @@ def generate_tetrahedral_connectivity(
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
Expand All @@ -56,7 +64,7 @@ def generate_tetrahedral_connectivity(
\f
Arguments:
h5files: List of H5 volume data files.
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.
Expand All @@ -65,91 +73,93 @@ def generate_tetrahedral_connectivity(
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.
force: Optional. Overwrite the existing tetrahedral connectivity.
Default: False
"""
h5files = [(h5py.File(filename, "a"), filename) for filename in h5files]
filename = h5file
h5file = h5py.File(filename, "a")

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,
)
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"

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))
# 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],
)

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()
)
# 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(),
)

for h5file in h5files:
h5file[0].close()
h5file.close()


@click.command(
name="generate-tetrahedral-connectivity",
help=generate_tetrahedral_connectivity.__doc__,
)
@click.argument(
"h5files",
"h5file",
type=click.Path(exists=True, file_okay=True, dir_okay=False, readable=True),
nargs=-1,
nargs=1,
required=True,
)
@click.option(
Expand Down Expand Up @@ -178,7 +188,7 @@ def generate_tetrahedral_connectivity(
type=float,
help=(
"The time at which to stop visualizing. The stop-time value is "
"not included."
"included."
),
)
@click.option(
Expand All @@ -188,13 +198,10 @@ def generate_tetrahedral_connectivity(
help="The coordinates to use for visualization",
)
@click.option(
"--delete-tetrahedral-connectivity",
"--force",
is_flag=True,
default=False,
help=(
"Delete all the tetrahedral connectivity between start and "
"stop time with the given stride."
),
help="Overwrite existing tetrahedral connectivity.",
)
def generate_tetrahedral_connectivity_command(**kwargs):
_rich_traceback_guard = True # Hide traceback until here
Expand Down
42 changes: 16 additions & 26 deletions tests/Unit/Visualization/Python/Test_GenerateTetrahedralConnectivity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,30 +48,31 @@ def test_generate_tetrahedral_connectivity(self):
)

generate_tetrahedral_connectivity(
h5files=[data_file],
h5file=data_file,
subfile_name="element_data",
start_time=0.0,
stop_time=1.0,
stride=1,
coordinates="InertialCoordinates",
delete_tetrahedral_connectivity=True,
force=True,
)

self.assertFalse(
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",
)
with self.assertRaisesRegex(ValueError, "Unable to create dataset"):
generate_tetrahedral_connectivity(
h5file=data_file,
subfile_name="element_data",
start_time=0.0,
stop_time=1.0,
stride=1,
coordinates="InertialCoordinates",
)

self.assertTrue(
"tetrahedral_connectivity"
Expand Down Expand Up @@ -111,13 +112,13 @@ def test_cli(self):
*data_files,
"-d",
"element_data",
"--delete-tetrahedral-connectivity",
"--force",
],
catch_exceptions=False,
)
self.assertEqual(result.exit_code, 0)

self.assertFalse(
self.assertTrue(
"tetrahedral_connectivity"
in h5py.File(data_file, "r")["element_data.vol"][
"ObservationId1090594013349131584"
Expand All @@ -131,9 +132,9 @@ def test_cli(self):
"-d",
"element_data",
],
catch_exceptions=False,
catch_exceptions=True,
)
self.assertEqual(result.exit_code, 0)
self.assertEqual(result.exit_code, 1)

self.assertTrue(
"tetrahedral_connectivity"
Expand All @@ -142,17 +143,6 @@ def test_cli(self):
]
)

# 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)
Expand Down

0 comments on commit fc9da37

Please sign in to comment.