refactor: BoundingBox.structured_grid returns a PropertyGrid

COMPAT.md

@@ -39,6 +39,7 @@ unchanged, but the new path is preferred going forward.
 | `LoopStructural.modelling.features._feature_converters.convert_feature_to_structural_frame` | `GeologicalModel.convert_feature_to_structural_frame` | 2026-07-24 |
 | `LoopStructural.geometry.StructuredGrid` / `loop_common.geometry.StructuredGrid` (results container with `properties`/`cell_properties`) | `PropertyGrid` (same module; `StructuredGrid` stays as an alias). `LoopStructural.geometry` now re-exports the `loop_common` class instead of carrying a duplicate copy. Requires `loop-common>=0.1.1,<0.2.0`; publish the updated shared package before releasing LoopStructural | 2026-10-29 |
 | `BoundingBox.cell_centres()` positions (`regular_grid(nsteps-1) + 0.5*step_vector`, which overshot the box by up to half a cell) | Midpoints of the `nsteps-1` cells between the `regular_grid()` points (always inside the box). Count unchanged. New `BoundingBox.node_step_vector` = `extent/(nsteps-1)`; `scalar_field`, `gradient_norm_scalar_field` and `get_block_model` grids now use it, so their values line up with their node and cell positions | 2026-10-29 |
+| `BoundingBox.structured_grid()` returned a `loop_common.supports.StructuredGrid` (the interpolation support: no `properties` dict, `cell_data` / `vertex_data` silently dropped, `nsteps + 1` nodes per axis) | Returns a `PropertyGrid` with `nsteps` nodes per axis spanning the box (the points of `regular_grid()`), keeping `cell_data` / `vertex_data`. `.nodes`, `.vtk()`, `.plot()`, `.properties`, `.cell_properties` work; the support's interpolation methods (`evaluate_value(points, array)`, `position_to_dof_coefs`, ...) are gone from this return value. For a grid to interpolate on build `loop_common.supports.StructuredGrid` from `step_vector` | 2026-10-29 |
 | `GeologicalModel.create_and_add_fault(..., faults=[])` | `GeologicalModel.create_and_add_fault(..., faults=None)` (list built internally, same effective default) | 2026-07-27 |
 | `GeologicalModel.create_and_add_intrusion(..., intrusion_frame_parameters={}, geometric_scaling_parameters={})` | `GeologicalModel.create_and_add_intrusion(..., intrusion_frame_parameters=None, geometric_scaling_parameters=None)` (dicts built internally, same effective default) | 2026-07-27 |
 | `GeologicalModel.get_fault_surfaces(faults=[])` | `GeologicalModel.get_fault_surfaces(faults=None)` (list built internally, same effective default) | 2026-07-27 |

LoopStructural/modelling/core/_model_exporter.py

