From 3d43f5ce1b0de3830ec93b814175414820867f4e Mon Sep 17 00:00:00 2001 From: "Jean Pierre Mandujano G." Date: Mon, 3 Aug 2026 14:15:28 -0500 Subject: [PATCH] fix: apply chain rule in StructuredGrid2D element gradient StructuredGrid2D.get_element_gradient_for_location returned shape-function derivatives with respect to local (0-1) cell coordinates, so gradients were inflated by step_vector on each axis and became resolution-dependent. Mirror the division by step_vector that StructuredGrid3D already does. Fixes #289 --- .../supports/_2d_structured_grid.py | 2 + .../interpolator/test_2d_discrete_support.py | 41 +++++++++++++++++++ 2 files changed, 43 insertions(+) diff --git a/packages/loop_common/src/loop_common/supports/_2d_structured_grid.py b/packages/loop_common/src/loop_common/supports/_2d_structured_grid.py index 8f611e10..974821b6 100644 --- a/packages/loop_common/src/loop_common/supports/_2d_structured_grid.py +++ b/packages/loop_common/src/loop_common/supports/_2d_structured_grid.py @@ -460,6 +460,8 @@ def get_element_gradient_for_location( T[:, 1, 2] = 1 - local_coords[:, 0] T[:, 1, 3] = local_coords[:, 0] + T[:, 0, :] /= self.step_vector[None, 0] + T[:, 1, :] /= self.step_vector[None, 1] return vertices, T, elements, inside def get_element_for_location( diff --git a/tests/unit/interpolator/test_2d_discrete_support.py b/tests/unit/interpolator/test_2d_discrete_support.py index ac9b1cdc..042c381e 100644 --- a/tests/unit/interpolator/test_2d_discrete_support.py +++ b/tests/unit/interpolator/test_2d_discrete_support.py @@ -80,3 +80,44 @@ def test_structured_grid2d_vtk_assigns_quad_cell_types(): triangulated = vtk_grid.triangulate() assert triangulated.n_cells == grid.n_elements * 2 + + +def test_evaluate_gradient_2d_world_units(): + """ + Gradient must be returned in world units, independent of the cell size. + + Regression test for #289: get_element_gradient_for_location returned the + shape-function derivatives with respect to local (0-1) cell coordinates, + so gradients were inflated by step_vector on each axis. + """ + grid = StructuredGrid2D( + origin=np.zeros(2), nsteps=np.array([10, 10]), step_vector=np.array([2.5, 0.5]) + ) + # f(x, y) = x and f(x, y) = y have unit gradients regardless of cell size + gradient_x = np.mean(grid.evaluate_gradient(grid.barycentre, grid.nodes[:, 0]), axis=0) + gradient_y = np.mean(grid.evaluate_gradient(grid.barycentre, grid.nodes[:, 1]), axis=0) + assert np.allclose(gradient_x, np.array([1.0, 0.0])) + assert np.allclose(gradient_y, np.array([0.0, 1.0])) + + +def test_fdi_2d_gradient_resolution_independent(): + """ + The interpolated gradient magnitude must not depend on nelements (#289). + """ + from LoopStructural.geometry import BoundingBox + from LoopStructural.interpolators import InterpolatorFactory + + sqrt2 = np.sqrt(2.0) + bbox = BoundingBox(dimensions=2, origin=np.array([0.0, 0.0]), maximum=np.array([100.0, 100.0])) + points = np.random.default_rng(0).uniform(10, 90, size=(60, 2)) + for nelements in (1e3, 4e3): + interpolator = InterpolatorFactory.create_interpolator( + interpolatortype="FDI", boundingbox=bbox, nelements=nelements + ) + interpolator.set_value_constraints( + np.column_stack([points, (points[:, 0] + points[:, 1]) / sqrt2]) + ) + interpolator.setup_interpolator() + interpolator.solve_system(solver="cg") + gradient_norm = np.linalg.norm(interpolator.evaluate_gradient(np.array([[50.0, 50.0]]))[0]) + assert np.isclose(gradient_norm, 1.0, atol=0.05)