From 2e1b5204b9088a7d802ece2cd0c80806991b15b7 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Tue, 20 Oct 2026 10:27:12 +1030 Subject: [PATCH 1/5] feat(interp): add centroid and gauss4 tet quadrature rules --- .../src/loop_interpolation/_quadrature.py | 75 +++++++++++++++++++ .../tests/test_quadrature.py | 58 ++++++++++++++ 2 files changed, 133 insertions(+) create mode 100644 packages/loop_interpolation/src/loop_interpolation/_quadrature.py create mode 100644 packages/loop_interpolation/tests/test_quadrature.py diff --git a/packages/loop_interpolation/src/loop_interpolation/_quadrature.py b/packages/loop_interpolation/src/loop_interpolation/_quadrature.py new file mode 100644 index 00000000..48038613 --- /dev/null +++ b/packages/loop_interpolation/src/loop_interpolation/_quadrature.py @@ -0,0 +1,75 @@ +"""Tetrahedron volume quadrature rules. + +Barycentric convention: c0 = 1 - r - s - t, c1..c3 = r, s, t. +Weights are normalised to sum to 1, so a rule gives the mean of a +function over the tetrahedron. Multiply by the element volume for +the integral. +""" + +import numpy as np + +_GAUSS4_A = 0.1381966011250105 +_GAUSS4_B = 0.5854101966249685 + +_RULES = ("centroid", "gauss4") + + +def tet_quadrature(rule): + """Return a tetrahedron quadrature rule. + + Parameters + ---------- + rule : str + 'centroid' (1 point, exact for linear functions) or + 'gauss4' (4 points, exact for quadratic functions). + + Returns + ------- + bary : np.ndarray + (n, 4) barycentric coordinates of the points. + weights : np.ndarray + (n,) weights that sum to 1. + + Raises + ------ + ValueError + If the rule name is not known. + """ + if rule == "centroid": + return np.full((1, 4), 0.25), np.ones(1) + if rule == "gauss4": + a, b = _GAUSS4_A, _GAUSS4_B + bary = np.full((4, 4), a) + bary[np.arange(4), np.arange(4)] = b + return bary, np.full(4, 0.25) + raise ValueError(f"Unknown quadrature rule '{rule}', use one of {_RULES}") + + +def tet_quadrature_points(nodes, elements, rule): + """Map a quadrature rule onto tetrahedra. + + Parameters + ---------- + nodes : np.ndarray + (n_nodes, 3) node coordinates. + elements : np.ndarray + (m, >=4) element node indices. Only the 4 corners are used. + rule : str + Rule name, see `tet_quadrature`. + + Returns + ------- + points : np.ndarray + (m * n, 3) points, grouped by element. + tet_idx : np.ndarray + (m * n,) index of the element (row of `elements`) for each point. + weights : np.ndarray + (m * n,) rule weights repeated for each element (sum to 1 per element). + """ + bary, w = tet_quadrature(rule) + corners = np.asarray(nodes)[np.asarray(elements)[:, :4]] # (m, 4, 3) + points = np.einsum("pk,mkd->mpd", bary, corners).reshape(-1, 3) + n = len(w) + tet_idx = np.repeat(np.arange(corners.shape[0]), n) + weights = np.tile(w, corners.shape[0]) + return points, tet_idx, weights diff --git a/packages/loop_interpolation/tests/test_quadrature.py b/packages/loop_interpolation/tests/test_quadrature.py new file mode 100644 index 00000000..29d6b30f --- /dev/null +++ b/packages/loop_interpolation/tests/test_quadrature.py @@ -0,0 +1,58 @@ +import itertools +from math import factorial + +import numpy as np +import pytest + +from loop_interpolation._quadrature import tet_quadrature, tet_quadrature_points + +RULES = ["centroid", "gauss4"] + + +def _exact_mean(a): + # mean of r^a0 s^a1 t^a2 over the reference tet: 3! * prod(a_i!) / (3 + sum a)! + return 6.0 * np.prod([factorial(i) for i in a]) / factorial(3 + sum(a)) + + +@pytest.mark.parametrize("rule", RULES) +def test_weights_sum_to_one(rule): + bary, w = tet_quadrature(rule) + assert np.isclose(w.sum(), 1.0) + assert np.allclose(bary.sum(axis=1), 1.0) + + +@pytest.mark.parametrize( + "rule,max_degree", [("centroid", 1), ("gauss4", 2)] +) +def test_exact_for_polynomials(rule, max_degree): + bary, w = tet_quadrature(rule) + rst = bary[:, 1:] + for a in itertools.product(range(max_degree + 1), repeat=3): + if sum(a) > max_degree: + continue + f = np.prod(rst ** np.array(a), axis=1) + assert np.isclose(w @ f, _exact_mean(a)), a + + +def test_unknown_rule(): + with pytest.raises(ValueError): + tet_quadrature("nope") + with pytest.raises(ValueError): + tet_quadrature_points(np.zeros((4, 3)), np.arange(4)[None], "nope") + + +def test_mapping_indices_and_points(): + nodes = np.array( + [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 1]], dtype=float + ) + elements = np.array([[0, 1, 2, 3], [1, 2, 3, 4]]) + pts, idx, w = tet_quadrature_points(nodes, elements, "gauss4") + assert pts.shape == (8, 3) + assert np.array_equal(idx, [0, 0, 0, 0, 1, 1, 1, 1]) + assert np.allclose(w.reshape(2, 4).sum(axis=1), 1.0) + for e in range(2): + mean = (w[idx == e, None] * pts[idx == e]).sum(axis=0) + assert np.allclose(mean, nodes[elements[e]].mean(axis=0)) + pc, ic, wc = tet_quadrature_points(nodes, elements, "centroid") + assert np.array_equal(ic, [0, 1]) + assert np.allclose(pc, nodes[elements].mean(axis=1)) From ec7752149a33cec56e2f2882879bb3e0ea4795d6 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Tue, 20 Oct 2026 10:35:22 +1030 Subject: [PATCH 2/5] fix(supports): P2 evaluate_shape_derivatives returns NaN for outside points Compute only for valid elements, NaN rows and element -1 elsewhere. Also return (dN, elements, inside) to match P1 (review item 17); update callers and fakes. --- .../supports/_2d_p2_unstructured.py | 31 +++++++++++++-- .../src/loop_common/supports/_3d_p2_tetra.py | 29 ++++++++++++-- .../supports/_p2_structured_tetra.py | 27 ++++++++++++- .../test_p2_shape_derivatives_outside.py | 39 +++++++++++++++++++ .../src/loop_interpolation/_p2interpolator.py | 19 +++++---- .../unit/interpolator/test_p2interpolator.py | 2 +- 6 files changed, 129 insertions(+), 18 deletions(-) create mode 100644 packages/loop_common/tests/test_p2_shape_derivatives_outside.py diff --git a/packages/loop_common/src/loop_common/supports/_2d_p2_unstructured.py b/packages/loop_common/src/loop_common/supports/_2d_p2_unstructured.py index 09e26d0f..fcd59450 100644 --- a/packages/loop_common/src/loop_common/supports/_2d_p2_unstructured.py +++ b/packages/loop_common/src/loop_common/supports/_2d_p2_unstructured.py @@ -211,7 +211,31 @@ def evaluate_shape_d2(self, indexes): ii += 1 return d2 - def evaluate_shape_derivatives(self, locations, elements=None): + def evaluate_shape_derivatives(self, locations: np.ndarray, elements=None): + """Evaluate P2 shape derivatives in global coordinates. + + Returns + ------- + dN, elements, inside + dN has shape (n, 2, n_nodes_per_element). Points outside the mesh + have element index -1, ``inside`` False and NaN rows in dN. They do + not affect other points in the batch. + """ + locations = np.asarray(locations) + if elements is None: + _verts, _c, elements, inside = self.get_element_for_location(locations) + elements = np.asarray(elements) + valid = (elements >= 0) & np.asarray(inside, dtype=bool) + else: + elements = np.asarray(elements) + valid = elements >= 0 + out_elements = np.where(valid, elements, -1) + dN = np.full((locations.shape[0], 2, self.elements.shape[1]), np.nan) + if np.any(valid): + dN[valid], _ = self._evaluate_shape_derivatives_valid(locations[valid], elements[valid]) + return dN, out_elements, valid + + def _evaluate_shape_derivatives_valid(self, locations, elements=None): """ Compute dN/ds (1st row), dN/dt(2nd row) """ @@ -221,7 +245,7 @@ def evaluate_shape_derivatives(self, locations, elements=None): else: tri = elements M = np.ones((elements.shape[0], 3, 3)) - M[:, :, 1:] = self.vertices[self.elements[elements], :][:, :3, :] + M[:, :, 1:] = self.nodes[self.elements[elements], :][:, :3, :] points_ = np.ones((locations.shape[0], 3)) points_[:, 1:] = locations minv = np.linalg.inv(M) @@ -334,8 +358,7 @@ def evaluate_gradient(self, pos: np.ndarray, property_array: np.ndarray) -> np.n raise ValueError("property array must have same length as nodes") values = np.zeros(pos.shape) values[:] = np.nan - element_gradients, tri = self.evaluate_shape_derivatives(pos[:, :2]) - inside = tri >= 0 + element_gradients, tri, inside = self.evaluate_shape_derivatives(pos[:, :2]) values[inside, :] = ( element_gradients[inside, :, :] * property_array[self.elements[tri[inside], None, :]] ).sum(2) diff --git a/packages/loop_common/src/loop_common/supports/_3d_p2_tetra.py b/packages/loop_common/src/loop_common/supports/_3d_p2_tetra.py index b73cf7d2..43e36fe0 100644 --- a/packages/loop_common/src/loop_common/supports/_3d_p2_tetra.py +++ b/packages/loop_common/src/loop_common/supports/_3d_p2_tetra.py @@ -191,7 +191,31 @@ def evaluate_shape_d2(self, indexes: np.ndarray) -> np.ndarray: ii += 1 return d2 - def evaluate_shape_derivatives( + def evaluate_shape_derivatives(self, locations: np.ndarray, elements=None): + """Evaluate P2 shape derivatives in global coordinates. + + Returns + ------- + dN, elements, inside + dN has shape (n, 3, n_nodes_per_element). Points outside the mesh + have element index -1, ``inside`` False and NaN rows in dN. They do + not affect other points in the batch. + """ + locations = np.asarray(locations) + if elements is None: + _verts, _c, elements, inside = self.get_element_for_location(locations) + elements = np.asarray(elements) + valid = (elements >= 0) & np.asarray(inside, dtype=bool) + else: + elements = np.asarray(elements) + valid = elements >= 0 + out_elements = np.where(valid, elements, -1) + dN = np.full((locations.shape[0], 3, self.elements.shape[1]), np.nan) + if np.any(valid): + dN[valid], _ = self._evaluate_shape_derivatives_valid(locations[valid], elements[valid]) + return dN, out_elements, valid + + def _evaluate_shape_derivatives_valid( self, locations: np.ndarray, elements: np.ndarray = None ) -> np.ndarray: """ @@ -372,8 +396,7 @@ def evaluate_gradient(self, pos: np.ndarray, property_array: np.ndarray) -> np.n raise ValueError("property array must have same length as nodes") values = np.zeros(pos.shape) values[:] = np.nan - element_gradients, tetra = self.evaluate_shape_derivatives(pos[:, :3]) - inside = tetra >= 0 + element_gradients, tetra, inside = self.evaluate_shape_derivatives(pos[:, :3]) values[inside, :] = ( element_gradients[:, :, :] * property_array[self.elements[tetra[inside], None, :]] ).sum(2) diff --git a/packages/loop_common/src/loop_common/supports/_p2_structured_tetra.py b/packages/loop_common/src/loop_common/supports/_p2_structured_tetra.py index 04e8ac60..f5cf26c2 100644 --- a/packages/loop_common/src/loop_common/supports/_p2_structured_tetra.py +++ b/packages/loop_common/src/loop_common/supports/_p2_structured_tetra.py @@ -260,6 +260,30 @@ def evaluate_shape(self, locations: np.ndarray): return N, elements, inside def evaluate_shape_derivatives(self, locations: np.ndarray, elements=None): + """Evaluate P2 shape derivatives in global coordinates. + + Returns + ------- + dN, elements, inside + dN has shape (n, 3, n_nodes_per_element). Points outside the mesh + have element index -1, ``inside`` False and NaN rows in dN. They do + not affect other points in the batch. + """ + locations = np.asarray(locations) + if elements is None: + _verts, _c, elements, inside = self.get_element_for_location(locations) + elements = np.asarray(elements) + valid = (elements >= 0) & np.asarray(inside, dtype=bool) + else: + elements = np.asarray(elements) + valid = elements >= 0 + out_elements = np.where(valid, elements, -1) + dN = np.full((locations.shape[0], 3, self.elements.shape[1]), np.nan) + if np.any(valid): + dN[valid], _ = self._evaluate_shape_derivatives_valid(locations[valid], elements[valid]) + return dN, out_elements, valid + + def _evaluate_shape_derivatives_valid(self, locations: np.ndarray, elements=None): """Evaluate quadratic tetrahedral shape derivatives at locations.""" locations = np.array(locations) if elements is None: @@ -480,8 +504,7 @@ def evaluate_gradient(self, pos: np.ndarray, property_array: np.ndarray) -> np.n values = np.zeros(pos.shape) values[:] = np.nan - element_gradients, tetras = self.evaluate_shape_derivatives(pos) - _, _, inside = self.evaluate_shape(pos) + element_gradients, tetras, inside = self.evaluate_shape_derivatives(pos) values[inside, :] = np.einsum( "ijk,ik->ij", element_gradients[inside, :, :], diff --git a/packages/loop_common/tests/test_p2_shape_derivatives_outside.py b/packages/loop_common/tests/test_p2_shape_derivatives_outside.py new file mode 100644 index 00000000..b030d091 --- /dev/null +++ b/packages/loop_common/tests/test_p2_shape_derivatives_outside.py @@ -0,0 +1,39 @@ +"""evaluate_shape_derivatives must handle points outside the mesh.""" + +import numpy as np +import pytest +from loop_common.supports._3d_p2_tetra import P2UnstructuredTetMesh +from loop_common.supports._p2_structured_tetra import P2TetMesh + + +def _unstructured(): + return P2UnstructuredTetMesh( + origin=np.zeros(3), step_vector=np.ones(3) / 4, nsteps_cells=np.array([4, 4, 4]) + ) + + +def _structured(): + return P2TetMesh(origin=np.zeros(3), nsteps=np.array([3, 3, 3]), step_vector=np.ones(3)) + + +@pytest.mark.parametrize( + "make,inside_pt,outside_pt", + [ + (_unstructured, [0.41, 0.52, 0.33], [5.0, 5.0, 5.0]), + (_structured, [1.1, 1.2, 1.3], [50.0, 50.0, 50.0]), + ], +) +def test_outside_point_gives_nan_and_does_not_poison_batch(make, inside_pt, outside_pt): + mesh = make() + ref = mesh.evaluate_shape_derivatives(np.array([inside_pt])) + pts = np.array([inside_pt, outside_pt, inside_pt]) + res = mesh.evaluate_shape_derivatives(pts) + assert len(res) == 3 + dN, elements, inside = res + assert dN.shape == (3, 3, 10) + assert elements[1] == -1 + assert not inside[1] + assert np.all(np.isnan(dN[1])) + assert inside[0] and inside[2] + assert np.all(np.isfinite(dN[0])) and np.all(np.isfinite(dN[2])) + np.testing.assert_allclose(dN[0], ref[0][0]) diff --git a/packages/loop_interpolation/src/loop_interpolation/_p2interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_p2interpolator.py index a3b765de..608b83cf 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_p2interpolator.py +++ b/packages/loop_interpolation/src/loop_interpolation/_p2interpolator.py @@ -109,8 +109,9 @@ def copy(self): def add_gradient_constraints(self, w: float = 1.0): points = self.get_gradient_constraints() if points.shape[0] > 0: - grad, elements = self.support.evaluate_shape_derivatives(points[:, : self.dimensions]) - inside = elements > -1 + grad, elements, inside = self.support.evaluate_shape_derivatives( + points[:, : self.dimensions] + ) # element size cancels in the row normalisation, w is applied after it wt = np.full(int(np.sum(inside)), float(w)) A = np.einsum( @@ -140,8 +141,9 @@ def add_gradient_orthogonal_constraints( """ if points.shape[0] > 0: - grad, elements = self.support.evaluate_shape_derivatives(points[:, : self.dimensions]) - inside = elements > -1 + grad, elements, inside = self.support.evaluate_shape_derivatives( + points[:, : self.dimensions] + ) wt = np.full(int(np.sum(inside)), float(w)) A = np.einsum("ijk,ij->ik", grad[inside, :], vector[inside, :]) B = np.zeros(A.shape[0]) @@ -153,8 +155,9 @@ def add_gradient_orthogonal_constraints( def add_norm_constraints(self, w: float = 1.0): points = self.get_norm_constraints() if points.shape[0] > 0: - grad, elements = self.support.evaluate_shape_derivatives(points[:, : self.dimensions]) - inside = elements > -1 + grad, elements, inside = self.support.evaluate_shape_derivatives( + points[:, : self.dimensions] + ) wt = np.full(int(np.sum(inside)), float(w)) elements = np.tile(self.support.elements[elements[inside]], (self.dimensions, 1, 1)) elements = elements.swapaxes(0, 1) @@ -255,10 +258,10 @@ def minimise_edge_jumps( if callable(vector_func): norm = vector_func(cp[:, i, :]) # evaluate the shape function for the edges for each neighbouring triangle - cp_Dt, cp_tri1 = self.support.evaluate_shape_derivatives( + cp_Dt, cp_tri1, _inside1 = self.support.evaluate_shape_derivatives( cp[:, i, :], elements=self.support.shared_element_relationships[:, 0] ) - cp_Dn, cp_tri2 = self.support.evaluate_shape_derivatives( + cp_Dn, cp_tri2, _inside2 = self.support.evaluate_shape_derivatives( cp[:, i, :], elements=self.support.shared_element_relationships[:, 1] ) # constraint for each cp is triangle - neighbour create a Nx12 matrix diff --git a/tests/unit/interpolator/test_p2interpolator.py b/tests/unit/interpolator/test_p2interpolator.py index b70401fe..27a8921a 100644 --- a/tests/unit/interpolator/test_p2interpolator.py +++ b/tests/unit/interpolator/test_p2interpolator.py @@ -23,7 +23,7 @@ def evaluate_shape_derivatives(self, points): n = points.shape[0] grad = np.ones((n, self.elements.shape[1], 3)) elements = np.zeros(n, dtype=int) - return grad, elements + return grad, elements, np.ones(n, dtype=bool) def evaluate_shape(self, points): n = points.shape[0] From 1ba9d39bb14866c2367b89d33f85c26931a54694 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Tue, 20 Oct 2026 10:35:27 +1030 Subject: [PATCH 3/5] refactor(supports): update P2 structured test for 3-value evaluate_shape_derivatives --- packages/loop_common/tests/test_p2_structured_tetra.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/loop_common/tests/test_p2_structured_tetra.py b/packages/loop_common/tests/test_p2_structured_tetra.py index e855906c..940e2761 100644 --- a/packages/loop_common/tests/test_p2_structured_tetra.py +++ b/packages/loop_common/tests/test_p2_structured_tetra.py @@ -183,7 +183,7 @@ def test_p2tetmesh_shape_function_derivatives_exist(self): mesh = P2TetMesh(nsteps=np.array([3, 3, 3])) centroid = np.array([[1.0, 1.0, 1.0]]) - dN, _elem_ids = mesh.evaluate_shape_derivatives(centroid) + dN, _elem_ids, _inside = mesh.evaluate_shape_derivatives(centroid) # dN should have shape (n_points, 3, 10) assert dN.shape == (1, 3, 10) From e471212f277941706712d60f723baa0419cb4698 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Tue, 20 Oct 2026 10:37:47 +1030 Subject: [PATCH 4/5] feat(interp): add P2FoldInterpolator with quadrature-sampled fold constraints --- .../_p2_fold_interpolator.py | 129 ++++++++++++++++++ .../tests/test_p2_fold_interpolator.py | 77 +++++++++++ 2 files changed, 206 insertions(+) create mode 100644 packages/loop_interpolation/src/loop_interpolation/_p2_fold_interpolator.py create mode 100644 packages/loop_interpolation/tests/test_p2_fold_interpolator.py diff --git a/packages/loop_interpolation/src/loop_interpolation/_p2_fold_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_p2_fold_interpolator.py new file mode 100644 index 00000000..503688ea --- /dev/null +++ b/packages/loop_interpolation/src/loop_interpolation/_p2_fold_interpolator.py @@ -0,0 +1,129 @@ +""" +Piecewise quadratic interpolator using fold constraints. + +Summary of findings_1_1.md (what exists / what was added) +--------------------------------------------------------- +Existing: +- ``P2Interpolator`` and the P2 supports. ``evaluate_shape_derivatives`` + gives dN of shape (n, 3, 10), node order [v0..v3, e23, e03, e01, e12, e13, e02]. +- ``setup_with_fold_constraints`` (_fold_setup.py) and + ``resolve_fold_norm_target`` (_fold_norm_alignment.py), both reused as is. +- ``add_constraints_to_least_squares`` drops NaN rows, normalises the rows of A, + then applies the weight ``w``. A per-row ``w`` is therefore applied after + normalisation. +- ``get_quadrature_points`` on the support is a triangle (face) rule for edge + jumps. It does not fit a volume rule. +Added (Phase 1): +- ``_quadrature.py``: ``tet_quadrature`` / ``tet_quadrature_points`` + (centroid, gauss4, weights sum to 1, tet index kept). +- P2 ``evaluate_shape_derivatives`` returns (dN, elements, inside) with NaN + rows for outside points. +- This class. ``copy()`` is overridden because the base returns P2Interpolator. +""" +from __future__ import annotations + +from typing import Callable + +import numpy as np +from loop_common.logging import get_logger as getLogger + +from ._fold_event import FoldEvent +from ._fold_norm_alignment import resolve_fold_norm_target +from ._fold_setup import setup_with_fold_constraints +from ._p2interpolator import P2Interpolator +from ._quadrature import tet_quadrature_points + +logger = getLogger(__name__) + +_SAMPLE_RULES = ("centroid", "gauss4") + + +class P2FoldInterpolator(P2Interpolator): + """P2 interpolator with fold constraints (Laurent et al., 2016).""" + + def __init__(self, support, fold: FoldEvent | None = None): + P2Interpolator.__init__(self, support) + self.fold = fold + + def copy(self): + return P2FoldInterpolator(self.support, fold=self.fold) + + def setup_interpolator(self, **kwargs): + report = setup_with_fold_constraints( + fold=self.fold, + kwargs=kwargs, + interpolator_name=self.__class__.__name__, + base_setup=super().setup_interpolator, + add_fold_constraints=self.add_fold_constraints, + finalize_report=self.finalize_setup_diagnostics_report, + ) + # set after the call: includes fold_weights, which the base setup does not see + self._setup_kwargs = dict(kwargs) + return report + + def add_fold_constraints( + self, + fold_orientation=10.0, + fold_axis_w=10.0, + fold_normalisation=1.0, + fold_norm=-1.0, + dgz_alignment="warn", + sample_rule="gauss4", + mask_fn: Callable | None = None, + ): + """Add the three fold constraint families at quadrature points. + + Parameters + ---------- + fold_orientation, fold_axis_w, fold_normalisation : float or None + weights of the families. None skips a family. + fold_norm : float + norm target, see ``resolve_fold_norm_target`` + dgz_alignment : {"none", "warn", "correct"} + sample_rule : {"gauss4", "centroid"} + quadrature rule used to place sample points in every tetrahedron + mask_fn : callable, optional + called with the sample points; points where it returns True get zero weight + """ + if sample_rule not in _SAMPLE_RULES: + raise ValueError(f"Unknown sample_rule '{sample_rule}', use one of {_SAMPLE_RULES}") + support = self.support + points, tet_idx, q = tet_quadrature_points(support.nodes, support.elements, sample_rule) + dN, _out_elements, inside = support.evaluate_shape_derivatives(points, elements=tet_idx) + deformed_orientation, fold_axis, dgz = self.fold.get_deformed_orientation(points) + # energy per tet must not depend on the rule: weight the residual by w*sqrt(q) + sq = np.sqrt(q) + if mask_fn is not None: + sq = sq * (~np.asarray(mask_fn(points), dtype=bool)) + # keep finite rows only, so A, B, idc and w stay aligned + valid = np.asarray(inside, dtype=bool) & np.all(np.isfinite(dN), axis=(1, 2)) + idc_all = support.elements[tet_idx] + + def _add(vectors, target, w, name): + v = valid & np.all(np.isfinite(vectors), axis=1) + # no volume factor: row normalisation would cancel it + A = np.einsum("ij,ijk->ik", vectors[v], dN[v]) + B = np.full(A.shape[0], float(target)) + self.add_constraints_to_least_squares(A, B, idc_all[v], w=float(w) * sq[v], name=name) + + if fold_orientation is not None: + logger.info(f"Adding fold orientation constraint to w = {fold_orientation}") + _add(deformed_orientation, 0.0, fold_orientation, "fold orientation") + if fold_axis_w is not None: + logger.info(f"Adding fold axis constraint to w = {fold_axis_w}") + _add(fold_axis, 0.0, fold_axis_w, "fold axis") + if fold_normalisation is not None: + logger.info(f"Adding fold normalisation constraint to w = {fold_normalisation}") + target_norm = resolve_fold_norm_target( + fold=self.fold, + normal_constraints=self.get_norm_constraints(), + fold_norm=fold_norm, + dgz_alignment=dgz_alignment, + logger=logger, + default_fold_norm=-1.0, + ) + # Decision: unlike DiscreteFoldInterpolator, B is NOT multiplied by + # fold_normalisation (nor by volume). That scaling changes the target + # norm with the weight, which is a quirk. Here the weight only sets + # the importance and the target stays fold_norm. + _add(dgz, target_norm, fold_normalisation, "fold normalisation") diff --git a/packages/loop_interpolation/tests/test_p2_fold_interpolator.py b/packages/loop_interpolation/tests/test_p2_fold_interpolator.py new file mode 100644 index 00000000..ae55fc35 --- /dev/null +++ b/packages/loop_interpolation/tests/test_p2_fold_interpolator.py @@ -0,0 +1,77 @@ +"""Tests for P2FoldInterpolator.""" +import numpy as np +import pytest +from loop_interpolation import P2Interpolator, P2UnstructuredTetMesh +from loop_interpolation._p2_fold_interpolator import P2FoldInterpolator +from loop_interpolation._quadrature import tet_quadrature_points + + +class SimpleFold: + """Fold with constant orientation (z), axis (y) and dgz (x).""" + + def get_deformed_orientation(self, points): + n = points.shape[0] + ori = np.tile([0.0, 0.0, 1.0], (n, 1)) + axis = np.tile([0.0, 1.0, 0.0], (n, 1)) + dgz = np.tile([1.0, 0.0, 0.0], (n, 1)) + return ori, axis, dgz + + +@pytest.fixture +def mesh(): + return P2UnstructuredTetMesh( + origin=np.zeros(3), step_vector=np.ones(3), nsteps_cells=np.array([2, 2, 2]) + ) + + +def _interp(mesh): + return P2FoldInterpolator(mesh, fold=SimpleFold()) + + +def _setup(mesh, **fw): + i = _interp(mesh) + i.set_value_constraints(np.array([[0.5, 0.5, 0.5, 0.0], [1.5, 1.5, 1.5, 1.0]])) + i.setup_interpolator(fold_weights=fw) + return i + + +@pytest.mark.parametrize("rule,n", [("gauss4", 4), ("centroid", 1)]) +def test_counts_and_weights(mesh, rule, n): + i = _setup(mesh, sample_rule=rule) + names = {"fold orientation": 10.0, "fold axis": 10.0, "fold normalisation": 1.0} + _, _, q = tet_quadrature_points(mesh.nodes, mesh.elements, rule) + for name, w in names.items(): + c = i.constraints[name] + assert c["matrix"].shape[0] == mesh.elements.shape[0] * n + np.testing.assert_allclose(c["w"], w * np.sqrt(q)) + + +def test_rows_have_no_volume_scaling(mesh): + i = _setup(mesh, sample_rule="centroid") + c = i.constraints["fold normalisation"] + A = c["matrix"].toarray() + np.testing.assert_allclose(np.linalg.norm(A, axis=1), 1.0) + pts, tet, _ = tet_quadrature_points(mesh.nodes, mesh.elements, "centroid") + dN = mesh.evaluate_shape_derivatives(pts, elements=tet)[0] + raw = dN[:, 0, :] # dgz = x + # B = target / |row|, with no volume factor + np.testing.assert_allclose(c["b"], -1.0 / np.linalg.norm(raw, axis=1)) + + +def test_invalid_rule_raises(mesh): + with pytest.raises(ValueError): + _interp(mesh).setup_interpolator(fold_weights={"sample_rule": "bad"}) + + +def test_copy_type(mesh): + c = _interp(mesh).copy() + assert type(c) is P2FoldInterpolator + assert c.fold is not None + assert isinstance(c, P2Interpolator) + + +def test_update_keeps_fold_weights(mesh): + i = _setup(mesh, sample_rule="centroid", fold_axis_w=3.0) + assert i._setup_kwargs["fold_weights"]["fold_axis_w"] == 3.0 + i.setup_interpolator(**i._setup_kwargs) + np.testing.assert_allclose(i.constraints["fold axis"]["w"], 3.0) From fe5321c41ae692e52dc49a77e046546ef0e29e6c Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Tue, 20 Oct 2026 10:40:03 +1030 Subject: [PATCH 5/5] fix(interp): repeat per-row weights for 3d constraint input; test P2 fold settings and weight alignment --- .../_discrete_interpolator.py | 3 +- .../tests/test_p2_fold_interpolator.py | 103 ++++++++++++++++++ 2 files changed, 105 insertions(+), 1 deletion(-) diff --git a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py index 834a4e40..77899ba5 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py +++ b/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py @@ -364,7 +364,8 @@ def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"): if len(A.shape) > 2: n_rows = A.shape[0] * A.shape[1] if isinstance(w, np.ndarray): - w = np.tile(w, (A.shape[1])) + # reshape below puts row (i, j) at i * A.shape[1] + j, so repeat, not tile + w = np.repeat(w, A.shape[1]) A = A.reshape((A.shape[0] * A.shape[1], A.shape[2])) idc = idc.reshape((idc.shape[0] * idc.shape[1], idc.shape[2])) B = B.reshape(A.shape[0]) diff --git a/packages/loop_interpolation/tests/test_p2_fold_interpolator.py b/packages/loop_interpolation/tests/test_p2_fold_interpolator.py index ae55fc35..65f504bb 100644 --- a/packages/loop_interpolation/tests/test_p2_fold_interpolator.py +++ b/packages/loop_interpolation/tests/test_p2_fold_interpolator.py @@ -75,3 +75,106 @@ def test_update_keeps_fold_weights(mesh): assert i._setup_kwargs["fold_weights"]["fold_axis_w"] == 3.0 i.setup_interpolator(**i._setup_kwargs) np.testing.assert_allclose(i.constraints["fold axis"]["w"], 3.0) + + +def test_update_replays_setup_kwargs(mesh): + i = _interp(mesh) + i.set_value_constraints(np.array([[0.5, 0.5, 0.5, 0.0], [1.5, 1.5, 1.5, 1.0]])) + fw = {"sample_rule": "centroid", "fold_orientation": 7.0, "fold_axis_w": 3.0} + i.setup_interpolator(fold_weights=fw, regularisation=0.3) + assert i.solve_system("cg") + i.up_to_date = False + # drop the constraints, then update() must rebuild them from the stored kwargs + i.reset() + assert i.update() + for name, w in (("fold orientation", 7.0), ("fold axis", 3.0)): + c = i.constraints[name] + assert c["matrix"].shape[0] == mesh.elements.shape[0] # centroid: 1 per tet + np.testing.assert_allclose(c["w"], w * np.sqrt(0.0 + _q_centroid(mesh))) + assert i._setup_kwargs["regularisation"] == 0.3 + assert i._setup_kwargs["fold_weights"] == fw + assert i.interpolation_weights["cgw"] == i.resolve_regularisation_config( + regularisation=0.3, directional_regularisation=None + ).isotropic + + +def _q_centroid(mesh): + return tet_quadrature_points(mesh.nodes, mesh.elements, "centroid")[2] + + +def test_copy_is_independent_fold_interpolator(mesh): + i = _setup(mesh, sample_rule="centroid") + c = i.copy() + assert c is not i and c.fold is i.fold and c.support is i.support + assert c.constraints is not i.constraints + # copy starts fresh: it can be set up with its own kwargs + c.set_value_constraints(np.array([[0.5, 0.5, 0.5, 0.0]])) + c.setup_interpolator(fold_weights={"sample_rule": "gauss4"}) + assert c.constraints["fold axis"]["matrix"].shape[0] == 4 * mesh.elements.shape[0] + assert i.constraints["fold axis"]["matrix"].shape[0] == mesh.elements.shape[0] + + +def test_weights_stay_aligned_with_rows_unequal_weights(mesh): + # a weight that differs per sample point, so a row/weight mix-up shows up + pts, tet, q = tet_quadrature_points(mesh.nodes, mesh.elements, "gauss4") + + def mask_fn(p): + return np.asarray(p)[:, 0] > 1.2 # zero weight on part of the mesh + + i = _interp(mesh) + i.set_value_constraints(np.array([[0.5, 0.5, 0.5, 0.0]])) + # a non constant dgz so rows differ and can be matched back to their sample + i.fold.get_deformed_orientation = lambda p: ( + np.c_[p[:, 1] + 1, p[:, 2] + 2, p[:, 0] + 3], + np.tile([0.0, 1.0, 0.0], (len(p), 1)), + np.c_[p[:, 0] + 1, np.ones(len(p)), p[:, 2] + 1], + ) + i.setup_interpolator( + fold_weights={"sample_rule": "gauss4", "fold_axis_w": None, "fold_normalisation": None} + ) + base = i.constraints["fold orientation"] + # rebuild with the mask, from scratch, and compare to the expected per-sample weights + i.reset() + i.add_fold_constraints(fold_axis_w=None, fold_normalisation=None, mask_fn=mask_fn) + c = i.constraints["fold orientation"] + dN = mesh.evaluate_shape_derivatives(pts, elements=tet)[0] + vec = i.fold.get_deformed_orientation(pts)[0] + raw = np.einsum("ij,ijk->ik", vec, dN) + raw /= np.linalg.norm(raw, axis=1)[:, None] + idc = mesh.elements[tet] + expected_w = 10.0 * np.sqrt(q) * (~mask_fn(pts)) + # unequal weights: zero on masked samples, in an irregular pattern + assert len(np.unique(np.round(expected_w, 12))) == 2 + A = c["matrix"].toarray() + assert A.shape[0] == len(pts) + exp_A = np.zeros_like(A) + for r in range(len(pts)): + exp_A[r, idc[r]] = raw[r] + np.testing.assert_allclose(A, exp_A, atol=1e-12) + np.testing.assert_allclose(c["w"], expected_w) + assert base["matrix"].shape == c["matrix"].shape + + +def test_weights_follow_rows_when_nan_rows_dropped(mesh): + from loop_interpolation._discrete_interpolator import DiscreteInterpolator # noqa: F401 + + i = _interp(mesh) + A = np.array([[1.0, 0.0], [np.nan, 1.0], [0.0, 2.0], [3.0, 4.0]]) + idc = np.array([[0, 1], [1, 2], [2, 3], [3, 4]]) + B = np.array([1.0, 2.0, 3.0, 4.0]) + w = np.array([10.0, 20.0, 30.0, 40.0]) + i.add_constraints_to_least_squares(A, B, idc, w=w, name="t") + c = i.constraints["t"] + np.testing.assert_allclose(c["w"], [10.0, 30.0, 40.0]) + np.testing.assert_allclose(c["b"], [1.0, 1.5, 0.8]) + M = c["matrix"].toarray() + assert M[1, 3] == 1.0 and M[2, 4] == 0.8 + + +def test_per_row_weight_with_3d_input_is_repeated_not_tiled(mesh): + # (n, k, m) input is flattened row-major: row (i, j) -> i*k + j + i = _interp(mesh) + A = np.ones((2, 3, 2)) + idc = np.tile(np.array([[0, 1]]), (2, 3, 1)) + i.add_constraints_to_least_squares(A, np.ones((2, 3)), idc, w=np.array([1.0, 2.0]), name="t3") + np.testing.assert_allclose(i.constraints["t3"]["w"], [1, 1, 1, 2, 2, 2])