From 47158ebd5f48c318fb7414c3114f1e784fae0a58 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 13:37:49 +1030 Subject: [PATCH] feat: trilinear interpolation on PropertyGrid, via StructuredGrid3DGeometry Move the pure-geometry pieces of the interpolation support into StructuredGrid3DGeometry so the support and PropertyGrid share them: trilinear weights, position_to_local_coordinates, position_to_dof_coefs, position_to_cell_corners / _vertices, node_indexes_to_position, element_size / element_scale and evaluate_value (now also accepting (n_nodes, k) data). PropertyGrid.evaluate_value(name, points) interpolates a node property at arbitrary points (NaN outside), which is what is needed to transfer a property onto another support. The structured-support reference test is unchanged and passes. Co-Authored-By: Claude Sonnet 5.5 --- .../loop_common/geometry/_structured_grid.py | 21 +++ .../geometry/_structured_grid_3d.py | 175 ++++++++++++++++++ .../supports/_3d_base_structured.py | 64 ------- .../supports/_3d_structured_grid.py | 103 ----------- .../loop_common/tests/test_property_grid.py | 72 +++++++ 5 files changed, 268 insertions(+), 167 deletions(-) diff --git a/packages/loop_common/src/loop_common/geometry/_structured_grid.py b/packages/loop_common/src/loop_common/geometry/_structured_grid.py index 3b2a8de2e..d0cb41b77 100644 --- a/packages/loop_common/src/loop_common/geometry/_structured_grid.py +++ b/packages/loop_common/src/loop_common/geometry/_structured_grid.py @@ -131,6 +131,27 @@ def cell_index_at(self, points): index[inside] = geometry.global_cell_indices(cell_ijk[inside]) return index + def evaluate_value(self, name, points): + """Trilinear interpolation of the node property ``name`` at the points. + + This is what is needed to transfer a property onto another support or mesh: pass the + other support's node coordinates as ``points``. + + Parameters + ---------- + name : str + key of ``properties``; values are per node, in node order (x varies fastest), either + ``(n_nodes,)`` or ``(n_nodes, k)`` + points : array-like + N x 3 locations (or one location) in the grid's own frame + + Returns + ------- + np.ndarray + ``(N,)`` (or ``(N, k)``) values; NaN for points outside the grid + """ + return self.geometry.evaluate_value(np.asarray(points, dtype=float), self.properties[name]) + def cell_value_at(self, name, points): """Value of the cell property ``name`` in the cell each point falls in. diff --git a/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py b/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py index 4dd39258c..b78ac033c 100644 --- a/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py +++ b/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py @@ -265,3 +265,178 @@ def global_node_indices(self, node_indexes): def global_cell_indices(self, cell_indexes): """Flat cell indexes from (..., 3) cell ijk indexes; the order cell data is stored in.""" return self._global_indicies(cell_indexes, self.nsteps_cells) + + def position_to_cell_corners(self, pos): + """Get the global indices of the vertices (corners) of the cell containing each point. + + Parameters + ---------- + pos : np.array + (N, 3) array of xyz coordinates representing the positions of N points. + + Returns + ------- + globalidx : np.array + (N, 8) array of global indices corresponding to the 8 corner nodes of the cell + each point lies in. If a point lies outside the support, its corresponding entry + will be set to -1. + inside : np.array + (N,) boolean array indicating whether each point is inside the support domain. + """ + cell_indexes, inside = self.position_to_cell_index(pos) + nx, ny = self.nsteps[0], self.nsteps[1] + offsets = np.array( + [0, 1, nx, nx + 1, nx * ny, nx * ny + 1, nx * ny + nx, nx * ny + nx + 1], + dtype=np.intp, + ) + g = cell_indexes[:, 0] + nx * cell_indexes[:, 1] + nx * ny * cell_indexes[:, 2] + globalidx = g[:, None] + offsets[None, :] # (N, 8) + globalidx[~inside] = -1 + return globalidx, inside + + def position_to_cell_vertices(self, pos): + """Get the vertices of the cell a point is in + + Parameters + ---------- + pos : np.array + Nx3 array of xyz locations + + Returns + ------- + np.array((N,3),dtype=float), np.array(N,dtype=int) + vertices, inside + """ + gi, inside = self.position_to_cell_corners(pos) + node_indexes = self.global_index_to_node_index(gi) + return self.node_indexes_to_position(node_indexes), inside + + def node_indexes_to_position(self, node_indexes: np.ndarray) -> np.ndarray: + original_shape = node_indexes.shape + node_indexes = node_indexes.reshape((-1, 3)) + xyz = np.zeros((node_indexes.shape[0], 3), dtype=float) + xyz[:, 0] = self.origin[0] + self.step_vector[0] * node_indexes[:, 0] + xyz[:, 1] = self.origin[1] + self.step_vector[1] * node_indexes[:, 1] + xyz[:, 2] = self.origin[2] + self.step_vector[2] * node_indexes[:, 2] + xyz = xyz.reshape(original_shape) + return xyz + + @property + def element_size(self): + return np.prod(self.step_vector) + + @property + def element_scale(self): + # all elements are the same size + return 1.0 + + def trilinear(self, local_coords): + """ + Returns the trilinear interpolation for the local coordinates + + Parameters + ---------- + x - double, array of doubles + y - double, array of doubles + z - double, array of doubles + + Returns + ------- + array of interpolation coefficients + + """ + interpolant = np.zeros((local_coords.shape[0], 8), dtype=np.float64) + # interpolant[:,0] = (1 - local_coords[:, 0])* (1 - local_coords[:, 1])* (1 - local_coords[:, 2]) + # interpolant[:,1] = local_coords[:, 0]* (1 - local_coords[:, 1])* (1 - local_coords[:, 2]) + # interpolant[:,2] = (1 - local_coords[:, 0])* local_coords[:, 1]* (1 - local_coords[:, 2]) + # interpolant[:,3] = (1 - local_coords[:, 0])* (1 - local_coords[:, 1])* local_coords[:, 2] + # interpolant[:,4] = local_coords[:, 0] * (1 - local_coords[:, 1]) * local_coords[:, 2] + # interpolant[:,5] = (1 - local_coords[:, 0]) * local_coords[:, 1] * local_coords[:, 2] + # interpolant[:,6] = local_coords[:, 0] * local_coords[:, 1] * (1 - local_coords[:, 2]) + # interpolant[:,7] = local_coords[:, 0] * local_coords[:, 1] * local_coords[:, 2] + interpolant[:, 0] = ( + (1 - local_coords[:, 0]) * (1 - local_coords[:, 1]) * (1 - local_coords[:, 2]) + ) + interpolant[:, 1] = local_coords[:, 0] * (1 - local_coords[:, 1]) * (1 - local_coords[:, 2]) + interpolant[:, 2] = (1 - local_coords[:, 0]) * local_coords[:, 1] * (1 - local_coords[:, 2]) + interpolant[:, 4] = (1 - local_coords[:, 0]) * (1 - local_coords[:, 1]) * local_coords[:, 2] + interpolant[:, 5] = local_coords[:, 0] * (1 - local_coords[:, 1]) * local_coords[:, 2] + interpolant[:, 6] = (1 - local_coords[:, 0]) * local_coords[:, 1] * local_coords[:, 2] + interpolant[:, 3] = local_coords[:, 0] * local_coords[:, 1] * (1 - local_coords[:, 2]) + interpolant[:, 7] = local_coords[:, 0] * local_coords[:, 1] * local_coords[:, 2] + return interpolant + + def position_to_local_coordinates(self, pos): + """ + Convert from global to local coordinates within a cel + + Parameters + ---------- + pos - array of positions inside + + Returns + ------- + localx, localy, localz + + """ + pos = self.check_position(pos) + cell_indexes, _inside = self.position_to_cell_index(pos) + + # Use cell origin (not modulo) so points on max faces/vertices map to local=1. + cell_origin = self.origin[None, :] + self.step_vector[None, :] * cell_indexes + local_coords = (pos - cell_origin) / self.step_vector[None, :] + local_coords = np.clip(local_coords, 0.0, 1.0) + return local_coords + + def position_to_dof_coefs(self, pos): + """ + Global posotion to interpolation coefficients + + Parameters + ---------- + pos + + Returns + ------- + + """ + local_coords = self.position_to_local_coordinates(pos) + weights = self.trilinear(local_coords) + return weights + + def evaluate_value(self, evaluation_points, property_array): + """Trilinear interpolation of node values at the evaluation points. + + Parameters + ---------- + evaluation_points : np.ndarray + N x 3 locations in the grid's own frame + property_array : np.ndarray + values at the nodes, shape ``(n_nodes,)`` or ``(n_nodes, k)``, in node order + (x varies fastest) + + Returns + ------- + np.ndarray + ``(N,)`` (or ``(N, k)``) interpolated values; NaN for points outside the grid + """ + property_array = np.asarray(property_array) + evaluation_points = self.check_position(evaluation_points) + if property_array.shape[0] != self.n_nodes: + logger.error("Property array does not match grid") + raise ValueError( + f"cannot assign {property_array.shape[0]} vlaues to array of shape {self.n_nodes}" + ) + idc, inside = self.position_to_cell_corners(evaluation_points) + if idc.shape[0] != inside.shape[0]: + raise ValueError("index does not match number of nodes") + weights = self.position_to_dof_coefs(evaluation_points[inside, :]) + if property_array.ndim == 1: + v = np.zeros(idc.shape) + v[:, :] = np.nan + v[inside, :] = weights + v[inside, :] *= property_array[idc[inside, :]] + return np.sum(v, axis=1) + result = np.full((idc.shape[0],) + property_array.shape[1:], np.nan) + result[inside] = np.einsum("nc,nc...->n...", weights, property_array[idc[inside, :]]) + return result diff --git a/packages/loop_common/src/loop_common/supports/_3d_base_structured.py b/packages/loop_common/src/loop_common/supports/_3d_base_structured.py index 55f606c1f..113797458 100644 --- a/packages/loop_common/src/loop_common/supports/_3d_base_structured.py +++ b/packages/loop_common/src/loop_common/supports/_3d_base_structured.py @@ -116,70 +116,6 @@ def check_position(self, pos: np.ndarray) -> np.ndarray: raise ValueError("Position array needs to be a list of points or a point") return pos - def position_to_cell_corners(self, pos): - """Get the global indices of the vertices (corners) of the cell containing each point. - - Parameters - ---------- - pos : np.array - (N, 3) array of xyz coordinates representing the positions of N points. - - Returns - ------- - globalidx : np.array - (N, 8) array of global indices corresponding to the 8 corner nodes of the cell - each point lies in. If a point lies outside the support, its corresponding entry - will be set to -1. - inside : np.array - (N,) boolean array indicating whether each point is inside the support domain. - """ - cell_indexes, inside = self.position_to_cell_index(pos) - nx, ny = self.nsteps[0], self.nsteps[1] - offsets = np.array( - [0, 1, nx, nx + 1, nx * ny, nx * ny + 1, nx * ny + nx, nx * ny + nx + 1], - dtype=np.intp, - ) - g = cell_indexes[:, 0] + nx * cell_indexes[:, 1] + nx * ny * cell_indexes[:, 2] - globalidx = g[:, None] + offsets[None, :] # (N, 8) - globalidx[~inside] = -1 - return globalidx, inside - - def position_to_cell_vertices(self, pos): - """Get the vertices of the cell a point is in - - Parameters - ---------- - pos : np.array - Nx3 array of xyz locations - - Returns - ------- - np.array((N,3),dtype=float), np.array(N,dtype=int) - vertices, inside - """ - gi, inside = self.position_to_cell_corners(pos) - node_indexes = self.global_index_to_node_index(gi) - return self.node_indexes_to_position(node_indexes), inside - - def node_indexes_to_position(self, node_indexes: np.ndarray) -> np.ndarray: - original_shape = node_indexes.shape - node_indexes = node_indexes.reshape((-1, 3)) - xyz = np.zeros((node_indexes.shape[0], 3), dtype=float) - xyz[:, 0] = self.origin[0] + self.step_vector[0] * node_indexes[:, 0] - xyz[:, 1] = self.origin[1] + self.step_vector[1] * node_indexes[:, 1] - xyz[:, 2] = self.origin[2] + self.step_vector[2] * node_indexes[:, 2] - xyz = xyz.reshape(original_shape) - return xyz - - @property - def element_size(self): - return np.prod(self.step_vector) - - @property - def element_scale(self): - # all elements are the same size - return 1.0 - def vtk(self, node_properties=None, cell_properties=None): try: import pyvista as pv diff --git a/packages/loop_common/src/loop_common/supports/_3d_structured_grid.py b/packages/loop_common/src/loop_common/supports/_3d_structured_grid.py index e11e5daf6..77d20e207 100644 --- a/packages/loop_common/src/loop_common/supports/_3d_structured_grid.py +++ b/packages/loop_common/src/loop_common/supports/_3d_structured_grid.py @@ -87,80 +87,6 @@ def cell_centres(self, global_index): + self.step_vector[None, :] * 0.5 ) - def trilinear(self, local_coords): - """ - Returns the trilinear interpolation for the local coordinates - - Parameters - ---------- - x - double, array of doubles - y - double, array of doubles - z - double, array of doubles - - Returns - ------- - array of interpolation coefficients - - """ - interpolant = np.zeros((local_coords.shape[0], 8), dtype=np.float64) - # interpolant[:,0] = (1 - local_coords[:, 0])* (1 - local_coords[:, 1])* (1 - local_coords[:, 2]) - # interpolant[:,1] = local_coords[:, 0]* (1 - local_coords[:, 1])* (1 - local_coords[:, 2]) - # interpolant[:,2] = (1 - local_coords[:, 0])* local_coords[:, 1]* (1 - local_coords[:, 2]) - # interpolant[:,3] = (1 - local_coords[:, 0])* (1 - local_coords[:, 1])* local_coords[:, 2] - # interpolant[:,4] = local_coords[:, 0] * (1 - local_coords[:, 1]) * local_coords[:, 2] - # interpolant[:,5] = (1 - local_coords[:, 0]) * local_coords[:, 1] * local_coords[:, 2] - # interpolant[:,6] = local_coords[:, 0] * local_coords[:, 1] * (1 - local_coords[:, 2]) - # interpolant[:,7] = local_coords[:, 0] * local_coords[:, 1] * local_coords[:, 2] - interpolant[:, 0] = ( - (1 - local_coords[:, 0]) * (1 - local_coords[:, 1]) * (1 - local_coords[:, 2]) - ) - interpolant[:, 1] = local_coords[:, 0] * (1 - local_coords[:, 1]) * (1 - local_coords[:, 2]) - interpolant[:, 2] = (1 - local_coords[:, 0]) * local_coords[:, 1] * (1 - local_coords[:, 2]) - interpolant[:, 4] = (1 - local_coords[:, 0]) * (1 - local_coords[:, 1]) * local_coords[:, 2] - interpolant[:, 5] = local_coords[:, 0] * (1 - local_coords[:, 1]) * local_coords[:, 2] - interpolant[:, 6] = (1 - local_coords[:, 0]) * local_coords[:, 1] * local_coords[:, 2] - interpolant[:, 3] = local_coords[:, 0] * local_coords[:, 1] * (1 - local_coords[:, 2]) - interpolant[:, 7] = local_coords[:, 0] * local_coords[:, 1] * local_coords[:, 2] - return interpolant - - def position_to_local_coordinates(self, pos): - """ - Convert from global to local coordinates within a cel - - Parameters - ---------- - pos - array of positions inside - - Returns - ------- - localx, localy, localz - - """ - pos = self.check_position(pos) - cell_indexes, _inside = self.position_to_cell_index(pos) - - # Use cell origin (not modulo) so points on max faces/vertices map to local=1. - cell_origin = self.origin[None, :] + self.step_vector[None, :] * cell_indexes - local_coords = (pos - cell_origin) / self.step_vector[None, :] - local_coords = np.clip(local_coords, 0.0, 1.0) - return local_coords - - def position_to_dof_coefs(self, pos): - """ - Global posotion to interpolation coefficients - - Parameters - ---------- - pos - - Returns - ------- - - """ - local_coords = self.position_to_local_coordinates(pos) - weights = self.trilinear(local_coords) - return weights - def neighbour_global_indexes(self, mask=None, **kwargs): """ Get neighbour indexes @@ -292,35 +218,6 @@ def neighbour_global_indexes(self, mask=None, **kwargs): + self.nsteps[0, None, None] * self.nsteps[1, None, None] * neighbours[2, :, :] ).astype(np.int64) - def evaluate_value(self, evaluation_points, property_array): - """ - Evaluate the value of of the property at the locations. - Trilinear interpolation dot corner values - - Parameters - ---------- - evaluation_points np array of locations - property_name string of property name - - Returns - ------- - - """ - if property_array.shape[0] != self.n_nodes: - logger.error("Property array does not match grid") - raise ValueError( - f"cannot assign {property_array.shape[0]} vlaues to array of shape {self.n_nodes}" - ) - idc, inside = self.position_to_cell_corners(evaluation_points) - # print(idc[inside,:], self.n_nodes,inside) - if idc.shape[0] != inside.shape[0]: - raise ValueError("index does not match number of nodes") - v = np.zeros(idc.shape) - v[:, :] = np.nan - v[inside, :] = self.position_to_dof_coefs(evaluation_points[inside, :]) - v[inside, :] *= property_array[idc[inside, :]] - return np.sum(v, axis=1) - def evaluate_gradient(self, evaluation_points, property_array) -> np.ndarray: """Evaluate the gradient at a location given node values diff --git a/packages/loop_common/tests/test_property_grid.py b/packages/loop_common/tests/test_property_grid.py index 797662981..03eb5ed88 100644 --- a/packages/loop_common/tests/test_property_grid.py +++ b/packages/loop_common/tests/test_property_grid.py @@ -89,3 +89,75 @@ def test_vector_valued_cell_data(grid): grid.cell_properties["v"] = np.column_stack([np.arange(grid.n_cells), -np.arange(grid.n_cells)]).astype(float) values = grid.cell_value_at("v", grid.cell_centres[[0, -1]]) np.testing.assert_allclose(values, [[0, 0], [grid.n_cells - 1, -(grid.n_cells - 1)]]) + + +# --------------------------------------------------------------------------- trilinear +LINEAR = np.array([1.5, -2.0, 0.75]) + + +def _linear(points): + return points @ LINEAR + 3.0 + + +def test_trilinear_reproduces_a_linear_field_exactly(grid): + grid.properties["f"] = _linear(grid.nodes) + points = np.random.default_rng(1).uniform(grid.origin, grid.maximum, (300, 3)) + np.testing.assert_allclose(grid.evaluate_value("f", points), _linear(points), atol=1e-9) + + +def test_trilinear_returns_the_node_values_at_the_nodes(grid): + values = np.random.default_rng(2).normal(size=grid.n_nodes) + grid.properties["f"] = values + np.testing.assert_allclose(grid.evaluate_value("f", grid.nodes), values, atol=1e-12) + + +def test_trilinear_is_nan_outside(grid): + grid.properties["f"] = _linear(grid.nodes) + outside = np.array([grid.maximum + grid.step_vector, grid.origin - grid.step_vector]) + assert np.isnan(grid.evaluate_value("f", outside)).all() + + +def test_trilinear_at_the_cell_midpoint_is_the_mean_of_its_corner_values(grid): + values = np.random.default_rng(3).normal(size=grid.n_nodes) + grid.properties["f"] = values + geometry = grid.geometry + centres = grid.cell_centres + corners, _ = geometry.position_to_cell_corners(centres) + np.testing.assert_allclose(grid.evaluate_value("f", centres), values[corners].mean(axis=1)) + + +def test_trilinear_multicomponent(grid): + grid.properties["v"] = np.column_stack([_linear(grid.nodes), -_linear(grid.nodes)]) + points = np.random.default_rng(4).uniform(grid.origin, grid.maximum, (50, 3)) + result = grid.evaluate_value("v", points) + assert result.shape == (50, 2) + np.testing.assert_allclose(result[:, 0], _linear(points), atol=1e-9) + np.testing.assert_allclose(result[:, 1], -_linear(points), atol=1e-9) + + +def test_trilinear_agrees_with_the_interpolation_support(grid): + from loop_common.supports import StructuredGrid + + support = StructuredGrid( + origin=grid.origin, nsteps_cells=grid.nsteps - 1, step_vector=grid.step_vector + ) + values = np.random.default_rng(5).normal(size=grid.n_nodes) + grid.properties["f"] = values + points = np.random.default_rng(6).uniform(grid.origin - grid.step_vector, grid.maximum + grid.step_vector, (200, 3)) + np.testing.assert_allclose( + grid.evaluate_value("f", points), support.evaluate_value(points, values), equal_nan=True + ) + + +def test_transfer_a_property_onto_a_finer_support(grid): + """The use case: evaluate a property at the nodes of another support.""" + from loop_common.supports import StructuredGrid + + grid.properties["f"] = _linear(grid.nodes) + finer = StructuredGrid( + origin=grid.origin, + nsteps_cells=(grid.nsteps - 1) * 3, + step_vector=grid.step_vector / 3.0, + ) + transferred = grid.evaluate_value("f", finer.nodes) + np.testing.assert_allclose(transferred, _linear(finer.nodes), atol=1e-9)