diff --git a/.gitignore b/.gitignore
index 9e73f54..289f404 100644
--- a/.gitignore
+++ b/.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
diff --git a/docs/api/bndf.rst b/docs/api/bndf.rst
index af1b1f3..f41e6db 100644
--- a/docs/api/bndf.rst
+++ b/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:
diff --git a/examples/bndf/boundary.ipynb b/examples/bndf/boundary.ipynb
index 50b4764..c7023a4 100644
--- a/examples/bndf/boundary.ipynb
+++ b/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",
diff --git a/fdsreader/bndf/obstruction.py b/fdsreader/bndf/obstruction.py
index cb6ca7c..531dd25 100644
--- a/fdsreader/bndf/obstruction.py
+++ b/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,13 +219,25 @@ 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)
@@ -226,12 +245,12 @@ def vmin(self, orientation: Literal[-3, -2, -1, 0, 1, 2, 3] = 0) -> float:
                 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)
@@ -239,7 +258,7 @@ def vmax(self, orientation: Literal[-3, -2, -1, 0, 1, 2, 3] = 0) -> float:
                 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
diff --git a/fdsreader/export/obst_exporter.py b/fdsreader/export/obst_exporter.py
index 1b1a508..582e68f 100644
--- a/fdsreader/export/obst_exporter.py
+++ b/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]
diff --git a/fdsreader/fds_classes/mesh.py b/fdsreader/fds_classes/mesh.py
index 468e5b3..9972016 100644
--- a/fdsreader/fds_classes/mesh.py
+++ b/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.
diff --git a/fdsreader/fds_classes/mesh_collection.py b/fdsreader/fds_classes/mesh_collection.py
index cc3c40f..91c3435 100644
--- a/fdsreader/fds_classes/mesh_collection.py
+++ b/fdsreader/fds_classes/mesh_collection.py
@@ -1,6 +1,10 @@
-from typing import Iterable
+from typing import Dict, Iterable, Tuple, Union
 
+import numpy as np
+
+from fdsreader.bndf.obstruction import _global_boundary_arrays
 from fdsreader.fds_classes import Mesh
+from fdsreader.utils import Quantity
 from fdsreader.utils.data import FDSDataCollection
 
 
@@ -16,5 +20,20 @@ def get_by_id(self, mesh_id: str):
         """Get the mesh with corresponding id if it exists."""
         return next((mesh for mesh in self if mesh.id == mesh_id), None)
 
+    def get_global_boundary_data_arrays(
+        self, quantity: Union[str, Quantity]
+    ) -> Dict[int, Union[np.ndarray, Tuple[np.ndarray, ...]]]:
+        """Combines the boundary data on the walls of the simulation domain, e.g. the floor, over all meshes into one
+        array per orientation (-3=-z, -2=-y, -1=-x, 1=x, 2=y, 3=z; 3 is the floor, facing +z).
+
+        Faces of obstructions are not included, see :meth:`Obstruction.get_global_boundary_data_arrays`, which also
+        describes the shape of the arrays and when an orientation maps to a tuple of arrays. FDS writes these data
+        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.
+        """
+        return _global_boundary_arrays([mesh.get_boundary_data for mesh in self], quantity)
+
     def __repr__(self):
         return "MeshCollection(" + super().__repr__() + ")"
diff --git a/fdsreader/smoke3d/smoke3d.py b/fdsreader/smoke3d/smoke3d.py
index a6ac4d0..bd86c9e 100644
--- a/fdsreader/smoke3d/smoke3d.py
+++ b/fdsreader/smoke3d/smoke3d.py
@@ -1,6 +1,5 @@
 import math
 import os
-import sys
 import warnings
 from typing import Dict, Literal, Optional, Tuple, Union
 
@@ -11,6 +10,7 @@
 from fdsreader.fds_classes import Mesh
 from fdsreader.utils import NumpyArrayMixin, Quantity
 from fdsreader.utils.data import resample_mesh_axis
+from fdsreader.utils.misc import caller_stacklevel
 
 _HANDLED_FUNCTIONS = {}
 
@@ -21,16 +21,6 @@
 KINDS = ("gas", "hrrpuv", "temperature")
 
 _raw_codes_warned = False
-_PACKAGE_DIR = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
-
-
-def _caller_stacklevel() -> int:
-    """The stacklevel for warnings.warn that points at the first caller outside fdsreader, seen
-    from the function that calls this one."""
-    frame, level = sys._getframe(1), 1
-    while frame is not None and os.path.abspath(frame.f_code.co_filename).startswith(_PACKAGE_DIR + os.sep):
-        frame, level = frame.f_back, level + 1
-    return level
 
 
 def _warn_raw_codes():
