diff --git a/fdsreader/bndf/obstruction.py b/fdsreader/bndf/obstruction.py index 531dd25..242e7a2 100644 --- a/fdsreader/bndf/obstruction.py +++ b/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): diff --git a/fdsreader/fds_classes/mesh.py b/fdsreader/fds_classes/mesh.py index 9972016..c06ed92 100644 --- a/fdsreader/fds_classes/mesh.py +++ b/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 diff --git a/fdsreader/fds_classes/mesh_collection.py b/fdsreader/fds_classes/mesh_collection.py index 91c3435..d26d25c 100644 --- a/fdsreader/fds_classes/mesh_collection.py +++ b/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): diff --git a/tests/cases/bndf_domain_fds611.tgz b/tests/cases/bndf_domain_fds611.tgz index 81e01e5..f9dcd71 100644 --- a/tests/cases/bndf_domain_fds611.tgz +++ b/tests/cases/bndf_domain_fds611.tgz @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:d284eac8e560535b654eff43220812ffb329953e7b58922e45e7d301cd4e2905 -size 10134 +oid sha256:20f1937e92081b7b9ae359a916bd5ec08a8291e9531fd88c416e194de2884adf +size 10097 diff --git a/tests/cases/bndf_global_fds611.tgz b/tests/cases/bndf_global_fds611.tgz index 027951a..dcdcc54 100644 --- a/tests/cases/bndf_global_fds611.tgz +++ b/tests/cases/bndf_global_fds611.tgz @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:6e502ffcf16979705835eb7fe011e7fba41bfae95ddf61390961a1a10fd331d1 -size 25790 +oid sha256:4f0c986d8bfd062b88e51ce04e35d85266359470ae23383ee2fbff88a009ec35 +size 25179 diff --git a/tests/cases/fds_inputs/bndf_domain.fds b/tests/cases/fds_inputs/bndf_domain.fds index d1f7493..08513e5 100644 --- a/tests/cases/fds_inputs/bndf_domain.fds +++ b/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 / diff --git a/tests/cases/fds_inputs/bndf_global.fds b/tests/cases/fds_inputs/bndf_global.fds index 73b9cc7..b61dde0 100644 --- a/tests/cases/fds_inputs/bndf_global.fds +++ b/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 / diff --git a/tests/test_boundary_domain.py b/tests/test_boundary_domain.py index 60c250d..0dd2e71 100644 --- a/tests/test_boundary_domain.py +++ b/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") == {} diff --git a/tests/test_boundary_global_arrays.py b/tests/test_boundary_global_arrays.py index d2b87ee..db21d8f 100644 --- a/tests/test_boundary_global_arrays.py +++ b/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)