[BUG] StructuredGrid2D.get_element_gradient_for_location omits the chain-rule division by step_vector

[Jimmy-KL]

Version: 1.6.28 (also present in 1.6.27), Python 3.11.10, Windows

Summary

StructuredGrid2D.get_element_gradient_for_location returns the bilinear shape-function
derivatives with respect to local (normalised 0–1) cell coordinates and never converts them to
world units. StructuredGrid3D.get_element_gradient_for_location does the conversion:

# LoopStructural/interpolators/supports/_3d_structured_grid.py:440-442
T[:, 0, :] /= self.step_vector[None, 0]
T[:, 1, :] /= self.step_vector[None, 1]
T[:, 2, :] /= self.step_vector[None, 2]

The 2D version (_2d_structured_grid.py:446-456) has no equivalent, so the returned operator is
inflated by a factor of step_vector per axis.

Impact

get_element_gradient_for_location is used in two places, so this affects both reading and solving:

1. StructuredGrid2D.evaluate_gradient returns step_vector * ∇f instead of ∇f.
2. FiniteDifferenceInterpolator.add_norm_constraints builds its least-squares rows from the
   same operator, so a norm constraint n effectively enforces ∇f = n / step_vector. The target
   therefore depends on nelements, which makes orientation data silently resolution-dependent.
   (add_gradient_constraints is unaffected, since scaling the operator does not change the zero
   set of an orthogonality constraint — but it is 3D-only, see the note below.)

For anisotropic cells the error is also directional, not just a magnitude scaling.

Reproducer

A field with a known unit gradient: f(x, y) = (x + y) / sqrt(2), so |∇f| = 1.

import numpy as np
from LoopStructural.datatypes 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]))

for nelements in (5e3, 2e4, 8e4):
    interp = InterpolatorFactory.create_interpolator(
        interpolatortype="FDI", boundingbox=bbox, nelements=nelements
    )
    pts = np.random.default_rng(0).uniform(10, 90, size=(60, 2))
    interp.set_value_constraints(
        np.column_stack([pts, (pts[:, 0] + pts[:, 1]) / SQRT2])
    )
    interp.setup_interpolator()
    interp.solve_system(solver="cg")

    p = np.array([[50.0, 50.0]])
    api = np.linalg.norm(interp.evaluate_gradient(p)[0])
    h = 1.0
    fd = np.hypot(
        (interp.evaluate_value(p + [[h, 0]]) - interp.evaluate_value(p - [[h, 0]]))[0] / (2 * h),
        (interp.evaluate_value(p + [[0, h]]) - interp.evaluate_value(p - [[0, h]]))[0] / (2 * h),
    )
    step = float(np.mean(interp.support.step_vector))
    print(f"nelements={nelements:7.0f} step={step:.3f}  api={api:.4f}  finite-diff={fd:.4f}")

Observed — evaluate_gradient tracks step_vector to four decimal places, while a finite
difference of the same value field stays at the correct 1.0:

nelements=   5000 step=1.408  api=1.4088  finite-diff=1.0003
nelements=  20000 step=0.704  api=0.7040  finite-diff=0.9995
nelements=  80000 step=0.353  api=0.3509  finite-diff=0.9932

That is api / step = 1.0006, 1.0000, 0.9940 — the reported gradient is the true gradient
multiplied by the cell size.

Expected: evaluate_gradient returns ≈1.0 at every resolution.

Suggested fix

Mirror the 3D implementation at the end of get_element_gradient_for_location
(_2d_structured_grid.py, after line 454):

T[:, 0, :] /= self.step_vector[None, 0]
T[:, 1, :] /= self.step_vector[None, 1]
return vertices, T, elements, inside

A regression test asserting |∇f| ≈ 1 for f = (x + y)/sqrt(2) across two or three nelements
values would pin this down; the existing tests/unit/interpolator/test_2d_discrete_support.py
looks like the natural home.

Related

FiniteDifferenceInterpolator.add_gradient_constraints cannot be used with 2D data. It calls
get_vectors → normal_vector_to_strike_and_dip, which indexes normal_vector[:, 2]:

IndexError: index 2 is out of bounds for axis 1 with size 2

That leaves norm constraints as the only orientation option in 2D, which is why the scaling bug
above is hard to work around.