@@ -43,7 +33,7 @@ def _warn_raw_codes():
             "to_physical() or to_global(physical=True) to convert them, see "
             "https://fdsreader.readthedocs.io/en/latest/api/smoke3d.html",
             UserWarning,
-            stacklevel=_caller_stacklevel(),
+            stacklevel=caller_stacklevel(),
         )
 
 
@@ -349,7 +339,7 @@ def _limits(self, keyword: str, low: float | None, high: float | None) -> Tuple[
                 f"The .smv file has no {keyword} (written from FDS 6.10 on), assuming {default} for "
                 f"{self.quantity.name}, {source}. Pass the values set on &DUMP if they were changed.",
                 UserWarning,
-                stacklevel=_caller_stacklevel(),
+                stacklevel=caller_stacklevel(),
             )
         base = smv if smv is not None else default
         return (base[0] if low is None else low, base[1] if high is None else high)
diff --git a/fdsreader/utils/misc.py b/fdsreader/utils/misc.py
index 554774c..69a0e17 100644
--- a/fdsreader/utils/misc.py
+++ b/fdsreader/utils/misc.py
@@ -1,8 +1,21 @@
 import logging
+import os
+import sys
 from functools import wraps
 
 from fdsreader import settings
 
+_PACKAGE_DIR = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+
+
+def caller_stacklevel() -> int:
+    """The stacklevel for warnings.warn that points at the first caller outside fdsreader, seen
+    from the function that calls this one."""
+    frame, level = sys._getframe(1), 1
+    while frame is not None and os.path.abspath(frame.f_code.co_filename).startswith(_PACKAGE_DIR + os.sep):
+        frame, level = frame.f_back, level + 1
+    return level
+
 
 def log_error(module):
     def decorated(f):
