fix: make scalar_field / block-model grids line up with their values

COMPAT.md

@@ -38,6 +38,7 @@ unchanged, but the new path is preferred going forward.
 | `LoopStructural.modelling.features._feature_converters.add_fold_to_feature` (private module, imported directly by the QGIS plugin) | `GeologicalModel.add_fold_to_feature` | 2026-07-24 |
 | `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 |
 | `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

@@ -76,7 +76,7 @@ def get_block_model(
         # the PropertyGrid result container for storing values.
         grid = PropertyGrid(
             origin=model.bounding_box.origin,
-            step_vector=model.bounding_box.step_vector,
+            step_vector=model.bounding_box.node_step_vector,
             nsteps=model.bounding_box.nsteps,
             name=name,
         )

LoopStructural/modelling/features/_base_geological_feature.py

@@ -389,7 +389,7 @@ def scalar_field(self, bounding_box=None):
         # the PropertyGrid result container for storing values.
         grid = PropertyGrid(
             origin=bounding_box.origin,
-            step_vector=bounding_box.step_vector,
+            step_vector=bounding_box.node_step_vector,
             nsteps=bounding_box.nsteps,
             name=self.name,
         )
@@ -419,7 +419,7 @@ def gradient_norm_scalar_field(self, bounding_box=None):
             bounding_box = self.model.bounding_box
         grid = PropertyGrid(
             origin=bounding_box.origin,
-            step_vector=bounding_box.step_vector,
+            step_vector=bounding_box.node_step_vector,
             nsteps=bounding_box.nsteps,
             name=self.name,
         )

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

@@ -30,6 +30,20 @@ def __init__(
         """A bounding box for a model, defined by the
         origin, maximum and number of steps in each direction
 
+        Notes
+        -----
+        ``nsteps`` is used in two ways, so two spacings are provided:
+
+        * Cell-based (interpolation): ``nsteps`` is the number of cells in each
+          direction and ``step_vector`` is ``(maximum - origin) / nsteps``. Interpolation
+          supports are built from these.
+        * Node-based (values on a grid): ``regular_grid`` returns ``nsteps`` points spanning
+          the full extent, so the spacing between them is ``node_step_vector``,
+          ``(maximum - origin) / (nsteps - 1)``. ``cell_centres`` returns the
+          ``nsteps - 1`` centres of the cells between those points. Use this
+          convention (as pyvista/VTK ``ImageData`` and ``numpy.linspace`` do)
+          when storing values computed at ``regular_grid`` or ``cell_centres``.
+
         Parameters
         ----------
         dimensions : int, optional
@@ -393,6 +407,18 @@ def step_vector(self):
             raise LoopValueError("Cannot compute step_vector: nsteps contains zero values")
         return (self.maximum - self.origin) / self.nsteps
 
+    @property
+    def node_step_vector(self):
+        """Spacing between the ``nsteps`` points returned by :meth:`regular_grid`.
+
+        ``(maximum - origin) / (nsteps - 1)``; compare ``step_vector``, which divides
+        by ``nsteps`` and is the cell size used to build interpolation supports.
+        """
+        nodes = np.asarray(self.nsteps, dtype=float)
+        if np.any(nodes < 2):
+            raise LoopValueError("Cannot compute node_step_vector: nsteps must be at least 2")
+        return (self.maximum - self.origin) / (nodes - 1)
+
     @property
     def length(self):
         return self.maximum - self.origin
@@ -573,21 +599,29 @@ def regular_grid(
         return locs
 
     def cell_centres(self, order: str = "F") -> np.ndarray:
-        """Get the cell centres of a regular grid
+        """Get the centres of the cells between the points of :meth:`regular_grid`
+
+        There are ``nsteps - 1`` cells in each direction, each ``node_step_vector`` wide,
+        so every centre is inside the bounding box.
 
         Parameters
         ----------
         order : str, optional
-            order of the grid, by default "C"
+            order of the grid, by default "F"
 
         Returns
         -------
         np.ndarray
-            array of cell centres
+            array of cell centres, in the local coordinate system
         """
-        locs = self.regular_grid(order=order, nsteps=self.nsteps - 1)
-
-        return locs + 0.5 * self.step_vector
+        half = 0.5 * self.node_step_vector
+        coordinates = [
+            np.linspace(self.origin[i] + half[i], self.maximum[i] - half[i], self.nsteps[i] - 1)
+            for i in range(self.dimensions)
+        ]
+        coordinate_grid = np.meshgrid(*coordinates, indexing="ij")
+        locs = np.array([coord.flatten(order=order) for coord in coordinate_grid]).T
+        return self.project(locs)
 
     def to_dict(self) -> dict:
         """Export the defining characteristics of the bounding

packages/loop_common/src/loop_common/geometry/_structured_grid.py

@@ -131,7 +131,7 @@ def cell_centres(self):
         """
         x = np.linspace(
             self.origin[0] + self.step_vector[0] * 0.5,
-            self.maximum[0] + self.step_vector[0] * 0.5,
+            self.maximum[0] - self.step_vector[0] * 0.5,
             self.nsteps[0] - 1,
         )
         y = np.linspace(

tests/unit/modelling/test_scalar_field_grid.py

@@ -0,0 +1,96 @@
+"""The grid returned by ``feature.scalar_field`` must be registered to the values stored on it.
+
+``scalar_field`` evaluates the feature at ``bounding_box.regular_grid()`` /
+``bounding_box.cell_centres()`` and stores the results on a grid, so the grid's node
+and cell positions must be the points the values were computed at. ``regular_grid`` puts
+``nsteps`` points across the full extent (spacing ``node_step_vector``); ``step_vector`` is the
+cell size used for interpolation supports and must not be used to describe that grid.
+
+The feature is the analytical plane ``f(x, y, z) = x``, so the right value at any location is
+its x coordinate and no interpolation is involved.
+"""
+import numpy as np
+import pytest
+
+from LoopStructural.geometry import BoundingBox
+from LoopStructural.modelling.features._analytical_feature import (
+    AnalyticalGeologicalFeature,
+)
+
+
+@pytest.fixture
+def feature():
+    # value == x coordinate
+    return AnalyticalGeologicalFeature(
+        name="x",
+        vector=np.array([1.0, 0.0, 0.0]),
+        origin=np.array([0.0, 0.0, 0.0]),
+    )
+
+
+@pytest.fixture
+def bounding_box():
+    return BoundingBox(origin=[0.0, 0.0, 0.0], maximum=[10.0, 10.0, 10.0], nsteps=[10, 10, 10])
+
+
+def test_scalar_field_values_are_at_the_node_positions(feature, bounding_box):
+    grid = feature.scalar_field(bounding_box)
+    nodes = np.asarray(grid.nodes)
+    np.testing.assert_allclose(grid.properties["x"], nodes[:, 0], atol=1e-8)
+
+
+def test_scalar_field_grid_spans_the_bounding_box(feature, bounding_box):
+    grid = feature.scalar_field(bounding_box)
+    nodes = np.asarray(grid.nodes)
+    np.testing.assert_allclose(nodes.min(axis=0), bounding_box.origin, atol=1e-8)
+    np.testing.assert_allclose(nodes.max(axis=0), bounding_box.maximum, atol=1e-8)
+
+
+def test_scalar_field_cell_values_are_at_the_cell_centres(feature, bounding_box):
+    grid = feature.scalar_field(bounding_box)
+    centres = np.asarray(grid.cell_centres)
+    np.testing.assert_allclose(grid.cell_properties["x"], centres[:, 0], atol=1e-8)
+
+
+def test_cell_centres_are_inside_the_bounding_box(bounding_box):
+    centres = bounding_box.cell_centres()
+    assert np.all(centres >= bounding_box.origin - 1e-8)
+    assert np.all(centres <= bounding_box.maximum + 1e-8)
+
+
+@pytest.fixture
+def offset_box():
+    # non-zero origin and different node counts per axis: catches mixed-up axes and offsets
+    return BoundingBox(origin=[100.0, -20.0, 5.0], maximum=[112.0, 4.0, 35.0], nsteps=[7, 5, 4])
+
+
+def test_node_step_vector_is_the_spacing_of_regular_grid(offset_box):
+    points = offset_box.regular_grid(local=False)
+    for axis in range(3):
+        coords = np.unique(np.round(points[:, axis], 10))
+        assert len(coords) == offset_box.nsteps[axis]
+        np.testing.assert_allclose(np.diff(coords), offset_box.node_step_vector[axis])
+
+
+def test_cell_centres_count_and_midpoints(offset_box):
+    centres = offset_box.cell_centres()
+    assert centres.shape == (np.prod(offset_box.nsteps - 1), 3)
+    for axis in range(3):
+        coords = np.unique(np.round(centres[:, axis], 10))
+        nodes = np.linspace(
+            offset_box.origin[axis], offset_box.maximum[axis], offset_box.nsteps[axis]
+        )
+        np.testing.assert_allclose(coords, 0.5 * (nodes[:-1] + nodes[1:]))
+
+
+def test_scalar_field_on_offset_anisotropic_box(feature, offset_box):
+    feature = AnalyticalGeologicalFeature(
+        name="f",
+        vector=np.array([1.0, 2.0, 3.0]),
+        origin=np.array([0.0, 0.0, 0.0]),
+    )
+    grid = feature.scalar_field(offset_box)
+    expected_nodes = feature.evaluate_value(np.asarray(grid.nodes))
+    np.testing.assert_allclose(grid.properties["f"], expected_nodes, atol=1e-8)
+    expected_cells = feature.evaluate_value(np.asarray(grid.cell_centres))
+    np.testing.assert_allclose(grid.cell_properties["f"], expected_cells, atol=1e-8)