From 4896833ed3c4eaf3ef44a749d016ea2a6902530d Mon Sep 17 00:00:00 2001
From: Marc Fehr <marc.fehr@hhu.de>
Date: Thu, 15 Oct 2026 22:27:57 +0200
Subject: [PATCH 1/3] feat: boundary data on the walls of the simulation domain

Boundary data on walls of the domain that are not part of an obstruction, e.g.
the ground of an outdoor case, are attached to the meshes. FDS may split such a
wall into several patches with the same orientation, e.g. where a vent begins,
but Boundary.data maps each orientation to one patch, so all but the last were
dropped; vmin/vmax per orientation saw only that patch, and orientations listed
duplicates.

- Boundary.patches lists all patches; orientations are unique and sorted;
  vmin/vmax cover all patches of an orientation. Boundary.data stays as it is
  but warns when it drops patches.
- MeshCollection.get_global_boundary_data_arrays combines the domain walls over
  all meshes, with the same grid code as the obstruction faces.
- Mesh.get_boundary_data has a docstring, the API page and the boundary
  notebook say where these data live.

The new case bndf_domain (FDS 6.11.1) has two meshes with 0.1 m and 0.2 m cells
and a burner vent splitting the floor of each into four patches.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .gitignore                               |  1 +
 docs/api/bndf.rst                        |  6 ++
 examples/bndf/boundary.ipynb             | 21 ++++++
 fdsreader/bndf/obstruction.py            | 84 ++++++++++++++++--------
 fdsreader/fds_classes/mesh.py            | 13 +++-
 fdsreader/fds_classes/mesh_collection.py | 29 +++++++-
 tests/cases/bndf_domain_fds611.tgz       |  3 +
 tests/cases/fds_inputs/bndf_domain.fds   | 16 +++++
 tests/test_boundary_domain.py            | 81 +++++++++++++++++++++++
 9 files changed, 223 insertions(+), 31 deletions(-)
 create mode 100644 tests/cases/bndf_domain_fds611.tgz
 create mode 100644 tests/cases/fds_inputs/bndf_domain.fds
 create mode 100644 tests/test_boundary_domain.py

diff --git a/.gitignore b/.gitignore
index 9e73f54a..289f404e 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 af1b1f3f..f41e6dbd 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 50b47642..04b27677 100644
--- a/examples/bndf/boundary.ipynb
+++ b/examples/bndf/boundary.ipynb
@@ -190,6 +190,27 @@
     "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 the room stands on the obstruction `Boden` and the outer faces are open, so\n",
+    "these data only hold the ambient temperature."
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "6013d7a5",
diff --git a/fdsreader/bndf/obstruction.py b/fdsreader/bndf/obstruction.py
index cb6ca7cc..479a31b2 100644
--- a/fdsreader/bndf/obstruction.py
+++ b/fdsreader/bndf/obstruction.py
@@ -1,5 +1,6 @@
 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
@@ -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,8 +219,20 @@ 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=2,
+            )
+        return data
 
     def vmin(self, orientation: Literal[-3, -2, -1, 0, 1, 2, 3] = 0) -> float:
         """Minimum value of all patches at any time.
@@ -226,7 +245,7 @@ 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))
 
     def vmax(self, orientation: Literal[-3, -2, -1, 0, 1, 2, 3] = 0) -> float:
         """Maximum value of all patches at any time.
@@ -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))
 
     def clear_cache(self):
         """Remove all data from the internal cache that has been loaded so far to free memory."""
@@ -298,7 +317,8 @@ def _add_patches(
                 patch._boundary_parent = self._boundary_data[bid]
 
         if not settings.LAZY_LOAD:
-            _ = self._boundary_data[bid].data
+            for patch in patches:
+                _ = patch.data
 
     @property
     def orientations(self):
@@ -712,33 +732,15 @@ 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()
+        patches = list()
         for subobst in self._all_subobstructions:
             try:
-                boundary = subobst.get_data(quantity)
+                patches.extend(subobst.get_data(quantity).patches)
             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:
+        if not patches 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
+        return _global_boundary_arrays(patches)
 
     def clear_cache(self):
         """Remove all data from the internal cache that has been loaded so far to free memory."""
@@ -802,6 +804,32 @@ def __repr__(self, *args, **kwargs):
         )
 
 
