diff --git a/COMPAT.md b/COMPAT.md index 7a5441b0..98834c54 100644 --- a/COMPAT.md +++ b/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 | diff --git a/LoopStructural/modelling/core/_model_exporter.py b/LoopStructural/modelling/core/_model_exporter.py index 7098e048..1fbadb35 100644 --- a/LoopStructural/modelling/core/_model_exporter.py +++ b/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()) diff --git a/LoopStructural/modelling/features/_base_geological_feature.py b/LoopStructural/modelling/features/_base_geological_feature.py index 73956c21..cddd98da 100644 --- a/LoopStructural/modelling/features/_base_geological_feature.py +++ b/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, diff --git a/ROADMAP.md b/ROADMAP.md index b8f7df71..814e3c47 100644 --- a/ROADMAP.md +++ b/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` diff --git a/packages/loop_common/src/loop_common/geometry/_bounding_box.py b/packages/loop_common/src/loop_common/geometry/_bounding_box.py index b1a284a6..3b5f5f5a 100644 --- a/packages/loop_common/src/loop_common/geometry/_bounding_box.py +++ b/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, diff --git a/packages/loop_common/tests/test_bounding_box_structured_grid.py b/packages/loop_common/tests/test_bounding_box_structured_grid.py new file mode 100644 index 00000000..89f35002 --- /dev/null +++ b/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)