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..3dbd41d5 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: @@ -48,51 +49,62 @@ 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, 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) - mask = np.ones((len(times), shape[0], shape[1], shape[2]), dtype=bool) + if region is None: + 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 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): """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] + region = list() + for axis, dim in enumerate(("x", "y", "z"), start=1): + 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 200c9bbd..1bef95da 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(co, self.extent[dim])].copy() return coords @@ -381,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 @@ -574,12 +555,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 @@ -591,14 +572,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: @@ -611,7 +590,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 @@ -623,6 +602,7 @@ def place(slc: SubSlice): ) mask = slc.mesh.get_obstruction_mask_slice(slc) if masked else None + slc_nodes = nodes[mesh] start_idx = dict() end_idx = dict() inner_axes = list() @@ -637,8 +617,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((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: @@ -648,7 +628,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(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: @@ -667,8 +647,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/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/fds_inputs/slcf_partial.fds b/tests/cases/fds_inputs/slcf_partial.fds new file mode 100644 index 00000000..eb0304f7 --- /dev/null +++ b/tests/cases/fds_inputs/slcf_partial.fds @@ -0,0 +1,30 @@ +&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 / +&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' / +&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=-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 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 new file mode 100644 index 00000000..90946827 --- /dev/null +++ b/tests/cases/slcf_partial_fds611.tgz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +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 new file mode 100644 index 00000000..881744e6 --- /dev/null +++ b/tests/test_slice_partial_xb.py @@ -0,0 +1,159 @@ +"""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 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 +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"] + +# 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", 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(request.param)) + 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 + ) + + +def _grids(slc, masked): + """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)] + 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): + 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) + dims = ("x", "y", "z") if partial.type == "3D" else ("x", "y") + 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) + 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)] + # 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.any() + 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:] == (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) + + +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) & (x < 2.0)]).any() # no obstruction in M2 + + +@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 + + +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) + 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)))