refactor: address the review of #147

fdsreader/bndf/obstruction.py

@@ -324,24 +324,33 @@ def orientations(self):
         return []
 
     def get_coordinates(
-        self, ignore_cell_centered: bool = False
+        self, ignore_cell_centered: bool = False, quantity: Union[str, Quantity, None] = None
     ) -> Dict[int, Dict[Literal["x", "y", "z"], np.ndarray]]:
         """Returns a dictionary containing a numpy ndarray with coordinates for each dimension.
         For cell-centered boundary data, the coordinates can be adjusted to represent cell-centered coordinates.
 
         :param ignore_cell_centered: Whether to shift the coordinates when the bndf is cell_centered or not.
+        :param quantity: The quantity whose data the coordinates belong to; node-based and cell-centred quantities have
+            different coordinates. Defaults to the first quantity loaded.
         """
         if self.has_boundary_data:
-            return {
-                patch.orientation: patch.get_coordinates(ignore_cell_centered)
-                for patch in next(iter(self._boundary_data.values())).patches
-            }
+            boundary = next(iter(self._boundary_data.values())) if quantity is None else self.get_data(quantity)
+            return {patch.orientation: patch.get_coordinates(ignore_cell_centered) for patch in boundary.patches}
         return {}
 
-    def get_nearest_index(self, dimension: Literal["x", "y", "z"], orientation: int, value: float) -> int:
-        """Get the nearest mesh coordinate index in a specific dimension for a specific orientation."""
+    def get_nearest_index(
+        self,
+        dimension: Literal["x", "y", "z"],
+        orientation: int,
+        value: float,
+        quantity: Union[str, Quantity, None] = None,
+    ) -> int:
+        """Get the nearest mesh coordinate index in a specific dimension for a specific orientation.
+
+        :param quantity: The quantity whose data the index refers to, see :meth:`get_coordinates`.
+        """
         if self.has_boundary_data:
-            coords = self.get_coordinates()[orientation][dimension]
+            coords = self.get_coordinates(quantity=quantity)[orientation][dimension]
             idx = np.searchsorted(coords, value, side="left")
             if idx > 0 and (idx == coords.size or math.fabs(value - coords[idx - 1]) < math.fabs(value - coords[idx])):
                 return idx - 1
@@ -531,12 +540,14 @@ def get_visible_times(self, times: Sequence[float]):
         return np.array([])
 
     def get_coordinates(
-        self, ignore_cell_centered: bool = False
+        self, ignore_cell_centered: bool = False, quantity: Union[str, Quantity, None] = None
     ) -> Dict[int, Dict[Literal["x", "y", "z"], np.ndarray]]:
         """Returns a dictionary containing a numpy ndarray with coordinates for each dimension.
         For cell-centered boundary data, the coordinates can be adjusted to represent cell-centered coordinates.
 
         :param ignore_cell_centered: Whether to shift the coordinates when the bndf is cell_centered or not.
+        :param quantity: The quantity whose data the coordinates belong to; node-based and cell-centred quantities have
+            different coordinates. Defaults to the first quantity loaded for each part of the obstruction.
         """
         if self.has_boundary_data:
             all_coords = dict()