+def _global_boundary_arrays(patches: List[Patch]) -> Dict[int, Union[np.ndarray, Tuple[np.ndarray, ...]]]:
+    """Combines patches 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`.
+    """
+    patches_by_orientation: Dict[int, List[Patch]] = dict()
+    for patch in patches:
+        patches_by_orientation.setdefault(patch.orientation, list()).append(patch)
+
+    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/fds_classes/mesh.py b/fdsreader/fds_classes/mesh.py
index 468e5b37..2808485f 100644
--- a/fdsreader/fds_classes/mesh.py
+++ b/fdsreader/fds_classes/mesh.py
@@ -229,9 +229,18 @@ def _add_patches(
                 patch._boundary_parent = self._boundary_data[bid]
 
         if not settings.LAZY_LOAD:
-            _ = self._boundary_data[bid].data
+            for patch in patches:
+                _ = patch.data
+
+    def get_boundary_data(self, quantity: Union[str, Quantity]) -> Boundary:
+        """The boundary data of a quantity on the faces of this mesh that are walls of the simulation domain, e.g.
+        the floor, rather than faces of obstructions (those belong to :class:`Obstruction`).
 
-    def get_boundary_data(self, quantity: Union[str, Quantity]):
+        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.
+        """
         if isinstance(quantity, Quantity):
             quantity = quantity.name
         b = next(
diff --git a/fdsreader/fds_classes/mesh_collection.py b/fdsreader/fds_classes/mesh_collection.py
index cc3c40f6..32a2983e 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,28 @@ 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.
+        """
+        patches = list()
+        for mesh in self:
+            try:
+                patches.extend(mesh.get_boundary_data(quantity).patches)
+            except KeyError:
+                continue
+        if not patches:
+            raise KeyError(f"No boundary data on the walls of the simulation domain for quantity '{quantity}'")
+        return _global_boundary_arrays(patches)
+
     def __repr__(self):
         return "MeshCollection(" + super().__repr__() + ")"
diff --git a/tests/cases/bndf_domain_fds611.tgz b/tests/cases/bndf_domain_fds611.tgz
new file mode 100644
index 00000000..81e01e57
--- /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 00000000..d1f74937
--- /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 00000000..4d2f99ad
--- /dev/null
+++ b/tests/test_boundary_domain.py
@@ -0,0 +1,81 @@
+"""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 sorted(
+        (p for p in boundary.patches if p.orientation == FLOOR), key=lambda p: (p.extent.x_start, p.extent.y_start)
+    )
+
+
+def test_all_patches_of_a_face_are_kept(sim):
+    boundary = sim.meshes.get_by_id("West").get_boundary_data(CELLS)
+    assert len([p for p in boundary.patches if p.orientation == FLOOR]) == 4
+    assert boundary.orientations == sorted(set(boundary.orientations))
+    assert boundary.vmax(FLOOR) == max(np.max(p.data) for p in _floor_patches(sim, "West", CELLS))
+    with pytest.warns(UserWarning, match="Boundary.patches"):
+        _ = boundary.data
+
+
+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)
+
+
+def test_the_floor_holds_each_patch(sim):
+    floor = sim.meshes.get_global_boundary_data_arrays(CELLS)[FLOOR]
+    west = _floor_patches(sim, "West", CELLS)[0]  # x 0..0.6, y 0..1
+    np.testing.assert_array_equal(floor[:, 0:6, 0:10], west.data)
+    east = _floor_patches(sim, "East", CELLS)[0]  # x 1..2, y 0..0.4, 0.2 m cells
+    np.testing.assert_array_equal(floor[:, 10:20, 0:4], east.data.repeat(2, axis=1).repeat(2, 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")

From 9710e4fbb84d238a7a073664668cfe95238e80a6 Mon Sep 17 00:00:00 2001
From: Marc Fehr <marc.fehr@hhu.de>
Date: Thu, 15 Oct 2026 22:35:31 +0200
Subject: [PATCH 2/3] refactor: address the review of the domain boundary data

- Boundary.vmin/vmax for an orientation without patches return NaN, and
  Obstruction.vmin/vmax skip parts without that face instead of failing; an
  obstruction across meshes has faces only some of its parts write.
- export_obst_raw writes the orientation metadata in the sorted order the data
  follow.
- Mesh and SubObstruction share one quantity lookup, which the global arrays of
  meshes and obstructions use as well.
- The warning of Boundary.data points at the first caller outside fdsreader,
  with the helper moved from smoke3d to utils.misc.
- The Mesh.get_boundary_data docstring mentions faces between meshes.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 fdsreader/bndf/obstruction.py            | 85 +++++++++++++-----------
 fdsreader/export/obst_exporter.py        |  6 +-
 fdsreader/fds_classes/mesh.py            | 20 ++----
 fdsreader/fds_classes/mesh_collection.py | 10 +--
 fdsreader/smoke3d/smoke3d.py             | 16 +----
 fdsreader/utils/misc.py                  | 13 ++++
 tests/test_boundary_domain.py            | 10 ++-
 tests/test_boundary_global_arrays.py     |  8 +++
 8 files changed, 88 insertions(+), 80 deletions(-)

diff --git a/fdsreader/bndf/obstruction.py b/fdsreader/bndf/obstruction.py
index 479a31b2..c6a9f070 100644
--- a/fdsreader/bndf/obstruction.py
+++ b/fdsreader/bndf/obstruction.py
@@ -11,6 +11,7 @@
 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
@@ -230,14 +231,14 @@ def data(self) -> Dict[int, Patch]:
                 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=2,
+                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)
@@ -245,12 +246,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(min(np.min(p.data) for p in self._patches if p.orientation == orientation))
+            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)
@@ -258,7 +259,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(max(np.max(p.data) for p in self._patches if p.orientation == orientation))
+            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."""
@@ -358,19 +359,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):
@@ -732,15 +721,9 @@ def get_global_boundary_data_arrays(
 
         :param quantity: The quantity's name or short name for which to load the global arrays.
         """
-        patches = list()
-        for subobst in self._all_subobstructions:
-            try:
-                patches.extend(subobst.get_data(quantity).patches)
-            except KeyError:
-                continue
-        if not patches and self.has_boundary_data:
-            raise KeyError(f"No boundary data found for quantity '{quantity}'")
-        return _global_boundary_arrays(patches)
+        if not self.has_boundary_data:
+            return dict()
+        return _global_boundary_arrays([subobst._boundary_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."""
@@ -756,10 +739,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 float(np.nanmin(values)) if not np.isnan(values).all() else np.nan
 
     def vmax(
         self,
@@ -770,9 +751,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 float(np.nanmax(values)) if not np.isnan(values).all() else 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."""
@@ -804,13 +784,40 @@ def __repr__(self, *args, **kwargs):
         )
 
 
-def _global_boundary_arrays(patches: List[Patch]) -> Dict[int, Union[np.ndarray, Tuple[np.ndarray, ...]]]:
-    """Combines patches 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`.
+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
+    b = next(
+        (
+            b
+            for b in 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 _global_boundary_arrays(
+    owners_boundary_data: List[Dict[int, 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`. Raises KeyError if none of them has the quantity.
     """
     patches_by_orientation: Dict[int, List[Patch]] = dict()
-    for patch in patches:
-        patches_by_orientation.setdefault(patch.orientation, list()).append(patch)
+    for boundary_data in owners_boundary_data:
+        try:
+            boundary = _find_boundary(boundary_data, 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):
diff --git a/fdsreader/export/obst_exporter.py b/fdsreader/export/obst_exporter.py
index 1b1a5084..c894354f 100644
--- a/fdsreader/export/obst_exporter.py
+++ b/fdsreader/export/obst_exporter.py
@@ -35,8 +35,10 @@ 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)
+    # In the order of get_global_boundary_data_arrays, which the data files follow
+    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 2808485f..0ac02359 100644
--- a/fdsreader/fds_classes/mesh.py
+++ b/fdsreader/fds_classes/mesh.py
@@ -6,6 +6,7 @@
 
 from fdsreader import settings
 from fdsreader.bndf import Boundary, Patch
+from fdsreader.bndf.obstruction import _find_boundary
 from fdsreader.utils import Dimension, Extent, Quantity
 
 
@@ -233,27 +234,16 @@ def _add_patches(
                 _ = patch.data
 
     def get_boundary_data(self, quantity: Union[str, Quantity]) -> Boundary:
-        """The boundary data of a quantity on the faces of this mesh that are walls of the simulation domain, e.g.
-        the floor, rather than faces of obstructions (those belong to :class:`Obstruction`).
+        """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.
         """
-        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, 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 32a2983e..4b47a099 100644
--- a/fdsreader/fds_classes/mesh_collection.py
+++ b/fdsreader/fds_classes/mesh_collection.py
@@ -33,15 +33,7 @@ def get_global_boundary_data_arrays(
 
         :param quantity: The quantity's name or short name.
         """
-        patches = list()
-        for mesh in self:
-            try:
-                patches.extend(mesh.get_boundary_data(quantity).patches)
-            except KeyError:
-                continue
-        if not patches:
-            raise KeyError(f"No boundary data on the walls of the simulation domain for quantity '{quantity}'")
-        return _global_boundary_arrays(patches)
+        return _global_boundary_arrays([mesh._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 a6ac4d0d..bd86c9e1 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 554774c6..69a0e17c 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/test_boundary_domain.py b/tests/test_boundary_domain.py
index 4d2f99ad..c1872308 100644
--- a/tests/test_boundary_domain.py
+++ b/tests/test_boundary_domain.py
@@ -32,10 +32,16 @@ def _floor_patches(sim, mesh_id, quantity):
 def test_all_patches_of_a_face_are_kept(sim):
     boundary = sim.meshes.get_by_id("West").get_boundary_data(CELLS)
     assert len([p for p in boundary.patches if p.orientation == FLOOR]) == 4
-    assert boundary.orientations == sorted(set(boundary.orientations))
+    assert boundary.orientations == [-2, 2, 3]  # x = 0 is open, x = 1 lies between the meshes, the top is open
     assert boundary.vmax(FLOOR) == max(np.max(p.data) for p in _floor_patches(sim, "West", CELLS))
-    with pytest.warns(UserWarning, match="Boundary.patches"):
+    with pytest.warns(UserWarning, match="Boundary.patches") as record:
         _ = boundary.data
+    assert record[0].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):
diff --git a/tests/test_boundary_global_arrays.py b/tests/test_boundary_global_arrays.py
index 9c102a2b..d2b87eee 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)

From 998ac01f365cc67f3fd29377202a9a787d07fa66 Mon Sep 17 00:00:00 2001
From: Marc Fehr <marc.fehr@hhu.de>
Date: Thu, 15 Oct 2026 22:44:41 +0200
Subject: [PATCH 3/3] refactor: address the second review of the domain
 boundary data

- LAZY_LOAD=False keeps its previous effect on boundary data (none): reading
  every patch at load time would have been a silent change.
- get_boundary_data(quantity, orientation) and SubObstruction.get_coordinates
  no longer go through Boundary.data, which may warn.
- The global arrays take the public lookups of meshes and subobstructions
  instead of their internal dicts; the quantity lookup uses Quantity.__eq__.
- Obstruction.vmin/vmax reduce without the NaN values in one expression.
- The domain tests no longer depend on how FDS splits the floor into patches,
  since they also run on the output of new FDS builds; the warning test picks
  the matching warning.
- The notebook states what the ambient-only data in the example are.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 examples/bndf/boundary.ipynb             |  5 +--
 fdsreader/bndf/obstruction.py            | 40 ++++++++++--------------
 fdsreader/export/obst_exporter.py        |  3 +-
 fdsreader/fds_classes/mesh.py            |  5 ---
 fdsreader/fds_classes/mesh_collection.py |  2 +-
 tests/test_boundary_domain.py            | 37 ++++++++++++++--------
 6 files changed, 46 insertions(+), 46 deletions(-)

diff --git a/examples/bndf/boundary.ipynb b/examples/bndf/boundary.ipynb
index 04b27677..c7023a4f 100644
--- a/examples/bndf/boundary.ipynb
+++ b/examples/bndf/boundary.ipynb
@@ -207,8 +207,9 @@
     "ground = walls[3]  # faces +z: (time, x, y)\n",
     "```\n",
     "\n",
-    "In this example the room stands on the obstruction `Boden` and the outer faces are open, so\n",
-    "these data only hold the ambient temperature."
+    "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)."
    ]
   },
   {
diff --git a/fdsreader/bndf/obstruction.py b/fdsreader/bndf/obstruction.py
index c6a9f070..531dd25d 100644
--- a/fdsreader/bndf/obstruction.py
+++ b/fdsreader/bndf/obstruction.py
@@ -3,13 +3,12 @@
 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
 
@@ -317,10 +316,6 @@ def _add_patches(
             for patch in patches:
                 patch._boundary_parent = self._boundary_data[bid]
 
-        if not settings.LAZY_LOAD:
-            for patch in patches:
-                _ = patch.data
-
     @property
     def orientations(self):
         """Return all orientations for which there is data available."""
@@ -338,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 {}
 
@@ -650,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."""
@@ -723,7 +718,7 @@ def get_global_boundary_data_arrays(
         """
         if not self.has_boundary_data:
             return dict()
-        return _global_boundary_arrays([subobst._boundary_data for subobst in self._all_subobstructions], quantity)
+        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."""
@@ -740,7 +735,7 @@ 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 float(np.nanmin(values)) if not np.isnan(values).all() else np.nan
+        return min((v for v in values if not np.isnan(v)), default=np.nan)
 
     def vmax(
         self,
@@ -752,7 +747,7 @@ 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 float(np.nanmax(values)) if not np.isnan(values).all() else np.nan
+        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."""
@@ -788,30 +783,27 @@ def _find_boundary(boundary_data: Dict[int, Boundary], quantity: Union[str, Quan
     """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
-    b = next(
-        (
-            b
-            for b in boundary_data.values()
-            if b.quantity.name.lower() == quantity.lower() or b.quantity.short_name.lower() == quantity.lower()
-        ),
-        None,
-    )
+    # 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(
-    owners_boundary_data: List[Dict[int, Boundary]], quantity: Union[str, Quantity]
+    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`. Raises KeyError if none of them has the quantity.
+    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 boundary_data in owners_boundary_data:
+    for lookup in lookups:
         try:
-            boundary = _find_boundary(boundary_data, quantity)
+            boundary = lookup(quantity)
         except KeyError:
             continue
         for patch in boundary.patches:
diff --git a/fdsreader/export/obst_exporter.py b/fdsreader/export/obst_exporter.py
index c894354f..582e68f8 100644
--- a/fdsreader/export/obst_exporter.py
+++ b/fdsreader/export/obst_exporter.py
@@ -37,7 +37,8 @@ def export_obst_raw(obst: Obstruction, output_dir: str, ordering: Literal["C", "
     meta["TimeSteps"] = len(random_bndf.times)
     faces = random_bndf.data
     meta["NumOrientations"] = len(faces)
-    # In the order of get_global_boundary_data_arrays, which the data files follow
+    # 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]
diff --git a/fdsreader/fds_classes/mesh.py b/fdsreader/fds_classes/mesh.py
index 0ac02359..99720166 100644
--- a/fdsreader/fds_classes/mesh.py
+++ b/fdsreader/fds_classes/mesh.py
@@ -4,7 +4,6 @@
 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
@@ -229,10 +228,6 @@ def _add_patches(
             for patch in patches:
                 patch._boundary_parent = self._boundary_data[bid]
 
-        if not settings.LAZY_LOAD:
-            for patch in patches:
-                _ = patch.data
-
     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 4b47a099..91c3435f 100644
--- a/fdsreader/fds_classes/mesh_collection.py
+++ b/fdsreader/fds_classes/mesh_collection.py
@@ -33,7 +33,7 @@ def get_global_boundary_data_arrays(
 
         :param quantity: The quantity's name or short name.
         """
-        return _global_boundary_arrays([mesh._boundary_data for mesh in self], quantity)
+        return _global_boundary_arrays([mesh.get_boundary_data for mesh in self], quantity)
 
     def __repr__(self):
         return "MeshCollection(" + super().__repr__() + ")"
diff --git a/tests/test_boundary_domain.py b/tests/test_boundary_domain.py
index c1872308..60c250de 100644
--- a/tests/test_boundary_domain.py
+++ b/tests/test_boundary_domain.py
@@ -24,19 +24,29 @@ def sim(request):
 
 def _floor_patches(sim, mesh_id, quantity):
     boundary = sim.meshes.get_by_id(mesh_id).get_boundary_data(quantity)
-    return sorted(
-        (p for p in boundary.patches if p.orientation == FLOOR), key=lambda p: (p.extent.x_start, p.extent.y_start)
-    )
+    return [p for p in boundary.patches if p.orientation == FLOOR]
 
 
-def test_all_patches_of_a_face_are_kept(sim):
+# 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)
-    assert len([p for p in boundary.patches if p.orientation == FLOOR]) == 4
-    assert boundary.orientations == [-2, 2, 3]  # x = 0 is open, x = 1 lies between the meshes, the top is open
-    assert boundary.vmax(FLOOR) == max(np.max(p.data) for p in _floor_patches(sim, "West", CELLS))
     with pytest.warns(UserWarning, match="Boundary.patches") as record:
         _ = boundary.data
-    assert record[0].filename == __file__  # points at the caller, not into fdsreader
+    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):
@@ -63,12 +73,13 @@ def test_the_floor_over_both_meshes(sim):
     np.testing.assert_array_equal(floor[-1] > 0, burning)
 
 
-def test_the_floor_holds_each_patch(sim):
+@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]
-    west = _floor_patches(sim, "West", CELLS)[0]  # x 0..0.6, y 0..1
-    np.testing.assert_array_equal(floor[:, 0:6, 0:10], west.data)
-    east = _floor_patches(sim, "East", CELLS)[0]  # x 1..2, y 0..0.4, 0.2 m cells
-    np.testing.assert_array_equal(floor[:, 10:20, 0:4], east.data.repeat(2, axis=1).repeat(2, axis=2))
+    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):