@@ -71,15 +71,7 @@ def get_stratigraphic_surfaces(
     def get_block_model(
         model: GeologicalModel, name: str = 'block model'
     ) -> tuple[PropertyGrid, list]:
-        # NOTE: bounding_box.structured_grid() returns loop_common's
-        # interpolation-support StructuredGrid (no properties dict); use
-        # the PropertyGrid result container for storing values.
-        grid = PropertyGrid(
-            origin=model.bounding_box.origin,
-            step_vector=model.bounding_box.node_step_vector,
-            nsteps=model.bounding_box.nsteps,
-            name=name,
-        )
+        grid = model.bounding_box.structured_grid(name=name)
 
         grid.cell_properties['stratigraphy'] = model.evaluate_model(
             model.rescale(model.bounding_box.cell_centres())

LoopStructural/modelling/features/_base_geological_feature.py

@@ -4,7 +4,7 @@
 
 import numpy as np
 
-from LoopStructural.geometry import PropertyGrid, VectorPoints
+from LoopStructural.geometry import VectorPoints
 from LoopStructural.modelling.features import FeatureType
 from LoopStructural.utils import LoopIsosurfacer, LoopValueError, getLogger, surface_list
 from LoopStructural.utils.typing import NumericInput
@@ -384,15 +384,7 @@ def scalar_field(self, bounding_box=None):
             if self.model is None:
                 raise ValueError("Must specify bounding box")
             bounding_box = self.model.bounding_box
-        # NOTE: bounding_box.structured_grid() returns loop_common's
-        # interpolation-support StructuredGrid (no properties dict); use
-        # the PropertyGrid result container for storing values.
-        grid = PropertyGrid(
-            origin=bounding_box.origin,
-            step_vector=bounding_box.node_step_vector,
-            nsteps=bounding_box.nsteps,
-            name=self.name,
-        )
+        grid = bounding_box.structured_grid(name=self.name)
         value = self.evaluate_value(bounding_box.regular_grid(local=False, order='F'))
         grid.properties[self.name] = value
 
@@ -417,12 +409,7 @@ def gradient_norm_scalar_field(self, bounding_box=None):
             if self.model is None:
                 raise ValueError("Must specify bounding box")
             bounding_box = self.model.bounding_box
-        grid = PropertyGrid(
-            origin=bounding_box.origin,
-            step_vector=bounding_box.node_step_vector,
-            nsteps=bounding_box.nsteps,
-            name=self.name,
-        )
+        grid = bounding_box.structured_grid(name=self.name)
         value = np.linalg.norm(
             self.evaluate_gradient(bounding_box.regular_grid(local=False, order='F')),
             axis=1,

ROADMAP.md

@@ -794,7 +794,8 @@ just at release time.
   `loop_common.supports.StructuredGrid` (an interpolation support object with
   no properties dict), not the LoopStructural geometry `StructuredGrid`
   dataclass `scalar_field`/`get_block_model` actually need -- those call
-  sites now construct `LoopStructural.geometry.StructuredGrid` directly.
+  sites now construct `LoopStructural.geometry.StructuredGrid` directly. (Superseded: `bounding_box.structured_grid()` now returns a `PropertyGrid`, which
+  holds the properties; the interpolation support is built separately from `step_vector`.)
   Fixed a 2D-bounding-box regression in the new local-frame support-building
   code: `BoundingBox.corners` is 3D-only, so `create_support_from_bbox` and
   `structured_grid(local_coordinates=True)` now project `origin`/`maximum`

packages/loop_common/src/loop_common/geometry/_bounding_box.py

@@ -9,6 +9,8 @@
 # from LoopStructural.utils.exceptions import LoopValueError
 from loop_common.math import rng
 
+from ._structured_grid import PropertyGrid
+
 logger = getLogger(__name__)
 
 
@@ -699,12 +701,31 @@ def structured_grid(
         vertex_data: dict | None = None,
         name: str = "bounding_box",
         local_coordinates: bool = False,
-    ):
-        # imported here because loop_common.supports builds on loop_common.geometry
-        from loop_common.supports import StructuredGrid
+    ) -> PropertyGrid:
+        """Create a :class:`PropertyGrid` whose nodes are the points of :meth:`regular_grid`.
+
+        The grid has ``nsteps`` nodes in each direction spanning the bounding box
+        (spacing :attr:`node_step_vector`), so values computed at ``regular_grid()`` /
+        ``cell_centres()`` can be stored on it directly. For a grid to interpolate on, use
+        ``loop_common.supports.StructuredGrid`` (built from ``step_vector``, the cell size).
 
+        Parameters
+        ----------
+        cell_data : dict, optional
+            cell properties to store on the grid (copied)
+        vertex_data : dict, optional
+            node properties to store on the grid (copied)
+        name : str, optional
+            name of the grid, by default "bounding_box"
+        local_coordinates : bool, optional
+            build the grid in the bounding box's local frame, by default False
+
+        Returns
+        -------
+        PropertyGrid
+        """
         # python is passing a reference to the cell_data, vertex_data dicts so we need to
-        # copy them to make sure that different instances of StructuredGrid are not sharing the same
+        # copy them to make sure that different instances of PropertyGrid are not sharing the same
         # underlying objects
         if cell_data is None:
             cell_data = {}
@@ -721,12 +742,12 @@ def structured_grid(
             local_points = self.project(self.corners)
             local_origin = np.min(local_points, axis=0)
             local_maximum = np.max(local_points, axis=0)
-            step_vector = (local_maximum - local_origin) / self.nsteps
+            step_vector = (local_maximum - local_origin) / (self.nsteps - 1)
             origin = local_origin
         else:
-            step_vector = self.step_vector
+            step_vector = self.node_step_vector
             origin = self.origin
-        return StructuredGrid(
+        return PropertyGrid(
             origin=origin,
             step_vector=step_vector,
             nsteps=self.nsteps,

packages/loop_common/tests/test_bounding_box_structured_grid.py

@@ -0,0 +1,71 @@
+"""BoundingBox.structured_grid returns a PropertyGrid whose nodes are the regular_grid points."""
+import numpy as np
+import pytest
+from loop_common.geometry import BoundingBox, PropertyGrid
+
+BOXES = {
+    "unit": dict(origin=[0.0, 0.0, 0.0], maximum=[1.0, 1.0, 1.0], nsteps=[10, 10, 10]),
+    "offset_anisotropic": dict(origin=[100.0, -20.0, 5.0], maximum=[112.0, 4.0, 35.0], nsteps=[7, 5, 4]),
+}
+
+
+@pytest.fixture(params=list(BOXES))
+def box(request):
+    return BoundingBox(**BOXES[request.param])
+
+
+def test_returns_a_property_grid(box):
+    assert isinstance(box.structured_grid(), PropertyGrid)
+
+
+def test_nodes_are_the_regular_grid_points(box):
+    grid = box.structured_grid()
+    np.testing.assert_allclose(grid.nodes, box.regular_grid(local=False, order="F"))
+    np.testing.assert_array_equal(grid.nsteps, box.nsteps)
+
+
+def test_grid_spans_the_bounding_box(box):
+    grid = box.structured_grid()
+    np.testing.assert_allclose(grid.origin, box.origin)
+    np.testing.assert_allclose(grid.maximum, box.maximum)
+
+
+def test_cell_centres_are_the_bounding_box_cell_centres(box):
+    grid = box.structured_grid()
+    np.testing.assert_allclose(grid.cell_centres, box.reproject(box.cell_centres(order="F")))
+
+
+def test_data_is_stored_and_copied(box):
+    nodes = np.arange(float(box.nsteps.prod()))
+    cells = np.arange(float((box.nsteps - 1).prod()))
+    node_data = {"n": nodes}
+    grid = box.structured_grid(cell_data={"c": cells}, vertex_data=node_data, name="g")
+    assert grid.name == "g"
+    np.testing.assert_array_equal(grid.properties["n"], nodes)
+    np.testing.assert_array_equal(grid.cell_properties["c"], cells)
+    grid.properties["other"] = nodes
+    assert "other" not in node_data
+    other = box.structured_grid()
+    assert other.properties == {} and other.cell_properties == {}
+
+
+def test_values_computed_on_the_box_land_on_the_grid_nodes(box):
+    """The property that matters for scalar_field: a value computed at regular_grid() is at that node."""
+    grid = box.structured_grid()
+    grid.properties["x"] = box.regular_grid(local=False, order="F")[:, 0]
+    np.testing.assert_allclose(grid.properties["x"], grid.nodes[:, 0])
+    np.testing.assert_allclose(grid.evaluate_value("x", grid.nodes), grid.nodes[:, 0])
+
+
+def test_local_coordinates_grid_contains_the_rotated_box():
+    box = BoundingBox(origin=[0.0, 0.0, 0.0], maximum=[10.0, 10.0, 10.0], nsteps=[5, 5, 5])
+    angle = np.radians(45.0)
+    rotation = np.array(
+        [[np.cos(angle), -np.sin(angle), 0.0], [np.sin(angle), np.cos(angle), 0.0], [0.0, 0.0, 1.0]]
+    )
+    box.set_local_transform(local_origin=[0.0, 0.0, 0.0], rotation_matrix=rotation)
+    corners = box.project(box.corners)
+    grid = box.structured_grid(local_coordinates=True)
+    np.testing.assert_allclose(grid.origin, corners.min(axis=0), atol=1e-8)
+    np.testing.assert_allclose(grid.maximum, corners.max(axis=0), atol=1e-8)
+    np.testing.assert_array_equal(grid.nsteps, box.nsteps)