@@ -550,7 +561,10 @@ def get_coordinates(
                         continue
                     min_delta = math.inf
                     for subobst in self._all_subobstructions:
-                        subobst_coords = subobst.get_coordinates(ignore_cell_centered)
+                        try:
+                            subobst_coords = subobst.get_coordinates(ignore_cell_centered, quantity)
+                        except KeyError:  # this part has no data for the quantity
+                            continue
                         if orientation_int not in subobst_coords:
                             continue
                         co = subobst_coords[orientation_int][dim]
@@ -591,10 +605,19 @@ def get_coordinates(
             return all_coords
         return dict()
 
-    def get_nearest_index(self, dimension: Literal["x", "y", "z"], orientation: int, value: float) -> int:
-        """Get the nearest mesh coordinate index in a specific dimension for a specific orientation."""
+    def get_nearest_index(
+        self,
+        dimension: Literal["x", "y", "z"],
+        orientation: int,
+        value: float,
+        quantity: Union[str, Quantity, None] = None,
+    ) -> int:
+        """Get the nearest mesh coordinate index in a specific dimension for a specific orientation.
+
+        :param quantity: The quantity whose data the index refers to, see :meth:`get_coordinates`.
+        """
         if self.has_boundary_data:
-            coords = self.get_coordinates()[orientation][dimension]
+            coords = self.get_coordinates(quantity=quantity)[orientation][dimension]
             idx = np.searchsorted(coords, value, side="left")
             if idx > 0 and (idx == coords.size or math.fabs(value - coords[idx - 1]) < math.fabs(value - coords[idx])):
                 return idx - 1
@@ -734,8 +757,9 @@ def vmin(
 
         :param orientation: Optionally filter by patches with a specific orientation.
         """
-        values = [s.vmin(quantity, orientation) for s in self._all_subobstructions if s.has_boundary_data]
-        return min((v for v in values if not np.isnan(v)), default=np.nan)
+        return min(
+            (v for v in self._part_values(quantity, lambda b: b.vmin(orientation)) if not np.isnan(v)), default=np.nan
+        )
 
     def vmax(
         self,
@@ -746,8 +770,22 @@ def vmax(
 
         :param orientation: Optionally filter by patches with a specific orientation.
         """
-        values = [s.vmax(quantity, orientation) for s in self._all_subobstructions if s.has_boundary_data]
-        return max((v for v in values if not np.isnan(v)), default=np.nan)
+        return max(
+            (v for v in self._part_values(quantity, lambda b: b.vmax(orientation)) if not np.isnan(v)), default=np.nan
+        )
+
+    def _part_values(self, quantity: Union[str, Quantity], value: Callable[[Boundary], float]) -> List[float]:
+        """A value of the boundary data of each part that has the quantity. Raises KeyError if no part has it, but
+        the obstruction has boundary data."""
+        values = list()
+        for subobst in self._all_subobstructions:
+            try:
+                values.append(value(subobst.get_data(quantity)))
+            except KeyError:
+                continue
+        if not values and self.has_boundary_data:
+            raise KeyError(f"No boundary data found for quantity '{quantity}'")
+        return values
 
     def get_subobstructions(self, mesh: "Mesh") -> List[SubObstruction]:
         """All parts of this obstruction in the given mesh; several OBST lines can share an ID."""
@@ -810,6 +848,10 @@ def _global_boundary_arrays(
             patches_by_orientation.setdefault(patch.orientation, list()).append(patch)
     if not patches_by_orientation:
         raise KeyError(f"No boundary data found for quantity '{quantity}'")
+    if len({patch.cell_centered for patches in patches_by_orientation.values() for patch in patches}) > 1:
+        raise ValueError(
+            f"'{quantity}' is written node-based in some meshes and cell-centred in others, which cannot be combined"
+        )
 
     result = dict()
     for orientation in sorted(patches_by_orientation):

fdsreader/fds_classes/mesh.py

@@ -228,6 +228,11 @@ def _add_patches(
             for patch in patches:
                 patch._boundary_parent = self._boundary_data[bid]
 
+    @property
+    def has_boundary_data(self) -> bool:
+        """Whether boundary data has been output on the faces of this mesh, see :meth:`get_boundary_data`."""
+        return len(self._boundary_data) != 0
+
     def get_boundary_data(self, quantity: Union[str, Quantity]) -> Boundary:
         """The boundary data of a quantity on the faces of this mesh that are not faces of obstructions (those belong
         to :class:`Obstruction`), e.g. the floor of the simulation domain. With ``MB`` vents in a multi-mesh case (see

fdsreader/fds_classes/mesh_collection.py

@@ -31,8 +31,11 @@ def get_global_boundary_data_arrays(
         per mesh face; with ``MB`` vents in a multi-mesh case (see the FDS User Guide), it may also write faces
         between two meshes, which then appear as further positions.
 
-        :param quantity: The quantity's name or short name.
+        :param quantity: The quantity's name or short name. Raises KeyError if no mesh has data for it, and returns an
+            empty dict if there are no such data at all, like for an obstruction without boundary data.
         """
+        if not any(mesh.has_boundary_data for mesh in self):
+            return dict()
         return _global_boundary_arrays([mesh.get_boundary_data for mesh in self], quantity)
 
     def __repr__(self):

tests/cases/fds_inputs/bndf_domain.fds

@@ -7,8 +7,8 @@
 &SURF ID='FIRE', HRRPUA=500.0 /
 The vent splits the floor of both meshes into several patches
 &VENT XB=0.6,1.6,0.3,0.7,0.0,0.0, SURF_ID='FIRE' /
-&VENT MB='XMIN', SURF_ID='OPEN' /
-&VENT MB='ZMAX', SURF_ID='OPEN' /
+&VENT DB='XMIN', SURF_ID='OPEN' /
+&VENT DB='ZMAX', SURF_ID='OPEN' /
 An obstruction standing on the floor
 &OBST ID='Box', XB=0.2,0.4,0.6,0.8,0.0,0.2 /
 &BNDF QUANTITY='BURNING RATE', CELL_CENTERED=T /

tests/cases/fds_inputs/bndf_global.fds

@@ -6,8 +6,9 @@
 &REAC FUEL='PROPANE' /
 &SURF ID='FIRE', HRRPUA=500.0 /
 &VENT XB=0.1,0.5,0.1,0.9,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' /
+&VENT DB='XMIN', SURF_ID='OPEN' /
+&VENT DB='XMAX', SURF_ID='OPEN' /
+&VENT DB='ZMAX', SURF_ID='OPEN' /
 
 Across both meshes: the top at z = 0.53 snaps to 0.5 in 'Fine' and to 0.6 in 'Coarse'
 &OBST ID='Table', XB=0.6,1.4,0.2,0.8,0.0,0.53 /

tests/test_boundary_domain.py

@@ -10,7 +10,7 @@
 
 import numpy as np
 import pytest
-from fds_cases import case_dirs, case_id, load
+from fds_cases import CASES_DIR, case_dirs, case_id, load
 
 CELLS = "BURNING RATE"
 NODES = "WALL TEMPERATURE"
@@ -96,3 +96,9 @@ def test_obstruction_faces_are_not_included(sim):
 def test_unknown_quantity(sim):
     with pytest.raises(KeyError):
         sim.meshes.get_global_boundary_data_arrays("VELOCITY")
+
+
+def test_a_case_without_boundary_data_gives_no_arrays():
+    sim = load(CASES_DIR / "slcf_partial_fds611")
+    assert not any(mesh.has_boundary_data for mesh in sim.meshes)
+    assert sim.meshes.get_global_boundary_data_arrays("WALL TEMPERATURE") == {}

tests/test_boundary_global_arrays.py

@@ -94,3 +94,36 @@ def test_vmin_and_vmax_of_a_face_only_some_parts_have(sim):
     coarse = next(sub for sub in table._all_subobstructions if sub.mesh.id == "Coarse")
     assert table.vmax(NODES, 1) == coarse.get_data(NODES).vmax(1)
     assert table.vmin(NODES, 1) == coarse.get_data(NODES).vmin(1)
+
+
+def test_coordinates_of_a_chosen_quantity(sim):
+    # The -x face of 'Table' at x = 0.6 spans y 0.2..0.8: 7 nodes, or 6 cell centres
+    fine = next(sub for sub in sim.obstructions.get_by_id("Table")._all_subobstructions if sub.mesh.id == "Fine")
+    np.testing.assert_allclose(fine.get_coordinates(quantity=NODES)[-1]["y"], np.linspace(0.2, 0.8, 7), atol=1e-6)
+    np.testing.assert_allclose(fine.get_coordinates(quantity=CELLS)[-1]["y"], np.linspace(0.25, 0.75, 6), atol=1e-6)
+    assert fine.get_nearest_index("y", -1, 0.75, quantity=CELLS) == 5
+
+
+def test_boundary_data_of_one_orientation(sim):
+    assert list(sim.obstructions.get_by_id("Table").get_boundary_data(NODES, 1)) == ["Coarse"]
+
+
+def test_parts_without_the_quantity_are_skipped(sim):
+    table = load(sim.root_path).obstructions.get_by_id("Table")  # a fresh copy, as the test removes data
+    fine = next(sub for sub in table._all_subobstructions if sub.mesh.id == "Fine")
+    coarse = next(sub for sub in table._all_subobstructions if sub.mesh.id == "Coarse")
+    bid = next(bid for bid, b in fine._boundary_data.items() if b.quantity == NODES)
+    del fine._boundary_data[bid]
+    assert table.vmax(NODES) == coarse.get_data(NODES).vmax()
+    np.testing.assert_array_equal(table.get_global_boundary_data_arrays(NODES)[1], coarse.get_data(NODES).data[1].data)
+    with pytest.raises(KeyError):
+        table.vmax("VELOCITY")
+
+
+def test_node_based_and_cell_centred_patches_are_not_combined(sim):
+    table = load(sim.root_path).obstructions.get_by_id("Table")  # a fresh copy, as the test changes the patches
+    coarse = next(sub for sub in table._all_subobstructions if sub.mesh.id == "Coarse")
+    for patch in coarse.get_data(NODES).patches:
+        patch.cell_centered = True
+    with pytest.raises(ValueError, match="node-based"):
+        table.get_global_boundary_data_arrays(NODES)