diff --git a/tests/cases/bndf_domain_fds611.tgz b/tests/cases/bndf_domain_fds611.tgz
new file mode 100644
index 0000000..81e01e5
--- /dev/null
+++ b/tests/cases/bndf_domain_fds611.tgz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d284eac8e560535b654eff43220812ffb329953e7b58922e45e7d301cd4e2905
+size 10134
diff --git a/tests/cases/fds_inputs/bndf_domain.fds b/tests/cases/fds_inputs/bndf_domain.fds
new file mode 100644
index 0000000..d1f7493
--- /dev/null
+++ b/tests/cases/fds_inputs/bndf_domain.fds
@@ -0,0 +1,16 @@
+&HEAD CHID='bndf_domain', TITLE='Boundary data on the walls of the domain for fdsreader (#81)' /
+&MESH ID='West', IJK=10,10,10, XB=0.0,1.0,0.0,1.0,0.0,1.0 /
+&MESH ID='East', IJK=5,5,5, XB=1.0,2.0,0.0,1.0,0.0,1.0 /
+&TIME T_END=2.0 /
+&DUMP NFRAMES=4, SMOKE3D=F /
+&REAC FUEL='PROPANE' /
+&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' /
+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 /
+&BNDF QUANTITY='WALL TEMPERATURE' /
+&TAIL /
diff --git a/tests/test_boundary_domain.py b/tests/test_boundary_domain.py
new file mode 100644
index 0000000..60c250d
--- /dev/null
+++ b/tests/test_boundary_domain.py
@@ -0,0 +1,98 @@
+"""Boundary data on the walls of the simulation domain (#81), on bndf_domain_fds611 (FDS 6.11.1; input in
+tests/cases/fds_inputs/bndf_domain.fds).
+
+Two meshes: 'West' (x 0..1, 0.1 m cells) and 'East' (x 1..2, 0.2 m cells), y 0..1. A burner vent at x 0.6..1.6,
+y 0.3..0.7 on the floor splits the floor of each mesh into four patches; in 'East', FDS snaps the vent to y 0.4..0.8.
+The obstruction 'Box' stands on the floor. BURNING RATE is cell-centred, WALL TEMPERATURE node-based.
+"""
+
+import warnings
+
+import numpy as np
+import pytest
+from fds_cases import case_dirs, case_id, load
+
+CELLS = "BURNING RATE"
+NODES = "WALL TEMPERATURE"
+FLOOR = 3
+
+
+@pytest.fixture(scope="module", params=case_dirs("bndf_domain_fds611", "bndf_domain"), ids=case_id)
+def sim(request):
+    return load(request.param)
+
+
+def _floor_patches(sim, mesh_id, quantity):
+    boundary = sim.meshes.get_by_id(mesh_id).get_boundary_data(quantity)
+    return [p for p in boundary.patches if p.orientation == FLOOR]
+
+
+# The tests also run on the output of new FDS builds, so they do not rely on how FDS splits the floor into patches
+# x = 0 and the top are open, x = 1 lies between the meshes, so 'West' has no faces facing +-x and 'East' only -x
+@pytest.mark.parametrize(
+    "mesh_id, n_cells, orientations", [("West", 10 * 10, [-2, 2, 3]), ("East", 5 * 5, [-2, -1, 2, 3])]
+)
+def test_all_patches_of_a_face_are_kept(sim, mesh_id, n_cells, orientations):
+    patches = _floor_patches(sim, mesh_id, CELLS)
+    assert len(patches) > 1  # the vent splits the floor
+    assert sum(p.data[0].size for p in patches) == n_cells  # together they cover the floor of the mesh
+    boundary = sim.meshes.get_by_id(mesh_id).get_boundary_data(CELLS)
+    assert boundary.orientations == orientations
+    assert boundary.vmax(FLOOR) == max(np.max(p.data) for p in patches)
+
+
+def test_boundary_data_warns_when_it_drops_patches(sim):
+    boundary = sim.meshes.get_by_id("West").get_boundary_data(CELLS)
+    with pytest.warns(UserWarning, match="Boundary.patches") as record:
+        _ = boundary.data
+    warning = next(w for w in record if "Boundary.patches" in str(w.message))
+    assert warning.filename == __file__  # points at the caller, not into fdsreader
+
+
+def test_vmin_and_vmax_of_a_missing_face_are_nan(sim):
+    boundary = sim.meshes.get_by_id("West").get_boundary_data(CELLS)
+    assert np.isnan(boundary.vmin(1)) and np.isnan(boundary.vmax(1))
+
+
+def test_obstruction_boundary_data_does_not_warn(sim):
+    box = sim.obstructions.get_by_id("Box")
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        for sub in box._all_subobstructions:
+            _ = sub.get_data(CELLS).data
+
+
+def test_the_floor_over_both_meshes(sim):
+    floor = sim.meshes.get_global_boundary_data_arrays(CELLS)[FLOOR]
+    assert floor.shape == (sim.meshes[0].get_boundary_data(CELLS).n_t(count_duplicates=False), 20, 10)
+    assert not np.isnan(floor).any()
+
+    burning = np.zeros((20, 10), dtype=bool)
+    burning[6:10, 3:7] = True  # 'West': x 0.6..1.0, y 0.3..0.7
+    burning[10:16, 4:8] = True  # 'East': x 1.0..1.6, y 0.4..0.8
+    np.testing.assert_array_equal(floor[-1] > 0, burning)
+
+
+@pytest.mark.parametrize("mesh_id, repeat", [("West", 1), ("East", 2)])  # 'East' has cells twice as large
+def test_the_floor_holds_each_patch(sim, mesh_id, repeat):
+    floor = sim.meshes.get_global_boundary_data_arrays(CELLS)[FLOOR]
+    for patch in _floor_patches(sim, mesh_id, CELLS):
+        x = slice(round(patch.extent.x_start / 0.1), round(patch.extent.x_end / 0.1))
+        y = slice(round(patch.extent.y_start / 0.1), round(patch.extent.y_end / 0.1))
+        np.testing.assert_array_equal(floor[:, x, y], patch.data.repeat(repeat, axis=1).repeat(repeat, axis=2))
+
+
+def test_node_based_floor(sim):
+    floor = sim.meshes.get_global_boundary_data_arrays(NODES)[FLOOR]
+    assert floor.shape[1:] == (21, 11)
+    assert not np.isnan(floor).any()
+
+
+def test_obstruction_faces_are_not_included(sim):
+    # The top of 'Box' at z = 0.2 also faces +z, but belongs to the obstruction
+    assert isinstance(sim.meshes.get_global_boundary_data_arrays(CELLS)[FLOOR], np.ndarray)
+
+
+def test_unknown_quantity(sim):
+    with pytest.raises(KeyError):
+        sim.meshes.get_global_boundary_data_arrays("VELOCITY")
diff --git a/tests/test_boundary_global_arrays.py b/tests/test_boundary_global_arrays.py
index 9c102a2..d2b87ee 100644
--- a/tests/test_boundary_global_arrays.py
+++ b/tests/test_boundary_global_arrays.py
@@ -86,3 +86,11 @@ def test_export_of_a_face_at_two_heights(sim, tmp_path):
     meta = export_obst_raw(sim.obstructions.get_by_id("Table"), str(tmp_path))
     assert (tmp_path / "WALL_TEMPERATURE" / "obst-Table_quantity-WALL_TEMPERATURE.dat").stat().st_size > 0
     assert meta.endswith("obst-Table.yaml")
+
+
+def test_vmin_and_vmax_of_a_face_only_some_parts_have(sim):
+    # Only the part of 'Table' in 'Coarse' has a +x face, the part in 'Fine' ends at the mesh boundary
+    table = sim.obstructions.get_by_id("Table")
+    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)