--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/experiments/tests/test_minimisation.py	Thu Feb 26 09:32:12 2026 -0500
@@ -0,0 +1,34 @@
+from dolfinx_access import eval_Function_f64, min_Function_f64_p2
+
+
+# A simple test that the finite-dimensional subproblem is solved correctly
+def test_minimisation(pde):
+    test = fem.Function(pde.V)
+
+    # Beale function
+    def beale(x):
+        return (
+            (1.5 - x[0] + x[0] * x[1]) ** 2
+            + (2.25 - x[0] + x[0] * (x[1] ** 2)) ** 2
+            + (2.625 - x[0] + x[0] * (x[1] ** 3)) ** 2
+        )
+
+    test.interpolate(beale)
+    m = min_Function_f64_p2(test._cpp_object)
+    print("Beale minimizer:", m.x, m.v)
+
+    def quadratic_testf(x):
+        return x[0] ** 2 - 5 * x[1] ** 2 + x[0] * x[1] - x[0] + x[1]
+
+    testf = quadratic_testf
+    test.interpolate(testf)
+    n = 256
+    dist = 0.0
+    for x in np.linspace(pde.domain_size[0], pde.domain_size[1], n):
+        for y in np.linspace(pde.domain_size[2], pde.domain_size[3], n):
+            p = [x, y, 0.0]
+            v = eval_Function_f64(test._cpp_object, p)
+            v0 = testf(p)
+            d = np.linalg.norm(v - v0)
+            dist = max(d, dist)
+    print("Quadratic function evaluation max discrepancy:", dist)