fix: to_global for slices that cover only part of their mesh

.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

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,

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

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)

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 /

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.