From 749928d76bcec54994d19ca894c82e2c6842d1a5 Mon Sep 17 00:00:00 2001 From: Marc Fehr Date: Thu, 15 Oct 2026 11:18:31 +0200 Subject: [PATCH 1/4] fix: to_global for slices that cover only part of their mesh Slice.to_global() placed each subslice using the coordinates of its whole mesh, while the data of a slice defined with XB only covers its own range, so every such slice failed with a reshape ValueError (2D and 3D, node-based and cell-centred, with and without masked=True). Mesh.get_obstruction_mask_slice() likewise returned the mask of the whole mesh plane. A shared helper, extent_index_range(), now gives the nodes or cells a slice covers; SubSlice.get_coordinates(), the obstruction mask and to_global() all use it. The mask now has the shape of the subslice's data. Tested on an FDS 6.11.1 case (slcf_partial_fds611): each partial slice equals the matching section of a full slice in the same plane, with and without the mask. Co-Authored-By: Claude Opus 5.5 --- .gitignore | 1 + fdsreader/fds_classes/mesh.py | 31 ++++++----- fdsreader/slcf/slice.py | 23 ++++----- fdsreader/utils/extent.py | 17 +++++++ tests/cases/slcf_partial_fds611.tgz | 3 ++ tests/test_slice_partial_xb.py | 79 +++++++++++++++++++++++++++++ 6 files changed, 125 insertions(+), 29 deletions(-) create mode 100644 tests/cases/slcf_partial_fds611.tgz create mode 100644 tests/test_slice_partial_xb.py diff --git a/.gitignore b/.gitignore index 121e1343..0f9df37b 100644 --- a/.gitignore +++ b/.gitignore @@ -532,6 +532,7 @@ tests/cases/pl3d_data_fds*/ tests/cases/prof_data_fds*/ tests/cases/geomslice_data_fds*/ tests/cases/cface_data_fds*/ +tests/cases/slcf_partial_fds*/ # Pickle cache files *.pickle diff --git a/fdsreader/fds_classes/mesh.py b/fdsreader/fds_classes/mesh.py index bf7fc84c..d784b262 100644 --- a/fdsreader/fds_classes/mesh.py +++ b/fdsreader/fds_classes/mesh.py @@ -7,6 +7,7 @@ from fdsreader import settings from fdsreader.bndf import Boundary, Patch from fdsreader.utils import Dimension, Extent, Quantity +from fdsreader.utils.extent import extent_index_range class Mesh: @@ -74,25 +75,23 @@ def get_obstruction_mask(self, times: Sequence[float], cell_centered=False) -> n def get_obstruction_mask_slice(self, subslice): """Marks all cells of a single subslice which are blocked by an obstruction. - :returns: A 4-dimensional array with time as first and x,y,z as last dimensions. - The array depends on time as obstructions may be hidden at specific points in time. - For a 2D slice the axis of its orientation has length 1, for a 3D slice it is the mask of the - whole mesh. + :returns: A 4-dimensional array with time as first and x,y,z as last dimensions, matching the shape of the + subslice's data. The array depends on time as obstructions may be hidden at specific points in time. + For a 2D slice the axis of its orientation has length 1. """ orientation = subslice.orientation cell_centered = subslice.cell_centered - if orientation == 0: - # A 3D slice has no plane to cut out; indexing with orientation 0 would cut the time axis instead - return self.get_obstruction_mask(subslice.times, cell_centered=cell_centered) - value = subslice.extent[orientation][0] - - slc_index = self.coordinate_to_index((value,), dimension=(orientation,), cell_centered=cell_centered)[0] - - mask_indices = [slice(None)] * 4 - mask_indices[orientation] = slice(slc_index, slc_index + 1, 1) - mask_indices = tuple(mask_indices) - - return self.get_obstruction_mask(subslice.times, cell_centered=cell_centered)[mask_indices] + mask = self.get_obstruction_mask(subslice.times, cell_centered=cell_centered) + crop = [slice(None)] * 4 + if orientation != 0: # a 3D slice has no plane to cut out + value = subslice.extent[orientation][0] + slc_index = self.coordinate_to_index((value,), dimension=(orientation,), cell_centered=cell_centered)[0] + crop[orientation] = slice(slc_index, slc_index + 1) + # A slice defined with XB may cover only part of the mesh + for axis, dim in enumerate(("x", "y", "z"), start=1): + if axis != orientation: + crop[axis] = extent_index_range(self.coordinates[dim], subslice.extent[dim], cell_centered) + return mask[tuple(crop)] def coordinate_to_index( self, diff --git a/fdsreader/slcf/slice.py b/fdsreader/slcf/slice.py index 200c9bbd..20adddf0 100644 --- a/fdsreader/slcf/slice.py +++ b/fdsreader/slcf/slice.py @@ -12,6 +12,7 @@ from fdsreader.fds_classes import Mesh from fdsreader.utils import Dimension, Extent, NumpyArrayMixin, Quantity from fdsreader.utils.data import resample_mesh_axis +from fdsreader.utils.extent import extent_index_range _HANDLED_FUNCTIONS = {} @@ -53,20 +54,16 @@ def get_coordinates(self, ignore_cell_centered: bool = False) -> Dict[Literal["x :param ignore_cell_centered: Whether to shift the coordinates when the slice is cell_centered or not. """ - # orientation = ('x', 'y', 'z')[self.orientation - 1] if self.orientation != 0 else '' - # coords = {'x': set(), 'y': set(), 'z': set()} + cell_centered = self.cell_centered and not ignore_cell_centered coords: Dict[Literal["x", "y", "z"], np.ndarray] = {} for dim in ("x", "y", "z"): - co = self.mesh.coordinates[dim].copy() + co = self.mesh.coordinates[dim] # In case the slice is cell-centered, we will shift the coordinates by half a cell # and remove the last coordinate. Compute the per-cell offset before trimming, since a # mesh with only a single cell along this axis leaves just one coordinate afterwards. - if self.cell_centered and not ignore_cell_centered: + if cell_centered: co = co[:-1] + np.diff(co) / 2 - - coords[dim] = co[np.where(np.logical_and(co >= self.extent[dim][0], co <= self.extent[dim][1]))] - if coords[dim].size == 0: - coords[dim] = np.array([co[np.argmin(np.abs(co - self.extent[dim][0]))]]) + coords[dim] = co[extent_index_range(self.mesh.coordinates[dim], self.extent[dim], cell_centered)].copy() return coords @@ -591,14 +588,12 @@ def to_global( } step_sizes_max = {"x": 0, "y": 0, "z": 0} steps = dict() - global_max = {"x": -math.inf, "y": -math.inf, "z": -math.inf} for dim in ("x", "y", "z"): for sslc in subslices.values(): step_size = sslc.mesh.coordinates[dim][1] - sslc.mesh.coordinates[dim][0] step_sizes_min[dim] = min(step_size, step_sizes_min[dim]) step_sizes_max[dim] = max(step_size, step_sizes_max[dim]) - global_max[dim] = max(sslc.mesh.coordinates[dim][-1], global_max[dim]) for dim in ("x", "y", "z"): if step_sizes_min[dim] == 0: @@ -623,6 +618,8 @@ def place(slc: SubSlice): ) mask = slc.mesh.get_obstruction_mask_slice(slc) if masked else None + # The nodes the subslice covers: a slice defined with XB may cover only part of its mesh + nodes = slc.get_coordinates(ignore_cell_centered=True) start_idx = dict() end_idx = dict() inner_axes = list() @@ -637,8 +634,8 @@ def place(slc: SubSlice): 1, ) - start_idx[dim] = int(round((slc.mesh.coordinates[dim][0] - coord_min[dim]) / step_sizes_min[dim])) - end_idx[dim] = int(round((slc.mesh.coordinates[dim][-1] - coord_min[dim]) / step_sizes_min[dim])) + start_idx[dim] = int(round((nodes[dim][0] - coord_min[dim]) / step_sizes_min[dim])) + end_idx[dim] = int(round((nodes[dim][-1] - coord_min[dim]) / step_sizes_min[dim])) if self.cell_centered: if n_repeat > 1: @@ -648,7 +645,7 @@ def place(slc: SubSlice): continue end_idx[dim] += 1 - if not np.isclose(slc.mesh.coordinates[dim][-1], global_max[dim]): + if not np.isclose(nodes[dim][-1], coord_max[dim]): inner_axes.append(axis) slc_data = resample_mesh_axis(slc_data, axis + 1, n_repeat, keep_last=True) if masked: diff --git a/fdsreader/utils/extent.py b/fdsreader/utils/extent.py index 15f5bc22..75b4a5d8 100644 --- a/fdsreader/utils/extent.py +++ b/fdsreader/utils/extent.py @@ -1,5 +1,6 @@ from typing import List, Tuple, Union +import numpy as np from typing_extensions import Literal @@ -100,3 +101,19 @@ def as_list(self, reduced=True) -> List[float]: :param reduced: Whether to leave out empty extents or not. """ return list(self.as_tuple(reduced)) + + +def extent_index_range(coordinates: np.ndarray, bounds: Tuple[float, float], cell_centered=False) -> slice: + """Indices of the nodes (or, if cell_centered, the cells) within the bounds along one axis of a mesh. + If there are none, the index of the one nearest to the lower bound. + + :param coordinates: The node coordinates of the mesh along the axis. + :param bounds: Lower and upper bound along the axis, e.g. ``extent["x"]``. + """ + if cell_centered: + coordinates = coordinates[:-1] + np.diff(coordinates) / 2 + inside = np.flatnonzero((coordinates >= bounds[0]) & (coordinates <= bounds[1])) + if inside.size == 0: + nearest = int(np.argmin(np.abs(coordinates - bounds[0]))) + return slice(nearest, nearest + 1) + return slice(int(inside[0]), int(inside[-1]) + 1) diff --git a/tests/cases/slcf_partial_fds611.tgz b/tests/cases/slcf_partial_fds611.tgz new file mode 100644 index 00000000..7f4583eb --- /dev/null +++ b/tests/cases/slcf_partial_fds611.tgz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:850f5209ce860762eaf5ef35598499dcf9a5a653da58362666f2286dae3b22d8 +size 121734 diff --git a/tests/test_slice_partial_xb.py b/tests/test_slice_partial_xb.py new file mode 100644 index 00000000..6b0b8d75 --- /dev/null +++ b/tests/test_slice_partial_xb.py @@ -0,0 +1,79 @@ +"""Slices defined with XB that cover only part of their mesh (slcf_partial_fds611, FDS 6.11.1). + +The case has two meshes (x 0..1 and 1..2, 0.1 m cells), an obstruction at x 0.6..0.9, y 0.4..0.6, +z 0..0.5, and in the same planes full reference slices and slices limited with XB: 2D at z = 0.4 and +3D over x 0.3..1.6, y 0.2..0.8, z 0.2..0.7, each node-based and cell-centred, plus a 2D slice inside +one mesh only. FDS writes the same values into both, so each partial slice must equal the matching +section of the full one. +""" + +from pathlib import Path + +import numpy as np +import pytest + +from fdsreader import Simulation, settings + +CASES_DIR = Path(__file__).parent / "cases" +PARTIAL_IDS = ["z_part", "z_part_cc", "vol_part", "vol_part_cc", "z_part_m1"] + + +@pytest.fixture(scope="module") +def sim(): + original = settings.ENABLE_CACHING + settings.ENABLE_CACHING = False + try: + yield Simulation(str(CASES_DIR / "slcf_partial_fds611")) + finally: + settings.ENABLE_CACHING = original + + +def _full_slice(sim, partial): + return next( + slc + for slc in sim.slices + if not slc.id and slc.type == partial.type and slc.cell_centered == partial.cell_centered + ) + + +@pytest.mark.parametrize("masked", [False, True]) +@pytest.mark.parametrize("slice_id", PARTIAL_IDS) +def test_partial_slice_equals_section_of_full_slice(sim, slice_id, masked): + partial = sim.slices.get_by_id(slice_id) + full = _full_slice(sim, partial) + values, coords = partial.to_global(masked=masked, fill=np.nan, return_coordinates=True) + reference, ref_coords = full.to_global(masked=masked, fill=np.nan, return_coordinates=True) + + section = [] + for dim in ("x", "y", "z") if partial.type == "3D" else ("x", "y"): + first, last = (int(np.argmin(np.abs(ref_coords[dim] - c))) for c in (coords[dim][0], coords[dim][-1])) + section.append(slice(first, last + 1)) + n_t = min(len(values), len(reference)) # the 3D reference slices may have one frame more + np.testing.assert_array_equal(values[:n_t], reference[(slice(0, n_t),) + tuple(section)]) + if masked: + assert np.isnan(values).any() == np.isnan(reference[(slice(0, n_t),) + tuple(section)]).any() + + +def test_partial_slice_grid_and_coordinates(sim): + values, coords = sim.slices.get_by_id("z_part").to_global(return_coordinates=True) + assert values.shape[1:] == (14, 7) # nodes x 0.3..1.6, y 0.2..0.8 at 0.1 m + np.testing.assert_allclose([coords["x"][0], coords["x"][-1]], [0.3, 1.6], atol=1e-6) + np.testing.assert_allclose([coords["y"][0], coords["y"][-1]], [0.2, 0.8], atol=1e-6) + + +def test_masked_partial_slice_masks_the_obstruction(sim): + values, coords = sim.slices.get_by_id("z_part").to_global(masked=True, fill=np.nan, return_coordinates=True) + x, y = np.meshgrid(coords["x"], coords["y"], indexing="ij") + inside = (x > 0.6 + 1e-6) & (x < 0.9 - 1e-6) & (y > 0.4 + 1e-6) & (y < 0.6 - 1e-6) + assert inside.any() + assert np.isnan(values[:, inside]).all() + assert not np.isnan(values[:, (x > 1.0)]).any() # no obstruction in the second mesh + + +@pytest.mark.parametrize("slice_id", PARTIAL_IDS) +def test_obstruction_mask_matches_subslice_data(sim, slice_id): + for subslice in sim.slices.get_by_id(slice_id).subslices: + mask = subslice.mesh.get_obstruction_mask_slice(subslice) + data = subslice.data + expected = data.shape if subslice.orientation == 0 else np.expand_dims(data, subslice.orientation).shape + assert mask.shape == expected From a58b3d3850ec772d03ccba9a6f5660d45e046f90 Mon Sep 17 00:00:00 2001 From: Marc Fehr Date: Thu, 15 Oct 2026 11:24:53 +0200 Subject: [PATCH 2/4] fix: address review findings on partial XB slices - to_global takes coord_min/coord_max from the same subslice nodes it places the data with, instead of filtering the mesh coordinates again. - get_obstruction_mask takes an optional region, so the mask of a slice is built only for the cells it covers instead of the whole mesh. - SubSlice.get_coordinates computes the cell centres once. - The test case gets a third mesh with 0.2 m cells; the tests compare against the raw data of every mesh and locate sections by node index instead of the returned coordinates. Co-Authored-By: Claude Opus 5.5 --- fdsreader/fds_classes/mesh.py | 41 ++++++++------ fdsreader/slcf/slice.py | 27 +++++---- tests/cases/slcf_partial_fds611.tgz | 4 +- tests/test_slice_partial_xb.py | 85 ++++++++++++++++++++++------- 4 files changed, 105 insertions(+), 52 deletions(-) diff --git a/fdsreader/fds_classes/mesh.py b/fdsreader/fds_classes/mesh.py index d784b262..dfc6e857 100644 --- a/fdsreader/fds_classes/mesh.py +++ b/fdsreader/fds_classes/mesh.py @@ -49,27 +49,36 @@ def __init__( self._boundary_data: Dict[int, Boundary] = dict() - def get_obstruction_mask(self, times: Sequence[float], cell_centered=False) -> np.ndarray: + def get_obstruction_mask( + self, times: Sequence[float], cell_centered=False, region: Tuple[slice, slice, slice] | None = None + ) -> np.ndarray: """Marks all cells which are blocked by an obstruction. :param times: All timesteps of the simulation. + :param region: Index ranges along x, y and z (with explicit start and stop) to which the mask is limited. + The whole mesh by default. :returns: A 4-dimensional array with time as first and x,y,z as last dimensions. The array depends on time as obstructions may be hidden as specific points in time. """ - shape = self.dimension.shape(cell_centered=cell_centered) - mask = np.ones((len(times), shape[0], shape[1], shape[2]), dtype=bool) + if region is None: + region = tuple(slice(0, n) for n in self.dimension.shape(cell_centered=cell_centered)) + mask = np.ones((len(times),) + tuple(r.stop - r.start for r in region), dtype=bool) c = 1 if cell_centered else 0 for obst in self.obstructions: subobst = obst[self] - x1, x2 = subobst.bound_indices["x"] - y1, y2 = subobst.bound_indices["y"] - z1, z2 = subobst.bound_indices["z"] + blocked = list() + for dim, r in zip(("x", "y", "z"), region): + low, high = subobst.bound_indices[dim] + high = max(high + c, low + 1) + blocked.append(slice(max(low, r.start) - r.start, min(high, r.stop) - r.start)) + if any(b.stop <= b.start for b in blocked): + continue # the obstruction lies outside of the region t_idx = 0 for t in subobst.get_visible_times(times): while not np.isclose(t, times[t_idx]): t_idx += 1 - mask[t_idx, x1 : max(x2 + c, x1 + 1), y1 : max(y2 + c, y1 + 1), z1 : max(z2 + c, z1 + 1)] = False + mask[(t_idx,) + tuple(blocked)] = False return mask def get_obstruction_mask_slice(self, subslice): @@ -81,17 +90,15 @@ def get_obstruction_mask_slice(self, subslice): """ orientation = subslice.orientation cell_centered = subslice.cell_centered - mask = self.get_obstruction_mask(subslice.times, cell_centered=cell_centered) - crop = [slice(None)] * 4 - if orientation != 0: # a 3D slice has no plane to cut out - value = subslice.extent[orientation][0] - slc_index = self.coordinate_to_index((value,), dimension=(orientation,), cell_centered=cell_centered)[0] - crop[orientation] = slice(slc_index, slc_index + 1) - # A slice defined with XB may cover only part of the mesh + region = list() for axis, dim in enumerate(("x", "y", "z"), start=1): - if axis != orientation: - crop[axis] = extent_index_range(self.coordinates[dim], subslice.extent[dim], cell_centered) - return mask[tuple(crop)] + if axis == orientation: + value = subslice.extent[orientation][0] + slc_index = self.coordinate_to_index((value,), dimension=(orientation,), cell_centered=cell_centered)[0] + region.append(slice(slc_index, slc_index + 1)) + else: # a slice defined with XB may cover only part of the mesh + region.append(extent_index_range(self.coordinates[dim], subslice.extent[dim], cell_centered)) + return self.get_obstruction_mask(subslice.times, cell_centered=cell_centered, region=tuple(region)) def coordinate_to_index( self, diff --git a/fdsreader/slcf/slice.py b/fdsreader/slcf/slice.py index 20adddf0..07dab8ef 100644 --- a/fdsreader/slcf/slice.py +++ b/fdsreader/slcf/slice.py @@ -63,7 +63,7 @@ def get_coordinates(self, ignore_cell_centered: bool = False) -> Dict[Literal["x # mesh with only a single cell along this axis leaves just one coordinate afterwards. if cell_centered: co = co[:-1] + np.diff(co) / 2 - coords[dim] = co[extent_index_range(self.mesh.coordinates[dim], self.extent[dim], cell_centered)].copy() + coords[dim] = co[extent_index_range(co, self.extent[dim])].copy() return coords @@ -571,12 +571,12 @@ def to_global( for subslices in subslice_sets: coord_min = {"x": math.inf, "y": math.inf, "z": math.inf} coord_max = {"x": -math.inf, "y": -math.inf, "z": -math.inf} + # The nodes each subslice covers: a slice defined with XB may cover only part of its mesh + nodes = {mesh: slc.get_coordinates(ignore_cell_centered=True) for mesh, slc in subslices.items()} for dim in ("x", "y", "z"): - for slc in subslices.values(): - co = slc.mesh.coordinates[dim] - co = co[np.where(np.logical_and(co >= slc.extent[dim][0], co <= slc.extent[dim][1]))] - coord_min[dim] = min(co[0], coord_min[dim]) - coord_max[dim] = max(co[-1], coord_max[dim]) + for mesh_nodes in nodes.values(): + coord_min[dim] = min(mesh_nodes[dim][0], coord_min[dim]) + coord_max[dim] = max(mesh_nodes[dim][-1], coord_max[dim]) # The global grid will use the finest mesh as base and duplicate values of the coarser # meshes. Therefore, we first find the finest mesh and calculate the step size in each @@ -606,7 +606,7 @@ def to_global( grid = np.full((self.n_t, steps["x"], steps["y"], steps["z"]), np.nan) - def place(slc: SubSlice): + def place(mesh: str, slc: SubSlice): """Maps a subslice onto the global grid, keeping the node rows on all of its upper faces. :returns: The data (masked if requested), the obstruction mask (None if not masked), the region of @@ -618,8 +618,7 @@ def place(slc: SubSlice): ) mask = slc.mesh.get_obstruction_mask_slice(slc) if masked else None - # The nodes the subslice covers: a slice defined with XB may cover only part of its mesh - nodes = slc.get_coordinates(ignore_cell_centered=True) + slc_nodes = nodes[mesh] start_idx = dict() end_idx = dict() inner_axes = list() @@ -634,8 +633,8 @@ def place(slc: SubSlice): 1, ) - start_idx[dim] = int(round((nodes[dim][0] - coord_min[dim]) / step_sizes_min[dim])) - end_idx[dim] = int(round((nodes[dim][-1] - coord_min[dim]) / step_sizes_min[dim])) + start_idx[dim] = int(round((slc_nodes[dim][0] - coord_min[dim]) / step_sizes_min[dim])) + end_idx[dim] = int(round((slc_nodes[dim][-1] - coord_min[dim]) / step_sizes_min[dim])) if self.cell_centered: if n_repeat > 1: @@ -645,7 +644,7 @@ def place(slc: SubSlice): continue end_idx[dim] += 1 - if not np.isclose(nodes[dim][-1], coord_max[dim]): + if not np.isclose(slc_nodes[dim][-1], coord_max[dim]): inner_axes.append(axis) slc_data = resample_mesh_axis(slc_data, axis + 1, n_repeat, keep_last=True) if masked: @@ -664,8 +663,8 @@ def place(slc: SubSlice): # simulation space, which were left to a neighbouring mesh (node-based slices only) written = None if self.cell_centered else np.zeros(grid.shape[1:], dtype=bool) inner_faces = list() - for slc in subslices.values(): - slc_data, mask, region, inner_axes = place(slc) + for mesh, slc in subslices.items(): + slc_data, mask, region, inner_axes = place(mesh, slc) # The node row on an upper face inside the simulation space coincides with the first row of the # neighbouring mesh, which supplies it, so it is left out here. A copy of it is kept for the case diff --git a/tests/cases/slcf_partial_fds611.tgz b/tests/cases/slcf_partial_fds611.tgz index 7f4583eb..8b81e057 100644 --- a/tests/cases/slcf_partial_fds611.tgz +++ b/tests/cases/slcf_partial_fds611.tgz @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:850f5209ce860762eaf5ef35598499dcf9a5a653da58362666f2286dae3b22d8 -size 121734 +oid sha256:3ccb804ee03100d695ce9893b1f08f37d3c9ab10b9de98ae3294b109dd0dbeb9 +size 133762 diff --git a/tests/test_slice_partial_xb.py b/tests/test_slice_partial_xb.py index 6b0b8d75..af424f68 100644 --- a/tests/test_slice_partial_xb.py +++ b/tests/test_slice_partial_xb.py @@ -1,10 +1,10 @@ """Slices defined with XB that cover only part of their mesh (slcf_partial_fds611, FDS 6.11.1). -The case has two meshes (x 0..1 and 1..2, 0.1 m cells), an obstruction at x 0.6..0.9, y 0.4..0.6, -z 0..0.5, and in the same planes full reference slices and slices limited with XB: 2D at z = 0.4 and -3D over x 0.3..1.6, y 0.2..0.8, z 0.2..0.7, each node-based and cell-centred, plus a 2D slice inside -one mesh only. FDS writes the same values into both, so each partial slice must equal the matching -section of the full one. +The case has three meshes along x: M1 (0..1) and M2 (1..2) with 0.1 m cells and M3 (2..3) with 0.2 m cells, +obstructions at x 0.6..0.9 and 2.4..2.8 (y 0.4..0.6), and in the same planes full reference slices and slices +limited with XB: 2D at z = 0.4 and 3D over x 0.3..2.6, y 0.2..0.8, z 0.2..0.6, each node-based and +cell-centred, plus a 2D slice inside M1 only. FDS writes the same values into both, so each partial slice must +equal the matching section of the full one, and the raw data of each mesh. """ from pathlib import Path @@ -36,28 +36,67 @@ def _full_slice(sim, partial): ) +def _grids(slc, masked): + """The global grids of a slice with the subslices placed on each; cell-centred 2D slices give one grid per + plane, here the 0.1 m meshes and the 0.2 m mesh, which FDS writes at different heights.""" + grids = slc.to_global(masked=masked, fill=np.nan) + if not isinstance(grids, tuple): + return [(grids, slc.subslices)] + fine = [sub for sub in slc.subslices if sub.mesh.id != "M3"] + coarse = [sub for sub in slc.subslices if sub.mesh.id == "M3"] + return [(grids[0], fine), (grids[1], coarse)] + + +def _origin_and_step(subslices, dim): + nodes = [sub.get_coordinates(ignore_cell_centered=True)[dim] for sub in subslices] + steps = [sub.mesh.coordinates[dim][1] - sub.mesh.coordinates[dim][0] for sub in subslices] + return min(float(n[0]) for n in nodes), min(steps) + + @pytest.mark.parametrize("masked", [False, True]) @pytest.mark.parametrize("slice_id", PARTIAL_IDS) def test_partial_slice_equals_section_of_full_slice(sim, slice_id, masked): partial = sim.slices.get_by_id(slice_id) full = _full_slice(sim, partial) - values, coords = partial.to_global(masked=masked, fill=np.nan, return_coordinates=True) - reference, ref_coords = full.to_global(masked=masked, fill=np.nan, return_coordinates=True) - - section = [] - for dim in ("x", "y", "z") if partial.type == "3D" else ("x", "y"): - first, last = (int(np.argmin(np.abs(ref_coords[dim] - c))) for c in (coords[dim][0], coords[dim][-1])) - section.append(slice(first, last + 1)) - n_t = min(len(values), len(reference)) # the 3D reference slices may have one frame more - np.testing.assert_array_equal(values[:n_t], reference[(slice(0, n_t),) + tuple(section)]) - if masked: - assert np.isnan(values).any() == np.isnan(reference[(slice(0, n_t),) + tuple(section)]).any() + dims = ("x", "y", "z") if partial.type == "3D" else ("x", "y") + unmasked = _grids(partial, masked=False) + for (values, partial_subs), (reference, full_subs), (written, _) in zip( + _grids(partial, masked), _grids(full, masked), unmasked + ): + section = [] + for axis, dim in enumerate(dims): + partial_origin, step = _origin_and_step(partial_subs, dim) + full_origin, _ = _origin_and_step(full_subs, dim) + start = int(round((partial_origin - full_origin) / step)) + section.append(slice(start, start + values.shape[axis + 1])) + reference = reference[(slice(None),) + tuple(section)] + # FDS may write the slice of a coarser mesh a little further than that of a finer one (vol_part_cc: + # y up to 1.0 in M3, 0.9 in M1 and M2), so the partial grid has gaps the full slice does not have + covered = ~np.isnan(written) + assert covered.mean() > 0.8 + np.testing.assert_array_equal(values[covered], reference[covered]) + + +@pytest.mark.parametrize("slice_id", ["z_part", "vol_part"]) +def test_partial_slice_holds_the_raw_data_of_every_mesh(sim, slice_id): + slc = sim.slices.get_by_id(slice_id) + grid, coords = slc.to_global(return_coordinates=True) + dims = ("x", "y", "z") if slc.type == "3D" else ("x", "y") + for sub in slc.subslices: + sub_coords = sub.get_coordinates() + index = [np.searchsorted(coords[dim], sub_coords[dim] - 1e-6) for dim in dims] + values = grid[(slice(None),) + np.ix_(*index)] + if sub.mesh.id != "M3": # the node row at the upper x face is supplied by the next mesh + values, data = values[:, :-1], sub.data[:, :-1] + else: + data = sub.data + np.testing.assert_array_equal(values, data) def test_partial_slice_grid_and_coordinates(sim): values, coords = sim.slices.get_by_id("z_part").to_global(return_coordinates=True) - assert values.shape[1:] == (14, 7) # nodes x 0.3..1.6, y 0.2..0.8 at 0.1 m - np.testing.assert_allclose([coords["x"][0], coords["x"][-1]], [0.3, 1.6], atol=1e-6) + assert values.shape[1:] == (24, 7) # nodes x 0.3..2.6, y 0.2..0.8 on the 0.1 m grid + np.testing.assert_allclose([coords["x"][0], coords["x"][-1]], [0.3, 2.6], atol=1e-6) np.testing.assert_allclose([coords["y"][0], coords["y"][-1]], [0.2, 0.8], atol=1e-6) @@ -67,7 +106,7 @@ def test_masked_partial_slice_masks_the_obstruction(sim): inside = (x > 0.6 + 1e-6) & (x < 0.9 - 1e-6) & (y > 0.4 + 1e-6) & (y < 0.6 - 1e-6) assert inside.any() assert np.isnan(values[:, inside]).all() - assert not np.isnan(values[:, (x > 1.0)]).any() # no obstruction in the second mesh + assert not np.isnan(values[:, (x > 1.0) & (x < 2.0)]).any() # no obstruction in M2 @pytest.mark.parametrize("slice_id", PARTIAL_IDS) @@ -77,3 +116,11 @@ def test_obstruction_mask_matches_subslice_data(sim, slice_id): data = subslice.data expected = data.shape if subslice.orientation == 0 else np.expand_dims(data, subslice.orientation).shape assert mask.shape == expected + + +def test_obstruction_mask_of_a_region_equals_the_cropped_mask_of_the_mesh(sim): + for subslice in sim.slices.get_by_id("vol_part").subslices: + mesh = subslice.mesh + full = mesh.get_obstruction_mask(subslice.times) + region = (slice(1, 4), slice(2, 5), slice(0, 3)) + np.testing.assert_array_equal(mesh.get_obstruction_mask(subslice.times, region=region), full[:, 1:4, 2:5, 0:3]) From e3acdf82bb47cdd6d1bd1717f8e42638c0ab837f Mon Sep 17 00:00:00 2001 From: Marc Fehr Date: Thu, 15 Oct 2026 11:33:00 +0200 Subject: [PATCH 3/4] test: check partial XB slices against each new FDS release The FDS compatibility workflow runs tests/cases/fds_inputs/*.fds through the newest FDS release. The input of the partial XB case is added there, and tests/test_slice_partial_xb.py also runs against that fresh output when FDS_FRESH_OUTPUT_DIR is set, so a change in how FDS writes partial slices shows up when a new release comes out. Co-Authored-By: Claude Opus 5.5 --- tests/cases/fds_inputs/slcf_partial.fds | 28 +++++++++++++++++++++++++ tests/test_slice_partial_xb.py | 17 +++++++++++---- 2 files changed, 41 insertions(+), 4 deletions(-) create mode 100644 tests/cases/fds_inputs/slcf_partial.fds diff --git a/tests/cases/fds_inputs/slcf_partial.fds b/tests/cases/fds_inputs/slcf_partial.fds new file mode 100644 index 00000000..4f52728c --- /dev/null +++ b/tests/cases/fds_inputs/slcf_partial.fds @@ -0,0 +1,28 @@ +&HEAD CHID='partial', TITLE='Partial XB slices for fdsreader' / +&MESH ID='M1', IJK=10,10,10, XB=0.0,1.0,0.0,1.0,0.0,1.0 / +&MESH ID='M2', IJK=10,10,10, XB=1.0,2.0,0.0,1.0,0.0,1.0 / +&MESH ID='M3', IJK=5,5,5, XB=2.0,3.0,0.0,1.0,0.0,1.0 / +&TIME T_END=2.0 / +&DUMP DT_SLCF=0.5 / +&REAC FUEL='PROPANE', SOOT_YIELD=0.01 / +&SURF ID='FIRE', HRRPUA=500.0 / +&OBST ID='Box', XB=0.6,0.9,0.4,0.6,0.0,0.5 / +&OBST ID='Coarse', XB=2.4,2.8,0.4,0.6,0.0,0.6 / +&VENT XB=1.2,1.4,0.4,0.6,0.0,0.0, SURF_ID='FIRE' / +&VENT MB='XMIN', SURF_ID='OPEN' / &VENT MB='XMAX', SURF_ID='OPEN' / +&VENT MB='ZMAX', SURF_ID='OPEN' / + +full reference slices +&SLCF PBZ=0.4, QUANTITY='TEMPERATURE' / +&SLCF PBZ=0.4, QUANTITY='TEMPERATURE', CELL_CENTERED=T / +&SLCF XB=0.0,3.0,0.0,1.0,0.0,1.0, QUANTITY='TEMPERATURE' / +&SLCF XB=0.0,3.0,0.0,1.0,0.0,1.0, QUANTITY='TEMPERATURE', CELL_CENTERED=T / + +partial slices over all three meshes (M3 has 0.2 m cells) +&SLCF XB=0.3,2.6,0.2,0.8,0.4,0.4, QUANTITY='TEMPERATURE', ID='z_part' / +&SLCF XB=0.3,2.6,0.2,0.8,0.4,0.4, QUANTITY='TEMPERATURE', CELL_CENTERED=T, ID='z_part_cc' / +&SLCF XB=0.3,2.6,0.2,0.8,0.2,0.6, QUANTITY='TEMPERATURE', ID='vol_part' / +&SLCF XB=0.3,2.6,0.2,0.8,0.2,0.6, QUANTITY='TEMPERATURE', CELL_CENTERED=T, ID='vol_part_cc' / +partial slice in one mesh only +&SLCF XB=0.2,0.8,0.0,1.0,0.4,0.4, QUANTITY='TEMPERATURE', ID='z_part_m1' / +&TAIL / diff --git a/tests/test_slice_partial_xb.py b/tests/test_slice_partial_xb.py index af424f68..0ef4ecdc 100644 --- a/tests/test_slice_partial_xb.py +++ b/tests/test_slice_partial_xb.py @@ -1,4 +1,5 @@ -"""Slices defined with XB that cover only part of their mesh (slcf_partial_fds611, FDS 6.11.1). +"""Slices defined with XB that cover only part of their mesh (slcf_partial_fds611, FDS 6.11.1; input in +tests/cases/fds_inputs/slcf_partial.fds). The case has three meshes along x: M1 (0..1) and M2 (1..2) with 0.1 m cells and M3 (2..3) with 0.2 m cells, obstructions at x 0.6..0.9 and 2.4..2.8 (y 0.4..0.6), and in the same planes full reference slices and slices @@ -7,6 +8,7 @@ equal the matching section of the full one, and the raw data of each mesh. """ +import os from pathlib import Path import numpy as np @@ -17,13 +19,20 @@ CASES_DIR = Path(__file__).parent / "cases" PARTIAL_IDS = ["z_part", "z_part_cc", "vol_part", "vol_part_cc", "z_part_m1"] +# The FDS compatibility workflow runs tests/cases/fds_inputs/slcf_partial.fds through the newest FDS release +# and points FDS_FRESH_OUTPUT_DIR at the output, so the same checks run against that release too. +_CASES = [CASES_DIR / "slcf_partial_fds611"] +_FRESH = Path(os.environ.get("FDS_FRESH_OUTPUT_DIR") or "/nonexistent") / "slcf_partial" +if _FRESH.is_dir(): + _CASES.append(_FRESH) -@pytest.fixture(scope="module") -def sim(): + +@pytest.fixture(scope="module", params=_CASES, ids=lambda case: "fresh" if case == _FRESH else "fds611") +def sim(request): original = settings.ENABLE_CACHING settings.ENABLE_CACHING = False try: - yield Simulation(str(CASES_DIR / "slcf_partial_fds611")) + yield Simulation(str(request.param)) finally: settings.ENABLE_CACHING = original From acbfcba295addb7652576aa1c5c354cb427f3984 Mon Sep 17 00:00:00 2001 From: Marc Fehr Date: Thu, 15 Oct 2026 12:50:48 +0200 Subject: [PATCH 4/4] fix: address the second review of partial XB slices - The test case gets a coarse mesh before the fine ones as well (four meshes: 0.2, 0.1, 0.1, 0.2 m), so both orders are covered. - get_obstruction_mask resolves open ends and negative indices of the region like numpy and rejects steps. - Slice.get_coordinates drops its own nearest-node fallback: every subslice now provides at least one coordinate per axis. - Tests: the region test cuts through each mesh's obstruction and checks that cells are masked, the section comparison checks that partial and full slices give the same number of grids, and the fallback test runs without patching SubSlice.get_coordinates. Co-Authored-By: Claude Opus 5.5 --- fdsreader/fds_classes/mesh.py | 10 +++- fdsreader/slcf/slice.py | 20 +------- tests/cases/fds_inputs/slcf_partial.fds | 16 +++--- tests/cases/slcf_partial_fds611.tgz | 4 +- tests/test_slice_coordinates.py | 30 ++++------- tests/test_slice_partial_xb.py | 66 +++++++++++++++++-------- 6 files changed, 75 insertions(+), 71 deletions(-) diff --git a/fdsreader/fds_classes/mesh.py b/fdsreader/fds_classes/mesh.py index dfc6e857..3dbd41d5 100644 --- a/fdsreader/fds_classes/mesh.py +++ b/fdsreader/fds_classes/mesh.py @@ -55,13 +55,19 @@ def get_obstruction_mask( """Marks all cells which are blocked by an obstruction. :param times: All timesteps of the simulation. - :param region: Index ranges along x, y and z (with explicit start and stop) to which the mask is limited. + :param region: Index ranges along x, y and z, as slices without a step, to which the mask is limited. The whole mesh by default. :returns: A 4-dimensional array with time as first and x,y,z as last dimensions. The array depends on time as obstructions may be hidden as specific points in time. """ + shape = self.dimension.shape(cell_centered=cell_centered) if region is None: - region = tuple(slice(0, n) for n in self.dimension.shape(cell_centered=cell_centered)) + region = (slice(None),) * 3 + if any(r.step not in (None, 1) for r in region): + raise ValueError(f"The region must not have a step, got {region}") + # Resolve open ends and negative indices like numpy indexing does + region = tuple(slice(*r.indices(n)[:2]) for r, n in zip(region, shape)) + region = tuple(slice(r.start, max(r.start, r.stop)) for r in region) mask = np.ones((len(times),) + tuple(r.stop - r.start for r in region), dtype=bool) c = 1 if cell_centered else 0 diff --git a/fdsreader/slcf/slice.py b/fdsreader/slcf/slice.py index 07dab8ef..1bef95da 100644 --- a/fdsreader/slcf/slice.py +++ b/fdsreader/slcf/slice.py @@ -378,24 +378,8 @@ def get_coordinates(self, ignore_cell_centered: bool = False) -> Dict[Literal["x # (sometimes only 6 decimal places, e.g. 1.0 and 1.000001) diff = coords[dim][1:] - coords[dim][:-1] coords[dim] = np.concatenate((coords[dim][:1], coords[dim][1:][diff > 0.000002]), axis=0) - - if len(coords[dim]) == 0: - slice_coordinate = self.extent[dim][0] - nearest_coordinate = np.inf - for subslice in self._subslices.values(): - mesh_coords = subslice.mesh.coordinates[dim] - idx = np.searchsorted(mesh_coords, slice_coordinate, side="left") - if idx > 0 and ( - idx == mesh_coords.size - or math.fabs(slice_coordinate - mesh_coords[idx - 1]) - < math.fabs(slice_coordinate - mesh_coords[idx]) - ): - idx = idx - 1 - if math.fabs(mesh_coords[idx] - slice_coordinate) < math.fabs( - nearest_coordinate - slice_coordinate - ): - nearest_coordinate = mesh_coords[idx] - coords[dim] = np.array([nearest_coordinate]) + # Every subslice has at least one coordinate on each axis (SubSlice.get_coordinates falls back to the + # nearest one), so this is only empty for a slice without subslices return coords diff --git a/tests/cases/fds_inputs/slcf_partial.fds b/tests/cases/fds_inputs/slcf_partial.fds index 4f52728c..eb0304f7 100644 --- a/tests/cases/fds_inputs/slcf_partial.fds +++ b/tests/cases/fds_inputs/slcf_partial.fds @@ -1,4 +1,5 @@ &HEAD CHID='partial', TITLE='Partial XB slices for fdsreader' / +&MESH ID='M0', IJK=5,5,5, XB=-1.0,0.0,0.0,1.0,0.0,1.0 / &MESH ID='M1', IJK=10,10,10, XB=0.0,1.0,0.0,1.0,0.0,1.0 / &MESH ID='M2', IJK=10,10,10, XB=1.0,2.0,0.0,1.0,0.0,1.0 / &MESH ID='M3', IJK=5,5,5, XB=2.0,3.0,0.0,1.0,0.0,1.0 / @@ -9,20 +10,21 @@ &OBST ID='Box', XB=0.6,0.9,0.4,0.6,0.0,0.5 / &OBST ID='Coarse', XB=2.4,2.8,0.4,0.6,0.0,0.6 / &VENT XB=1.2,1.4,0.4,0.6,0.0,0.0, SURF_ID='FIRE' / +&OBST ID='Front', XB=-0.6,-0.2,0.4,0.6,0.0,0.6 / &VENT MB='XMIN', SURF_ID='OPEN' / &VENT MB='XMAX', SURF_ID='OPEN' / &VENT MB='ZMAX', SURF_ID='OPEN' / full reference slices &SLCF PBZ=0.4, QUANTITY='TEMPERATURE' / &SLCF PBZ=0.4, QUANTITY='TEMPERATURE', CELL_CENTERED=T / -&SLCF XB=0.0,3.0,0.0,1.0,0.0,1.0, QUANTITY='TEMPERATURE' / -&SLCF XB=0.0,3.0,0.0,1.0,0.0,1.0, QUANTITY='TEMPERATURE', CELL_CENTERED=T / +&SLCF XB=-1.0,3.0,0.0,1.0,0.0,1.0, QUANTITY='TEMPERATURE' / +&SLCF XB=-1.0,3.0,0.0,1.0,0.0,1.0, QUANTITY='TEMPERATURE', CELL_CENTERED=T / -partial slices over all three meshes (M3 has 0.2 m cells) -&SLCF XB=0.3,2.6,0.2,0.8,0.4,0.4, QUANTITY='TEMPERATURE', ID='z_part' / -&SLCF XB=0.3,2.6,0.2,0.8,0.4,0.4, QUANTITY='TEMPERATURE', CELL_CENTERED=T, ID='z_part_cc' / -&SLCF XB=0.3,2.6,0.2,0.8,0.2,0.6, QUANTITY='TEMPERATURE', ID='vol_part' / -&SLCF XB=0.3,2.6,0.2,0.8,0.2,0.6, QUANTITY='TEMPERATURE', CELL_CENTERED=T, ID='vol_part_cc' / +partial slices over all four meshes (M0 and M3 have 0.2 m cells) +&SLCF XB=-0.6,2.6,0.2,0.8,0.4,0.4, QUANTITY='TEMPERATURE', ID='z_part' / +&SLCF XB=-0.6,2.6,0.2,0.8,0.4,0.4, QUANTITY='TEMPERATURE', CELL_CENTERED=T, ID='z_part_cc' / +&SLCF XB=-0.6,2.6,0.2,0.8,0.2,0.6, QUANTITY='TEMPERATURE', ID='vol_part' / +&SLCF XB=-0.6,2.6,0.2,0.8,0.2,0.6, QUANTITY='TEMPERATURE', CELL_CENTERED=T, ID='vol_part_cc' / partial slice in one mesh only &SLCF XB=0.2,0.8,0.0,1.0,0.4,0.4, QUANTITY='TEMPERATURE', ID='z_part_m1' / &TAIL / diff --git a/tests/cases/slcf_partial_fds611.tgz b/tests/cases/slcf_partial_fds611.tgz index 8b81e057..90946827 100644 --- a/tests/cases/slcf_partial_fds611.tgz +++ b/tests/cases/slcf_partial_fds611.tgz @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:3ccb804ee03100d695ce9893b1f08f37d3c9ab10b9de98ae3294b109dd0dbeb9 -size 133762 +oid sha256:b4c222a4e3856af40bedc78562c55b708c62fb61db01e7ca6445820e7e37d554 +size 145376 diff --git a/tests/test_slice_coordinates.py b/tests/test_slice_coordinates.py index 16c7c866..7cd193c9 100644 --- a/tests/test_slice_coordinates.py +++ b/tests/test_slice_coordinates.py @@ -63,17 +63,14 @@ def test_slice_get_coordinates_cell_centered_oriented(): assert coords["y"] == pytest.approx([1.95]) -def test_slice_get_coordinates_empty_dim_fallback_uses_mesh_object(monkeypatch): - """Regression test: the "no coordinates found" fallback in Slice.get_coordinates() - (fdsreader/slcf/slice.py) had the same self._subslices.keys()/values() mix-up as - GitHub issue #118 and would raise AttributeError: 'str' object has no attribute - 'coordinates' if it were ever triggered, e.g. if a subslice's coordinates for a - dimension ended up not overlapping the slice's extent at all. - - The extent's x-coordinate (5.0) is placed beyond the mesh's last x-coordinate - (3.0) on purpose: this drives np.searchsorted's index to mesh_coords.size, which - also exercised a second, independent bug in the same fallback (it incremented the - index instead of decrementing it, causing an IndexError one past the array end).""" +def test_slice_get_coordinates_extent_outside_the_mesh_falls_back_to_nearest_node(): + """Regression test: a subslice whose extent does not overlap its mesh's coordinates along an axis used to + leave Slice.get_coordinates() without coordinates there, and its own fallback had the same + self._subslices.keys()/values() mix-up as GitHub issue #118 and an IndexError one past the array end. + SubSlice.get_coordinates() now falls back to the nearest node itself (extent_index_range), so the slice + gets that node. + + The extent's x-coordinate (5.0) lies beyond the mesh's last x-coordinate (3.0) on purpose.""" mesh = SimpleNamespace( id="mesh1", coordinates={ @@ -87,16 +84,7 @@ def test_slice_get_coordinates_empty_dim_fallback_uses_mesh_object(monkeypatch): slc.cell_centered = False slc.orientation = 0 # 3D slice: every dimension goes through the general branch slc.extent = Extent(5, 5, 0, 4, 0, 3) - subslice = SubSlice(slc, "", None, slc.extent, mesh) - slc._subslices = {mesh.id: subslice} - - # Simulate the "no overlapping coordinate found" case for the x dimension, which - # is what forces Slice.get_coordinates() into its buggy fallback branch. - monkeypatch.setattr( - subslice, - "get_coordinates", - lambda ignore_cell_centered=False: {"x": np.array([]), "y": np.array([2.0]), "z": np.array([1.0])}, - ) + slc._subslices = {mesh.id: SubSlice(slc, "", None, slc.extent, mesh)} coords = slc.get_coordinates() # The nearest available mesh coordinate is the last one (3.0), not an out-of-bounds index. diff --git a/tests/test_slice_partial_xb.py b/tests/test_slice_partial_xb.py index 0ef4ecdc..881744e6 100644 --- a/tests/test_slice_partial_xb.py +++ b/tests/test_slice_partial_xb.py @@ -1,11 +1,12 @@ """Slices defined with XB that cover only part of their mesh (slcf_partial_fds611, FDS 6.11.1; input in tests/cases/fds_inputs/slcf_partial.fds). -The case has three meshes along x: M1 (0..1) and M2 (1..2) with 0.1 m cells and M3 (2..3) with 0.2 m cells, -obstructions at x 0.6..0.9 and 2.4..2.8 (y 0.4..0.6), and in the same planes full reference slices and slices -limited with XB: 2D at z = 0.4 and 3D over x 0.3..2.6, y 0.2..0.8, z 0.2..0.6, each node-based and -cell-centred, plus a 2D slice inside M1 only. FDS writes the same values into both, so each partial slice must -equal the matching section of the full one, and the raw data of each mesh. +The case has four meshes along x: M0 (-1..0) and M3 (2..3) with 0.2 m cells, M1 (0..1) and M2 (1..2) with +0.1 m cells, so a coarse mesh lies both before and after the fine ones. Obstructions are at x -0.6..-0.2, +0.6..0.9 and 2.4..2.8 (y 0.4..0.6). In the same planes there are full reference slices and slices limited +with XB: 2D at z = 0.4 and 3D over x -0.6..2.6, y 0.2..0.8, z 0.2..0.6, each node-based and cell-centred, +plus a 2D slice inside M1 only. FDS writes the same values into both, so each partial slice must equal the +matching section of the full one, and the raw data of each mesh. """ import os @@ -46,14 +47,17 @@ def _full_slice(sim, partial): def _grids(slc, masked): - """The global grids of a slice with the subslices placed on each; cell-centred 2D slices give one grid per - plane, here the 0.1 m meshes and the 0.2 m mesh, which FDS writes at different heights.""" + """The global grids of a slice with the subslices placed on each. A cell-centred 2D slice gives one grid + per plane, here the 0.2 m and the 0.1 m meshes, which FDS writes at different heights; as in to_global, + the first grid holds the subslices in the plane of the first one.""" grids = slc.to_global(masked=masked, fill=np.nan) if not isinstance(grids, tuple): return [(grids, slc.subslices)] - fine = [sub for sub in slc.subslices if sub.mesh.id != "M3"] - coarse = [sub for sub in slc.subslices if sub.mesh.id == "M3"] - return [(grids[0], fine), (grids[1], coarse)] + dim = ("x", "y", "z")[slc.orientation - 1] + base = slc.subslices[0].get_coordinates()[dim][0] + first = [sub for sub in slc.subslices if np.isclose(sub.get_coordinates()[dim][0], base)] + second = [sub for sub in slc.subslices if not np.isclose(sub.get_coordinates()[dim][0], base)] + return [(grids[0], first), (grids[1], second)] def _origin_and_step(subslices, dim): @@ -68,10 +72,9 @@ def test_partial_slice_equals_section_of_full_slice(sim, slice_id, masked): partial = sim.slices.get_by_id(slice_id) full = _full_slice(sim, partial) dims = ("x", "y", "z") if partial.type == "3D" else ("x", "y") - unmasked = _grids(partial, masked=False) - for (values, partial_subs), (reference, full_subs), (written, _) in zip( - _grids(partial, masked), _grids(full, masked), unmasked - ): + partial_grids, full_grids, unmasked = _grids(partial, masked), _grids(full, masked), _grids(partial, False) + assert len(partial_grids) == len(full_grids) + for (values, partial_subs), (reference, full_subs), (written, _) in zip(partial_grids, full_grids, unmasked): section = [] for axis, dim in enumerate(dims): partial_origin, step = _origin_and_step(partial_subs, dim) @@ -79,10 +82,11 @@ def test_partial_slice_equals_section_of_full_slice(sim, slice_id, masked): start = int(round((partial_origin - full_origin) / step)) section.append(slice(start, start + values.shape[axis + 1])) reference = reference[(slice(None),) + tuple(section)] - # FDS may write the slice of a coarser mesh a little further than that of a finer one (vol_part_cc: - # y up to 1.0 in M3, 0.9 in M1 and M2), so the partial grid has gaps the full slice does not have + # The partial grid has gaps the full slice does not have: between M0 and M3 in the plane of the coarse + # meshes, and where FDS writes the slice of a coarser mesh a little further than that of a finer one + # (vol_part_cc: y up to 1.0 in M0 and M3, 0.9 in M1 and M2) covered = ~np.isnan(written) - assert covered.mean() > 0.8 + assert covered.any() np.testing.assert_array_equal(values[covered], reference[covered]) @@ -104,8 +108,8 @@ def test_partial_slice_holds_the_raw_data_of_every_mesh(sim, slice_id): def test_partial_slice_grid_and_coordinates(sim): values, coords = sim.slices.get_by_id("z_part").to_global(return_coordinates=True) - assert values.shape[1:] == (24, 7) # nodes x 0.3..2.6, y 0.2..0.8 on the 0.1 m grid - np.testing.assert_allclose([coords["x"][0], coords["x"][-1]], [0.3, 2.6], atol=1e-6) + assert values.shape[1:] == (33, 7) # nodes x -0.6..2.6, y 0.2..0.8 on the 0.1 m grid + np.testing.assert_allclose([coords["x"][0], coords["x"][-1]], [-0.6, 2.6], atol=1e-6) np.testing.assert_allclose([coords["y"][0], coords["y"][-1]], [0.2, 0.8], atol=1e-6) @@ -128,8 +132,28 @@ def test_obstruction_mask_matches_subslice_data(sim, slice_id): def test_obstruction_mask_of_a_region_equals_the_cropped_mask_of_the_mesh(sim): + meshes_with_obstruction = 0 for subslice in sim.slices.get_by_id("vol_part").subslices: mesh = subslice.mesh full = mesh.get_obstruction_mask(subslice.times) - region = (slice(1, 4), slice(2, 5), slice(0, 3)) - np.testing.assert_array_equal(mesh.get_obstruction_mask(subslice.times, region=region), full[:, 1:4, 2:5, 0:3]) + if not mesh.obstructions: + continue + # A region around the obstruction with a margin, so it does not start at index 0 + bounds = mesh.obstructions[0][mesh].bound_indices + region = tuple(slice(max(bounds[dim][0] - 1, 0), bounds[dim][1] + 2) for dim in ("x", "y", "z")) + assert region[0].start > 0 + mask = mesh.get_obstruction_mask(subslice.times, region=region) + np.testing.assert_array_equal(mask, full[(slice(None),) + region]) + assert (~mask).any() and mask.any() + meshes_with_obstruction += 1 + assert meshes_with_obstruction == 3 # M0, M1 and M3 + + +def test_obstruction_mask_region_resolves_open_ends_and_rejects_steps(sim): + mesh = sim.slices.get_by_id("vol_part").subslices[0].mesh + times = sim.slices.get_by_id("vol_part").times + full = mesh.get_obstruction_mask(times) + open_region = (slice(None, 3), slice(-2, None), slice(None)) + np.testing.assert_array_equal(mesh.get_obstruction_mask(times, region=open_region), full[:, :3, -2:, :]) + with pytest.raises(ValueError): + mesh.get_obstruction_mask(times, region=(slice(0, 4, 2), slice(None), slice(None)))