feat: boundary data on the walls of the simulation domain

.gitignore

@@ -535,6 +535,7 @@ tests/cases/cface_data_fds*/
 tests/cases/slcf_partial_fds*/
 tests/cases/slcf_masks_fds*/
 tests/cases/bndf_global_fds*/
+tests/cases/bndf_domain_fds*/
 
 # Pickle cache files
 *.pickle

docs/api/bndf.rst

@@ -1,6 +1,12 @@
 Boundary data
 =============
 
+Boundary files (``&BNDF``) hold a quantity on solid surfaces. Data on the faces of an obstruction belong to the
+:class:`~fdsreader.bndf.obstruction.Obstruction`, see ``Obstruction.get_global_boundary_data_arrays``. Data on the
+walls of the simulation domain that are not part of an obstruction, e.g. the ground of an outdoor case, belong to
+the meshes: ``Mesh.get_boundary_data`` for one mesh, and ``sim.meshes.get_global_boundary_data_arrays`` for all
+meshes combined.
+
 .. automodule:: fdsreader.bndf.obstruction_collection
    :autosummary:
    :members:

examples/bndf/boundary.ipynb

@@ -190,6 +190,28 @@
     "ax.set(xlabel=\"Time [s]\", ylabel=f\"max. {quantity} [C]\", title=f\"{wall.id}, +x face\");"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "bf244d54",
