From 8dc57aec1830be7cc35bc167a8a89cafe9eea738 Mon Sep 17 00:00:00 2001 From: Marc Fehr Date: Thu, 15 Oct 2026 21:44:34 +0200 Subject: [PATCH 1/2] fix: global boundary arrays keep the time axis first and all face positions get_global_boundary_data_arrays squeezed every axis of length one, so faces one cell wide came out as (1, n_t, n) instead of (n_t, 1, n). It also decided whether the data is cell-centred from whichever quantity came first, which made a node-based quantity fail when the obstruction also had a cell-centred one, and it dropped the first of two face positions (an obstruction spanning meshes with different cell sizes, whose bounds FDS snaps to each grid). The arrays are now built from the patches of the requested quantity only, grouped by face position, and always have the shape (n_t, a, b). Several positions give a tuple ordered by position; export_obst_raw writes such faces part by part. The new case bndf_global (FDS 6.11.1) has two meshes with 0.1 m and 0.2 m cells, an obstruction across both and one a cell wide, with a node-based and a cell-centred quantity. On bndf_data the results equal the previous ones except for the swapped axes. Co-Authored-By: Claude Opus 5.5 --- .gitignore | 1 + fdsreader/bndf/obstruction.py | 232 +++++++++---------------- fdsreader/export/obst_exporter.py | 14 +- tests/cases/bndf_global_fds611.tgz | 3 + tests/cases/fds_inputs/bndf_global.fds | 19 ++ tests/test_boundary_global_arrays.py | 88 ++++++++++ 6 files changed, 201 insertions(+), 156 deletions(-) create mode 100644 tests/cases/bndf_global_fds611.tgz create mode 100644 tests/cases/fds_inputs/bndf_global.fds create mode 100644 tests/test_boundary_global_arrays.py diff --git a/.gitignore b/.gitignore index e0d6391d..9e73f54a 100644 --- a/.gitignore +++ b/.gitignore @@ -534,6 +534,7 @@ tests/cases/geomslice_data_fds*/ tests/cases/cface_data_fds*/ tests/cases/slcf_partial_fds*/ tests/cases/slcf_masks_fds*/ +tests/cases/bndf_global_fds*/ # Pickle cache files *.pickle diff --git a/fdsreader/bndf/obstruction.py b/fdsreader/bndf/obstruction.py index d82bd315..62131525 100644 --- a/fdsreader/bndf/obstruction.py +++ b/fdsreader/bndf/obstruction.py @@ -696,163 +696,45 @@ def has_boundary_data(self): def get_global_boundary_data_arrays( self, quantity: Union[str, Quantity] - ) -> Dict[int, Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]]: + ) -> Dict[int, Union[np.ndarray, Tuple[np.ndarray, ...]]]: """Creates a global numpy ndarray from all subobstruction's boundary data for each orientation. + Each array has the shape (n_t, a, b), where a and b are the two axes of the face in the order x, y, z (e.g. y + and z for the orientations -1 and 1). The grid uses the finest cell size of the meshes the obstruction spans; + values of coarser meshes are repeated, and points not covered by any patch are NaN. + + If the face lies at different positions in different meshes, e.g. because the obstruction spans meshes with + different cell sizes and FDS snaps its bounds to each grid, the orientation maps to a tuple with one array per + position, ordered by position. + :param quantity: The quantity's name or short name for which to load the global arrays. """ - if not self.has_boundary_data: - return dict() + patches_by_orientation: Dict[int, List[Patch]] = dict() for subobst in self._all_subobstructions: - if subobst.has_boundary_data: - cell_centered = next(iter(subobst._boundary_data.values())).cell_centered + if not subobst.has_boundary_data: + continue + 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_int in self.orientations: - subobst_sets = [list(), list()] - dim = ["x", "y", "z"][abs(orientation_int) - 1] - random_subobst = next( - (subobst for subobst in self._all_subobstructions if orientation_int in subobst.get_coordinates()), - None, - ) - if random_subobst is None: - raise KeyError(f"No sub-obstruction found for orientation {orientation_int}") - base_coord = random_subobst.get_coordinates(ignore_cell_centered=False)[orientation_int][dim][0] - - for subobst in self._all_subobstructions: - if not subobst.has_boundary_data: - continue - coords = subobst.get_coordinates(ignore_cell_centered=False) - if orientation_int not in coords: - continue - if coords[orientation_int][dim][0] == base_coord: - subobst_sets[0].append(subobst) - else: - subobst_sets[1].append(subobst) - if len(subobst_sets[1]) == 0: - subobst_sets.pop(1) - - first_result_grid = None - for subobsts in subobst_sets: - coord_min = {"x": math.inf, "y": math.inf, "z": math.inf} - coord_max = {"x": -math.inf, "y": -math.inf, "z": -math.inf} - for dim in ("x", "y", "z"): - for subobst in subobsts: - patch_extent = subobst.get_data(quantity).data[orientation_int].extent - co = subobst.get_coordinates(ignore_cell_centered=True)[orientation_int][dim] - co = co[ - np.where( - np.logical_and( - co >= patch_extent[dim][0], - co <= patch_extent[dim][1], - ) - ) - ] - coord_min[dim] = min(co[0], coord_min[dim]) - coord_max[dim] = max(co[-1], coord_max[dim]) - - # The global grid will use the finest mesh as base and duplicate values of the coarser meshes - # Therefore we first find the finest mesh and calculate the step size in each dimension - step_sizes_min = { - "x": coord_max["x"] - coord_min["x"], - "y": coord_max["y"] - coord_min["y"], - "z": coord_max["z"] - coord_min["z"], - } - step_sizes_max = {"x": 0, "y": 0, "z": 0} - steps = dict() - global_max = {"x": -math.inf, "y": -math.inf, "z": -math.inf} - - for dim in ("x", "y", "z"): - for subobst in subobsts: - subobst_coords = subobst.get_coordinates(ignore_cell_centered=True)[orientation_int] - if len(subobst_coords[dim]) <= 1: - step_size = 0 - else: - step_size = subobst_coords[dim][1] - subobst_coords[dim][0] - step_sizes_min[dim] = min(step_size, step_sizes_min[dim]) - step_sizes_max[dim] = max(step_size, step_sizes_max[dim]) - global_max[dim] = max(subobst_coords[dim][-1], global_max[dim]) - - for dim in ("x", "y", "z"): - if step_sizes_min[dim] == 0: - step_sizes_min[dim] = math.inf - steps[dim] = 1 - else: - steps[dim] = max(int(round((coord_max[dim] - coord_min[dim]) / step_sizes_min[dim])), 1) + ( - 0 if cell_centered else 1 - ) - - grid = np.full((self.n_t(count_duplicates=False), steps["x"], steps["y"], steps["z"]), np.nan) - - start_coordinates = {"x": coord_min["x"], "y": coord_min["y"], "z": coord_min["z"]} - start_idx = dict() - end_idx = dict() - for subobst in subobsts: - patch_data = np.expand_dims( - subobst.get_data(quantity).data[orientation_int].data, axis=abs(orientation_int) - ) - subobst_coords = subobst.get_coordinates(ignore_cell_centered=True)[orientation_int] - for axis in (0, 1, 2): - dim = ("x", "y", "z")[axis] - if axis == abs(orientation_int) - 1: - start_idx[dim] = 0 - end_idx[dim] = 1 - continue - n_repeat = max( - int(round((subobst_coords[dim][1] - subobst_coords[dim][0]) / step_sizes_min[dim])), 1 - ) - - start_idx[dim] = int( - round((subobst_coords[dim][0] - start_coordinates[dim]) / step_sizes_min[dim]) - ) - end_idx[dim] = int( - round((subobst_coords[dim][-1] - start_coordinates[dim]) / step_sizes_min[dim]) - ) - - # We ignore border points unless they are actually on the border of the simulation space as all - # other border points actually appear twice, as the subobstructions overlap. This only - # applies for face_centered data, as cell_centered data will not overlap. - if not cell_centered: - if axis != abs(orientation_int) - 1: - reduced_shape = list(patch_data.shape) - reduced_shape[axis + 1] -= 1 - reduced_data_slices = tuple(slice(s) for s in reduced_shape) - patch_data = patch_data[reduced_data_slices] - - # Temporarily save border points to add them back to the array again later - if subobst_coords[dim][-1] == global_max[dim]: - end_idx[dim] += 1 - temp_data_slices = [slice(s) for s in patch_data.shape] - temp_data_slices[axis + 1] = slice(patch_data.shape[axis + 1] - 1, None) - temp_data = patch_data[tuple(temp_data_slices)] - - if n_repeat > 1: - patch_data = np.repeat(patch_data, n_repeat, axis=axis + 1) - - # Add border points back again if needed - if not cell_centered and subobst_coords[dim][-1] == global_max[dim]: - patch_data = np.concatenate((patch_data, temp_data), axis=axis + 1) - - grid[ - :, start_idx["x"] : end_idx["x"], start_idx["y"] : end_idx["y"], start_idx["z"] : end_idx["z"] - ] = patch_data.reshape( - ( - self.n_t(count_duplicates=False), - end_idx["x"] - start_idx["x"], - end_idx["y"] - start_idx["y"], - end_idx["z"] - start_idx["z"], - ) - ) - - # Remove empty dimensions, but make sure to note remove the time dimension if there is only a single - # timestep - grid = np.squeeze(grid) - if len(grid.shape) == 2: - grid = grid[np.newaxis, :, :] - if first_result_grid is not None: - result[orientation_int] = (first_result_grid, grid) + for orientation, patches in patches_by_orientation.items(): + normal = ("x", "y", "z")[abs(orientation) - 1] + levels: List[Tuple[float, List[Patch]]] = list() + for patch in patches: + position = patch.extent[normal][0] + level = next((lvl for lvl in levels if math.isclose(lvl[0], position, abs_tol=1e-5)), None) + if level is None: + levels.append((position, [patch])) else: - result[orientation_int] = grid + level[1].append(patch) + grids = tuple(_global_patch_grid(level_patches, normal) for _, level_patches in sorted(levels)) + result[orientation] = grids[0] if len(grids) == 1 else grids return result def clear_cache(self): @@ -915,3 +797,55 @@ def __repr__(self, *args, **kwargs): + (f", Quantities={[q.short_name for q in self.quantities]}" if self.has_boundary_data else "") + ")" ) + + +def _patch_nodes(patch: Patch, dim: Literal["x", "y", "z"]) -> np.ndarray: + """The mesh nodes a patch spans along one axis.""" + co = patch.mesh.coordinates[dim] + lower, upper = patch.extent[dim] + tol = 1e-3 * float(np.min(np.diff(co))) if co.size > 1 else 1e-6 + return co[(co >= lower - tol) & (co <= upper + tol)] + + +def _global_patch_grid(patches: List[Patch], normal: Literal["x", "y", "z"]) -> np.ndarray: + """Combines patches lying in one plane into one array of shape (n_t, a, b) on the finest grid of the patches.""" + cell_centered = patches[0].cell_centered + axes = [dim for dim in ("x", "y", "z") if dim != normal] + nodes = {dim: [_patch_nodes(patch, dim) for patch in patches] for dim in axes} + + positions = dict() + for dim in axes: + lower = min(n[0] for n in nodes[dim]) + upper = max(n[-1] for n in nodes[dim]) + steps = [np.min(np.diff(n)) for n in nodes[dim] if n.size > 1] + step = min(steps) if steps else 0.0 + n_nodes = int(round((upper - lower) / step)) + 1 if step > 0 else 1 + positions[dim] = lower + step * np.arange(n_nodes) + if cell_centered: + positions[dim] = positions[dim][:-1] + step / 2 + + n_t = patches[0].data.shape[0] + grid = np.full((n_t,) + tuple(positions[dim].size for dim in axes), np.nan) + # On a node shared by two patches, the patch further up in x, y and z wins, as each is written after those below + order = sorted(range(len(patches)), key=lambda i: tuple(nodes[dim][i][0] for dim in axes)) + for i in order: + patch = patches[i] + global_indices = list() + patch_indices = list() + for dim in axes: + patch_nodes = nodes[dim][i] + tol = 1e-3 * (np.min(np.diff(patch_nodes)) if patch_nodes.size > 1 else 1.0) + pos = positions[dim] + if cell_centered: + inside = (pos > patch_nodes[0]) & (pos < patch_nodes[-1]) + n_values = patch_nodes.size - 1 + else: + inside = (pos >= patch_nodes[0] - tol) & (pos <= patch_nodes[-1] + tol) + n_values = patch_nodes.size + # Each point of the global grid takes the value of the patch's cell it lies in, or of the patch's last + # node at or below it, so values of coarser meshes are repeated + idx = np.searchsorted(patch_nodes, pos[inside] + tol, side="right") - 1 + global_indices.append(np.flatnonzero(inside)) + patch_indices.append(np.clip(idx, 0, n_values - 1)) + grid[(slice(None),) + np.ix_(*global_indices)] = patch.data[(slice(None),) + np.ix_(*patch_indices)] + return grid diff --git a/fdsreader/export/obst_exporter.py b/fdsreader/export/obst_exporter.py index c2b75a54..1b1a5084 100644 --- a/fdsreader/export/obst_exporter.py +++ b/fdsreader/export/obst_exporter.py @@ -56,9 +56,7 @@ def export_obst_raw(obst: Obstruction, output_dir: str, ordering: Literal["C", " } ) - def worker( - quantity: str, faces: Dict[int, Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]], vmin: float, vmax: float - ): + def worker(quantity: str, faces: Dict[int, Union[np.ndarray, Tuple[np.ndarray, ...]]], vmin: float, vmax: float): # Abort if no useful data is available if vmax <= 0: return @@ -76,10 +74,12 @@ def worker( with open(os.path.join(output_dir, quantity_name, filename), "wb") as rawfile: for face in faces.values(): # face for each orientation - for d in (face * out["ScaleFactor"]).astype(np.uint8): - if ordering == "F": - d = d.T - d.tofile(rawfile) + # A face lying at several positions (see get_global_boundary_data_arrays) is written part by part + for part in face if isinstance(face, tuple) else (face,): + for d in (part * out["ScaleFactor"]).astype(np.uint8): + if ordering == "F": + d = d.T + d.tofile(rawfile) with lock: meta["Quantities"].append(out) diff --git a/tests/cases/bndf_global_fds611.tgz b/tests/cases/bndf_global_fds611.tgz new file mode 100644 index 00000000..027951a7 --- /dev/null +++ b/tests/cases/bndf_global_fds611.tgz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6e502ffcf16979705835eb7fe011e7fba41bfae95ddf61390961a1a10fd331d1 +size 25790 diff --git a/tests/cases/fds_inputs/bndf_global.fds b/tests/cases/fds_inputs/bndf_global.fds new file mode 100644 index 00000000..73b9cc74 --- /dev/null +++ b/tests/cases/fds_inputs/bndf_global.fds @@ -0,0 +1,19 @@ +&HEAD CHID='bndf_global', TITLE='Global boundary data arrays for fdsreader' / +&MESH ID='Fine', IJK=10,10,10, XB=0.0,1.0,0.0,1.0,0.0,1.0 / +&MESH ID='Coarse', IJK=5,5,5, XB=1.0,2.0,0.0,1.0,0.0,1.0 / +&TIME T_END=2.0 / +&DUMP DT_BNDF=0.5, SMOKE3D=F / +&REAC FUEL='PROPANE' / +&SURF ID='FIRE', HRRPUA=500.0 / +&VENT XB=0.1,0.5,0.1,0.9,0.0,0.0, SURF_ID='FIRE' / +&VENT MB='XMIN', SURF_ID='OPEN' / &VENT MB='XMAX', SURF_ID='OPEN' / +&VENT MB='ZMAX', SURF_ID='OPEN' / + +Across both meshes: the top at z = 0.53 snaps to 0.5 in 'Fine' and to 0.6 in 'Coarse' +&OBST ID='Table', XB=0.6,1.4,0.2,0.8,0.0,0.53 / +One cell wide in x +&OBST ID='Strip', XB=0.2,0.3,0.2,0.8,0.0,0.3 / + +&BNDF QUANTITY='WALL TEMPERATURE' / +&BNDF QUANTITY='NET HEAT FLUX', CELL_CENTERED=T / +&TAIL / diff --git a/tests/test_boundary_global_arrays.py b/tests/test_boundary_global_arrays.py new file mode 100644 index 00000000..045abe53 --- /dev/null +++ b/tests/test_boundary_global_arrays.py @@ -0,0 +1,88 @@ +"""Obstruction.get_global_boundary_data_arrays on bndf_global_fds611 (FDS 6.11.1; input in +tests/cases/fds_inputs/bndf_global.fds). + +Two meshes: 'Fine' (x 0..1, 0.1 m cells) and 'Coarse' (x 1..2, 0.2 m cells). 'Table' spans both; its top at +z = 0.53 snaps to 0.5 in 'Fine' and to 0.6 in 'Coarse'. 'Strip' is one cell wide in x. Each face has the node-based +WALL TEMPERATURE and the cell-centred TOTAL HEAT FLUX. +""" + +import numpy as np +import pytest +from fds_cases import case_dirs, case_id, load + +NODES = "WALL TEMPERATURE" +CELLS = "TOTAL HEAT FLUX" + + +@pytest.fixture(scope="module", params=case_dirs("bndf_global_fds611", "bndf_global"), ids=case_id) +def sim(request): + return load(request.param) + + +def _patch(obst, mesh_id, quantity, orientation): + sub = next(sub for sub in obst._all_subobstructions if sub.mesh.id == mesh_id) + return sub.get_data(quantity).data[orientation].data + + +@pytest.mark.parametrize( + "orientation, nodes, cells", + [(-1, (7, 4), (6, 3)), (1, (7, 4), (6, 3)), (-2, (2, 4), (1, 3)), (2, (2, 4), (1, 3)), (3, (2, 7), (1, 6))], +) +def test_time_comes_first_on_faces_one_cell_wide(sim, orientation, nodes, cells): + strip = sim.obstructions.get_by_id("Strip") + n_t = strip.n_t(count_duplicates=False) + assert strip.get_global_boundary_data_arrays(NODES)[orientation].shape == (n_t,) + nodes + assert strip.get_global_boundary_data_arrays(CELLS)[orientation].shape == (n_t,) + cells + + +@pytest.mark.parametrize("quantity", [NODES, CELLS]) +def test_a_face_in_one_mesh_is_its_patch(sim, quantity): + strip = sim.obstructions.get_by_id("Strip") + for orientation, grid in strip.get_global_boundary_data_arrays(quantity).items(): + np.testing.assert_array_equal(grid, _patch(strip, "Fine", quantity, orientation)) + + +@pytest.mark.parametrize("quantity", [NODES, CELLS]) +def test_a_face_at_two_heights_keeps_both(sim, quantity): + table = sim.obstructions.get_by_id("Table") + lower, upper = table.get_global_boundary_data_arrays(quantity)[3] + np.testing.assert_array_equal(lower, _patch(table, "Fine", quantity, 3)) + np.testing.assert_array_equal(upper, _patch(table, "Coarse", quantity, 3)) + + +def test_coarse_values_are_repeated_on_the_finest_grid(sim): + # The face at y = 0.2 spans x 0.6..1.4 and z 0..0.5 in 'Fine', z 0..0.6 in 'Coarse': nodes every 0.1 m + table = sim.obstructions.get_by_id("Table") + grid = table.get_global_boundary_data_arrays(NODES)[-2] + fine = _patch(table, "Fine", NODES, -2) + coarse = _patch(table, "Coarse", NODES, -2) + assert grid.shape == (table.n_t(count_duplicates=False), 9, 7) + np.testing.assert_array_equal(grid[:, :4, :6], fine[:, :4, :]) + np.testing.assert_array_equal(grid[:, 5:7, 2], coarse[:, :2, 1]) # x 1.1 and 1.2 take the nodes 1.0 and 1.2 + np.testing.assert_array_equal(grid[:, 6, :], coarse[:, 1, :].repeat(2, axis=1)[:, :7]) # x 1.2, z 0.0..0.6 + assert np.isnan(grid[:, :4, 6]).all() # 'Fine' ends at z = 0.5 + + grid = table.get_global_boundary_data_arrays(CELLS)[-2] + coarse = _patch(table, "Coarse", CELLS, -2) + assert grid.shape == (table.n_t(count_duplicates=False), 8, 6) + np.testing.assert_array_equal(grid[:, 4:8, :], coarse.repeat(2, axis=1).repeat(2, axis=2)[:, :, :6]) + + +def test_node_based_and_cell_centred_quantities_on_one_obstruction(sim): + table = sim.obstructions.get_by_id("Table") + for quantity in (NODES, CELLS): + assert set(table.get_global_boundary_data_arrays(quantity)) == {-2, -1, 1, 2, 3} # no face on the floor + + +def test_unknown_quantity(sim): + with pytest.raises(KeyError): + sim.obstructions.get_by_id("Table").get_global_boundary_data_arrays("VELOCITY") + + +@pytest.mark.filterwarnings("ignore:invalid value encountered in cast:RuntimeWarning") +def test_export_of_a_face_at_two_heights(sim, tmp_path): + from fdsreader.export import export_obst_raw + + 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") From 567f78a930dc9e071c62cdf7d38991e8f1f2ca2d Mon Sep 17 00:00:00 2001 From: Marc Fehr Date: Thu, 15 Oct 2026 21:48:58 +0200 Subject: [PATCH 2/2] refactor: address the review of the global boundary arrays - Orientations come in sorted order, as before. - Meshes with fewer written frames leave the missing frames NaN instead of failing on the assignment. - Face positions and patch nodes are compared with a tolerance relative to the cell size; Patch.get_coordinates uses it too, so the grid code reuses it. - Document that the global grid is uniform. Co-Authored-By: Claude Opus 5.5 --- fdsreader/bndf/obstruction.py | 37 ++++++++++++++-------------- tests/test_boundary_global_arrays.py | 2 +- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/fdsreader/bndf/obstruction.py b/fdsreader/bndf/obstruction.py index 62131525..cb6ca7cc 100644 --- a/fdsreader/bndf/obstruction.py +++ b/fdsreader/bndf/obstruction.py @@ -92,7 +92,8 @@ def get_coordinates(self, ignore_cell_centered: bool = False) -> Dict[Literal["x if self.cell_centered and not ignore_cell_centered: co = co[:-1] + np.diff(co) / 2 - coords[dim] = co[np.where(np.logical_and(co >= self.extent[dim][0], co <= self.extent[dim][1]))] + tol = _tolerance(self.mesh.coordinates[dim]) + coords[dim] = co[(co >= self.extent[dim][0] - tol) & (co <= self.extent[dim][1] + tol)] if coords[dim].size == 0: coords[dim] = np.array([co[np.argmin(np.abs(co - self.extent[dim][0]))]]) @@ -701,7 +702,9 @@ def get_global_boundary_data_arrays( Each array has the shape (n_t, a, b), where a and b are the two axes of the face in the order x, y, z (e.g. y and z for the orientations -1 and 1). The grid uses the finest cell size of the meshes the obstruction spans; - values of coarser meshes are repeated, and points not covered by any patch are NaN. + values of coarser meshes are repeated, and points not covered by any patch are NaN. The grid is uniform, so + on stretched meshes values are placed at the nearest grid point below their node. Frames a mesh did not write + are NaN as well. If the face lies at different positions in different meshes, e.g. because the obstruction spans meshes with different cell sizes and FDS snaps its bounds to each grid, the orientation maps to a tuple with one array per @@ -711,8 +714,6 @@ def get_global_boundary_data_arrays( """ patches_by_orientation: Dict[int, List[Patch]] = dict() for subobst in self._all_subobstructions: - if not subobst.has_boundary_data: - continue try: boundary = subobst.get_data(quantity) except KeyError: @@ -723,17 +724,19 @@ def get_global_boundary_data_arrays( raise KeyError(f"No boundary data found for quantity '{quantity}'") result = dict() - for orientation, patches in patches_by_orientation.items(): + 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: + for patch in patches_by_orientation[orientation]: position = patch.extent[normal][0] - level = next((lvl for lvl in levels if math.isclose(lvl[0], position, abs_tol=1e-5)), None) + 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) - grids = tuple(_global_patch_grid(level_patches, normal) for _, level_patches in sorted(levels)) + 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 @@ -799,19 +802,16 @@ def __repr__(self, *args, **kwargs): ) -def _patch_nodes(patch: Patch, dim: Literal["x", "y", "z"]) -> np.ndarray: - """The mesh nodes a patch spans along one axis.""" - co = patch.mesh.coordinates[dim] - lower, upper = patch.extent[dim] - tol = 1e-3 * float(np.min(np.diff(co))) if co.size > 1 else 1e-6 - return co[(co >= lower - tol) & (co <= upper + tol)] +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 def _global_patch_grid(patches: List[Patch], normal: Literal["x", "y", "z"]) -> np.ndarray: """Combines patches lying in one plane into one array of shape (n_t, a, b) on the finest grid of the patches.""" cell_centered = patches[0].cell_centered axes = [dim for dim in ("x", "y", "z") if dim != normal] - nodes = {dim: [_patch_nodes(patch, dim) for patch in patches] for dim in axes} + nodes = {dim: [patch.get_coordinates(ignore_cell_centered=True)[dim] for patch in patches] for dim in axes} positions = dict() for dim in axes: @@ -824,7 +824,7 @@ def _global_patch_grid(patches: List[Patch], normal: Literal["x", "y", "z"]) -> if cell_centered: positions[dim] = positions[dim][:-1] + step / 2 - n_t = patches[0].data.shape[0] + n_t = max(patch.data.shape[0] for patch in patches) grid = np.full((n_t,) + tuple(positions[dim].size for dim in axes), np.nan) # On a node shared by two patches, the patch further up in x, y and z wins, as each is written after those below order = sorted(range(len(patches)), key=lambda i: tuple(nodes[dim][i][0] for dim in axes)) @@ -834,7 +834,7 @@ def _global_patch_grid(patches: List[Patch], normal: Literal["x", "y", "z"]) -> patch_indices = list() for dim in axes: patch_nodes = nodes[dim][i] - tol = 1e-3 * (np.min(np.diff(patch_nodes)) if patch_nodes.size > 1 else 1.0) + tol = _tolerance(patch_nodes) pos = positions[dim] if cell_centered: inside = (pos > patch_nodes[0]) & (pos < patch_nodes[-1]) @@ -847,5 +847,6 @@ def _global_patch_grid(patches: List[Patch], normal: Literal["x", "y", "z"]) -> idx = np.searchsorted(patch_nodes, pos[inside] + tol, side="right") - 1 global_indices.append(np.flatnonzero(inside)) patch_indices.append(np.clip(idx, 0, n_values - 1)) - grid[(slice(None),) + np.ix_(*global_indices)] = patch.data[(slice(None),) + np.ix_(*patch_indices)] + data = patch.data + grid[(slice(data.shape[0]),) + np.ix_(*global_indices)] = data[(slice(None),) + np.ix_(*patch_indices)] return grid diff --git a/tests/test_boundary_global_arrays.py b/tests/test_boundary_global_arrays.py index 045abe53..9c102a2b 100644 --- a/tests/test_boundary_global_arrays.py +++ b/tests/test_boundary_global_arrays.py @@ -71,7 +71,7 @@ def test_coarse_values_are_repeated_on_the_finest_grid(sim): def test_node_based_and_cell_centred_quantities_on_one_obstruction(sim): table = sim.obstructions.get_by_id("Table") for quantity in (NODES, CELLS): - assert set(table.get_global_boundary_data_arrays(quantity)) == {-2, -1, 1, 2, 3} # no face on the floor + assert list(table.get_global_boundary_data_arrays(quantity)) == [-2, -1, 1, 2, 3] # no face on the floor def test_unknown_quantity(sim):