diff --git a/tests/create_legacy_data.py b/tests/create_legacy_data.py
index ba39985..015665b 100644
--- a/tests/create_legacy_data.py
+++ b/tests/create_legacy_data.py
@@ -28,15 +28,13 @@ def create_reference_data(
     mesh = dolfin.UnitCubeMesh(1, 1, 1)
     V = dolfin.FunctionSpace(mesh, family, degree)
     W = dolfin.VectorFunctionSpace(mesh, family, degree)
-    x = dolfin.SpatialCoordinate(mesh)
-
-    # NOTE: f0 must be exactly representable in the (family, degree) function space above.
-    # Legacy dolfin (FIAT-based quadrature) and current dolfinx (basix-based quadrature)
-    # use different quadrature rules, so a non-polynomial (e.g. conditional) integrand
-    # would project to slightly different values in each toolchain, making bit-exact
-    # comparisons in test_legacy_readers.py inherently flaky across dolfinx versions.
-    f0 = x[0] ** 2 + 2 * x[1] * x[2] - x[1]
-    v0 = dolfin.project(f0, V)
+
+    # NOTE: f0 must be exactly representable in the (family, degree) function space above,
+    # and written via interpolation (not dolfin.project, which goes through a linear solve
+    # and does not reproduce exact nodal values here) so that comparisons against a fresh
+    # dolfinx solve/interpolation in test_legacy_readers.py are not sensitive to quadrature
+    # or solver differences between legacy dolfin and current dolfinx.
+    v0 = dolfin.interpolate(dolfin.Expression("x[0]*x[0] + 2*x[1]*x[2] - x[1]", degree=2), V)
     w0 = dolfin.interpolate(dolfin.Expression(("x[0]", "3*x[2]", "7*x[1]"), degree=1), W)
 
     v1 = dolfin.interpolate(dolfin.Expression("x[0]", degree=1), V)