+   "metadata": {},
+   "source": [
+    "## Walls of the simulation domain\n",
+    "\n",
+    "Data on walls of the simulation domain that are not part of an obstruction, for example the\n",
+    "burning rate on the ground of an outdoor case, belong to the meshes. FDS may split such a wall into\n",
+    "several patches, e.g. where a vent begins; `sim.meshes.get_global_boundary_data_arrays` combines\n",
+    "them over all meshes in the same way as for an obstruction:\n",
+    "\n",
+    "```python\n",
+    "walls = sim.meshes.get_global_boundary_data_arrays(\"BURNING RATE\")\n",
+    "ground = walls[3]  # faces +z: (time, x, y)\n",
+    "```\n",
+    "\n",
+    "In this example these data only hold the ambient temperature: they lie below the obstruction\n",
+    "`Boden` or on faces between meshes. The input sets its open boundaries with `MB`, which applies\n",
+    "to the faces of every mesh; `DB` would apply them to the domain only (see the FDS User Guide)."
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "6013d7a5",

fdsreader/bndf/obstruction.py

@@ -1,15 +1,16 @@
 import math
 import os
+import warnings
 
 # Unfortunately, this is necessary due to a cyclic reference. "Mesh" is only needed for static type hints anyway
-from typing import TYPE_CHECKING, Dict, List, Sequence, Tuple, Union
+from typing import TYPE_CHECKING, Callable, Dict, List, Sequence, Tuple, Union
 
 import numpy as np
 from typing_extensions import Literal
 
 import fdsreader.utils.fortran_data as fdtype
-from fdsreader import settings
 from fdsreader.utils import Dimension, Extent, Quantity
+from fdsreader.utils.misc import caller_stacklevel
 
 if TYPE_CHECKING:
     from fdsreader.fds_classes import Mesh
@@ -195,10 +196,16 @@ def n_t(self, count_duplicates=True) -> int:
         """
         return self._patches[0].n_t(count_duplicates=count_duplicates)
 
+    @property
+    def patches(self) -> List[Patch]:
+        """All patches, in the order FDS wrote them. A face of the simulation domain may consist of several patches
+        with the same orientation, e.g. where a vent begins."""
+        return list(self._patches)
+
     @property
     def orientations(self):
         """Return all orientations for which there is data available."""
-        return [p.orientation for p in self._patches]
+        return sorted({p.orientation for p in self._patches})
 
     def get_nearest_timestep(self, time: float) -> int:
         """Calculates the nearest timestep for which data has been output for this obstruction."""
@@ -212,34 +219,46 @@ def get_nearest_timestep(self, time: float) -> int:
 
     @property
     def data(self) -> Dict[int, Patch]:
-        """The :class:`Patch` in each direction (-3=-z, -2=-y, -1=-x, 1=x, 2=y, 3=y)."""
-        return {p.orientation: p for p in self._patches}
+        """The :class:`Patch` in each direction (-3=-z, -2=-y, -1=-x, 1=x, 2=y, 3=z).
+
+        A face of the simulation domain may consist of several patches with the same orientation; this only holds
+        the last of them and warns. Use :attr:`patches` or ``get_global_boundary_data_arrays`` for those faces.
+        """
+        data = {p.orientation: p for p in self._patches}
+        if len(data) < len(self._patches):
+            warnings.warn(
+                f"{self.quantity.name}: several patches face the same direction, Boundary.data only holds the last "
+                "of each. Use Boundary.patches, or get_global_boundary_data_arrays to combine them.",
+                UserWarning,
+                stacklevel=caller_stacklevel(),
+            )
+        return data
 
     def vmin(self, orientation: Literal[-3, -2, -1, 0, 1, 2, 3] = 0) -> float:
         """Minimum value of all patches at any time.
 
-        :param orientation: Optionally filter by patches with a specific orientation.
+        :param orientation: Optionally filter by patches with a specific orientation; NaN if there is none.
         """
         if orientation == 0:
             curr_min = np.min(self.lower_bounds)
             if curr_min == np.float32(+1e33):
                 return float(min(np.min(p.data) for p in self._patches))
             return float(curr_min)
         else:
-            return float(np.min(self.data[orientation].data))
+            return float(min((np.min(p.data) for p in self._patches if p.orientation == orientation), default=np.nan))
 
     def vmax(self, orientation: Literal[-3, -2, -1, 0, 1, 2, 3] = 0) -> float:
         """Maximum value of all patches at any time.
 
-        :param orientation: Optionally filter by patches with a specific orientation.
+        :param orientation: Optionally filter by patches with a specific orientation; NaN if there is none.
         """
         if orientation == 0:
             curr_max = np.max(self.upper_bounds)
             if curr_max == np.float32(-1e33):
                 return float(max(np.max(p.data) for p in self._patches))
             return float(curr_max)
         else:
-            return float(np.max(self.data[orientation].data))
+            return float(max((np.max(p.data) for p in self._patches if p.orientation == orientation), default=np.nan))
 
     def clear_cache(self):
         """Remove all data from the internal cache that has been loaded so far to free memory."""
@@ -297,9 +316,6 @@ def _add_patches(
             for patch in patches:
                 patch._boundary_parent = self._boundary_data[bid]
 
-        if not settings.LAZY_LOAD:
-            _ = self._boundary_data[bid].data
-
     @property
     def orientations(self):
         """Return all orientations for which there is data available."""
@@ -317,8 +333,8 @@ def get_coordinates(
         """
         if self.has_boundary_data:
             return {
-                orientation: patch.get_coordinates(ignore_cell_centered)
-                for orientation, patch in next(iter(self._boundary_data.values())).data.items()
+                patch.orientation: patch.get_coordinates(ignore_cell_centered)
+                for patch in next(iter(self._boundary_data.values())).patches
             }
         return {}
 
@@ -338,19 +354,7 @@ def has_boundary_data(self):
         return len(self._boundary_data) != 0
 
     def get_data(self, quantity: Union[str, Quantity]):
-        if isinstance(quantity, Quantity):
-            quantity = quantity.name
-        b = next(
-            (
-                b
-                for b in self._boundary_data.values()
-                if b.quantity.name.lower() == quantity.lower() or b.quantity.short_name.lower() == quantity.lower()
-            ),
-            None,
-        )
-        if b is None:
-            raise KeyError(f"No boundary data found for quantity '{quantity}'")
-        return b
+        return _find_boundary(self._boundary_data, quantity)
 
     def __getitem__(self, item):
         if isinstance(item, int):
@@ -641,7 +645,7 @@ def get_boundary_data(
         }
         if orientation == 0:
             return ret
-        return {mesh: bndf for mesh, bndf in ret.items() if orientation in bndf.data.keys()}
+        return {mesh: bndf for mesh, bndf in ret.items() if orientation in bndf.orientations}
 
     def get_nearest_timestep(self, time: float, visible_only: bool = False) -> int:
         """Calculates the nearest timestep for which data has been output for this obstruction."""
@@ -712,33 +716,9 @@ def get_global_boundary_data_arrays(
 
         :param quantity: The quantity's name or short name for which to load the global arrays.
         """
-        patches_by_orientation: Dict[int, List[Patch]] = dict()
-        for subobst in self._all_subobstructions:
-            try:
-                boundary = subobst.get_data(quantity)
-            except KeyError:
-                continue
-            for orientation, patch in boundary.data.items():
-                patches_by_orientation.setdefault(orientation, list()).append(patch)
-        if not patches_by_orientation and self.has_boundary_data:
-            raise KeyError(f"No boundary data found for quantity '{quantity}'")
-
-        result = dict()
-        for orientation in sorted(patches_by_orientation):
-            normal = ("x", "y", "z")[abs(orientation) - 1]
-            levels: List[Tuple[float, List[Patch]]] = list()
-            for patch in patches_by_orientation[orientation]:
-                position = patch.extent[normal][0]
-                tol = _tolerance(patch.mesh.coordinates[normal])
-                level = next((lvl for lvl in levels if abs(lvl[0] - position) <= tol), None)
-                if level is None:
-                    levels.append((position, [patch]))
-                else:
-                    level[1].append(patch)
-            levels.sort(key=lambda lvl: lvl[0])
-            grids = tuple(_global_patch_grid(level_patches, normal) for _, level_patches in levels)
-            result[orientation] = grids[0] if len(grids) == 1 else grids
-        return result
+        if not self.has_boundary_data:
+            return dict()
+        return _global_boundary_arrays([subobst.get_data for subobst in self._all_subobstructions], quantity)
 
     def clear_cache(self):
         """Remove all data from the internal cache that has been loaded so far to free memory."""
@@ -754,10 +734,8 @@ def vmin(
 
         :param orientation: Optionally filter by patches with a specific orientation.
         """
-        if self.has_boundary_data:
-            return min(s.vmin(quantity, orientation) for s in self._all_subobstructions)
-        else:
-            return np.nan
+        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)
 
     def vmax(
         self,
@@ -768,9 +746,8 @@ def vmax(
 
         :param orientation: Optionally filter by patches with a specific orientation.
         """
-        if self.has_boundary_data:
-            return max(s.vmax(quantity, orientation) for s in self._all_subobstructions)
-        return np.nan
+        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)
 
     def get_subobstructions(self, mesh: "Mesh") -> List[SubObstruction]:
         """All parts of this obstruction in the given mesh; several OBST lines can share an ID."""
@@ -802,6 +779,56 @@ def __repr__(self, *args, **kwargs):
         )
 
 
+def _find_boundary(boundary_data: Dict[int, Boundary], quantity: Union[str, Quantity]) -> Boundary:
+    """The boundary data of a quantity, by its name or short name, among those of a mesh or subobstruction."""
+    if isinstance(quantity, Quantity):
+        quantity = quantity.name
+    # Quantity == str matches the name or the short name, ignoring case
+    b = next((b for b in boundary_data.values() if b.quantity == quantity), None)
+    if b is None:
+        raise KeyError(f"No boundary data found for quantity '{quantity}'")
+    return b
+
+
+def _global_boundary_arrays(
+    lookups: List[Callable[[Union[str, Quantity]], Boundary]], quantity: Union[str, Quantity]
+) -> Dict[int, Union[np.ndarray, Tuple[np.ndarray, ...]]]:
+    """Combines the patches of a quantity from several meshes or subobstructions into one array of shape (n_t, a, b)
+    per orientation, or a tuple of them, one per position, if the faces of an orientation lie at several positions.
+    See :meth:`Obstruction.get_global_boundary_data_arrays`.
+
+    :param lookups: The boundary data lookup of each mesh or subobstruction, which raises KeyError if it does not have
+        the quantity. Raises KeyError if none of them has it.
+    """
+    patches_by_orientation: Dict[int, List[Patch]] = dict()
+    for lookup in lookups:
+        try:
+            boundary = lookup(quantity)
+        except KeyError:
+            continue
+        for patch in boundary.patches:
+            patches_by_orientation.setdefault(patch.orientation, list()).append(patch)
+    if not patches_by_orientation:
+        raise KeyError(f"No boundary data found for quantity '{quantity}'")
+
+    result = dict()
+    for orientation in sorted(patches_by_orientation):
+        normal = ("x", "y", "z")[abs(orientation) - 1]
+        levels: List[Tuple[float, List[Patch]]] = list()
+        for patch in patches_by_orientation[orientation]:
+            position = patch.extent[normal][0]
+            tol = _tolerance(patch.mesh.coordinates[normal])
+            level = next((lvl for lvl in levels if abs(lvl[0] - position) <= tol), None)
+            if level is None:
+                levels.append((position, [patch]))
+            else:
+                level[1].append(patch)
+        levels.sort(key=lambda lvl: lvl[0])
+        grids = tuple(_global_patch_grid(level_patches, normal) for _, level_patches in levels)
+        result[orientation] = grids[0] if len(grids) == 1 else grids
+    return result
+
+
 def _tolerance(coordinates: np.ndarray) -> float:
     """Tolerance for comparing positions on a grid: a thousandth of its smallest cell."""
     return 1e-3 * float(np.min(np.diff(coordinates))) if coordinates.size > 1 else 1e-6

fdsreader/export/obst_exporter.py

@@ -35,8 +35,11 @@ def export_obst_raw(obst: Obstruction, output_dir: str, ordering: Literal["C", "
 
     random_bndf = next(iter(obst.get_boundary_data(obst.quantities[0]).values()))
     meta["TimeSteps"] = len(random_bndf.times)
-    meta["NumOrientations"] = len(random_bndf.data)
-    for orientation, face in random_bndf.data.items():
+    faces = random_bndf.data
+    meta["NumOrientations"] = len(faces)
+    # Sorted like get_global_boundary_data_arrays. These are the faces of one part of the obstruction; for an
+    # obstruction spanning several meshes, the data files may hold faces that are not listed here.
+    for orientation, face in sorted(faces.items()):
         if abs(orientation) == 1:
             spacing1 = (bounding_box[3] - bounding_box[2]) / face.shape[0]
             spacing2 = (bounding_box[5] - bounding_box[4]) / face.shape[1]

fdsreader/fds_classes/mesh.py

@@ -4,8 +4,8 @@
 import numpy as np
 from typing_extensions import Literal
 
-from fdsreader import settings
 from fdsreader.bndf import Boundary, Patch
+from fdsreader.bndf.obstruction import _find_boundary
 from fdsreader.utils import Dimension, Extent, Quantity
 
 
@@ -228,23 +228,17 @@ def _add_patches(
             for patch in patches:
                 patch._boundary_parent = self._boundary_data[bid]
 
-        if not settings.LAZY_LOAD:
-            _ = self._boundary_data[bid].data
-
-    def get_boundary_data(self, quantity: Union[str, Quantity]):
-        if isinstance(quantity, Quantity):
-            quantity = quantity.name
-        b = next(
-            (
-                b
-                for b in self._boundary_data.values()
-                if b.quantity.name.lower() == quantity.lower() or b.quantity.short_name.lower() == quantity.lower()
-            ),
-            None,
-        )
-        if b is None:
-            raise KeyError(f"No boundary data found for quantity '{quantity}'")
-        return b
+    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
+        the FDS User Guide), FDS may also write faces between two meshes.
+
+        A face may consist of several patches with the same orientation (see :attr:`Boundary.patches`); to combine
+        them over all meshes, use :meth:`MeshCollection.get_global_boundary_data_arrays`.
+
+        :param quantity: The quantity's name or short name.
+        """
+        return _find_boundary(self._boundary_data, quantity)
 
     def __getitem__(self, dimension: Literal[0, 1, 2, "x", "y", "z"]) -> np.ndarray:
         """Get all values in given dimension.