diff --git a/setup.py b/setup.py
index a2827dcc..bdd0bee0 100644
--- a/setup.py
+++ b/setup.py
@@ -261,6 +261,16 @@ def configure():
         extra_compile_args=extra_compile_args,
         **conf,
     ),
+    Extension(
+        "underworld3.function._dmfield_wrapper",
+        sources=[
+            "src/underworld3/function/_dmfield_wrapper.pyx",
+        ],
+        extra_compile_args=extra_compile_args
+        + ["-Werror=incompatible-pointer-types", "-Werror=implicit-function-declaration"],
+        define_macros=[("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")],
+        **conf,
+    ),
 ]
 
 
diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi
index 99297f9e..4c827d2e 100644
--- a/src/underworld3/cython/petsc_extras.pxi
+++ b/src/underworld3/cython/petsc_extras.pxi
@@ -115,3 +115,19 @@ cdef extern from "petsc.h" nogil:
 
     # Not wrapped at this point
     PetscErrorCode VecConcatenate(PetscInt nx, const PetscVec X[], PetscVec *, PetscIS *)
+
+# ── DMField evaluation (PR #436) ──────────────────────────────────────
+# Cython externs for PETSc DMField C API — FE-exact field, gradient, and
+# Hessian evaluation at arbitrary points. Used by _dmfield_wrapper.pyx.
+
+cdef extern from "petscdmfield.h" nogil:
+    ctypedef struct _p_DMField "struct _p_DMField":
+        pass
+    PetscErrorCode DMFieldCreateDS(PetscDM dm, PetscInt field, PetscVec lvec, _p_DMField **dmf)
+    PetscErrorCode DMFieldEvaluate(_p_DMField *dmf, PetscVec pts, PetscInt dt,
+                                    void *B, void *D, void *H)
+    PetscErrorCode DMFieldDestroy(_p_DMField **dmf)
+
+cdef extern from "petscsys.h" nogil:
+    cdef enum:
+        PETSC_REAL       # PetscDataType value for double-precision real
diff --git a/src/underworld3/function/_dmfield_evaluate.py b/src/underworld3/function/_dmfield_evaluate.py
new file mode 100644
index 00000000..093f2f98
--- /dev/null
+++ b/src/underworld3/function/_dmfield_evaluate.py
@@ -0,0 +1,130 @@
+"""
+.. currentmodule:: underworld3.function._dmfield_evaluate
+
+**Internal primitive — user code should call ``uw.function.evaluate``
+for general evaluation.**
+
+``dmfield_evaluate`` uses PETSc's ``DMFieldEvaluate`` directly, bypassing
+several contracts that ``uw.function.evaluate`` provides:
+
+- Units and non-dimensionalisation frame
+- Point-location fallback ladder (RBF rescue for unlocatable points)
+- Sympy expression composition (this takes a single MeshVariable)
+- Collective ``global_evaluate`` semantics
+
+**Sanctioned uses** (internal / expert only):
+
+- Computing exact derivatives of solved FE fields (velocity strain rates,
+  temperature gradients) without L2 projection
+- Hessians for adaptivity metrics (no other UW3 source exists)
+- One-sided element-boundary evaluation for diagnostic probes
+- Internal machinery that needs to bypass the expression stack
+
+For smooth fields NOT in the FE space, the Clement recovery path
+(``uw.function.evaluate``) is superconvergent and more accurate.
+"""
+import numpy as np
+
+from ._dmfield_wrapper import DMFieldEvaluator
+
+
+def dmfield_evaluate(var, coords, gradient=True, hessian=False):
+    """Internal primitive — FE-exact field and derivative evaluation.
+
+    Evaluates a single ``MeshVariable`` at arbitrary coordinates using
+    PETSc's ``DMFieldEvaluate`` (FE basis interpolation, no L2 projection).
+
+    **User code should call ``uw.function.evaluate``**, which provides
+    units handling, location fallback, expression composition, and
+    collective parallel evaluation. This function bypasses those
+    contracts and is intended for internal machinery and expert use.
+
+    Parameters
+    ----------
+    var : MeshVariable
+        The field to evaluate (e.g. ``velocity``, ``temperature``).
+    coords : ndarray (n_points, dim)
+        Coordinates at which to evaluate.  Accepted forms:
+        - ``mesh.X.coords`` (non-dimensional [0-1])
+        - plain numpy array (assumed non-dimensional)
+        - ``UnitAwareArray`` (auto-converted to non-dimensional)
+    gradient : bool, optional
+        Return first spatial derivatives D.  Default ``True``.
+    hessian : bool, optional
+        Return second spatial derivatives H (Hessian).  Default ``False``.
+
+    Returns
+    -------
+    B : ndarray (n_points, nc) or None
+        Field values at each point.
+    D : ndarray (n_points, nc, dim) or None
+        First derivatives.  ``D[k, j, i]`` = partial(component *j*) /
+        partial x_i at point *k*.  ``None`` when ``gradient=False``.
+    H : ndarray (n_points, nc, dim, dim) or None
+        Second derivatives.  ``H[k, j, i, l]`` = partial^2(component *j*) /
+        partial x_i partial x_l at point *k*.  ``None`` when ``hessian=False``.
+
+    Notes
+    -----
+    **Parallel** — ``DMLocatePoints`` (called internally by
+    ``DMFieldEvaluate``) is COLLECTIVE on the mesh DM's communicator.
+    **All ranks must call this function**, even with zero local points.
+    Each point in the query set must lie in the **calling rank's local
+    subdomain** — a point owned by another rank may give ``NaN`` or
+    PETSc error 62.  Results are local to each rank and must be gathered
+    explicitly if a global result is needed.
+
+    **Unlocated points** — Points that ``DMLocatePoints`` cannot place
+    in any cell may return ``NaN`` or raise PETSc error 62 (the outcome
+    depends on batching and whether the point enters the location-hash
+    candidate set but fails the interior test).  Pre-filter points with
+    ``mesh.points_in_domain()`` or catch the error.  Use
+    ``uw.function.evaluate`` for RBF-interpolated fallback values at
+    domain boundaries.
+
+    **DMField lifecycle** — A fresh DMField is created and destroyed on
+    each call.  The overhead (~0.5 us) is negligible vs. evaluation cost
+    (~50-200 us).  No cache, no manual cleanup needed.
+
+    **Mesh variable data** — ``mesh.update_lvec()`` is called internally
+    to sync the variable's data to the PETSc local vector before
+    evaluation.  No manual sync is required.
+
+    Examples
+    --------
+    >>> B, D, H = dmfield_evaluate(velocity, mesh.X.coords)
+    >>>
+    >>> # Gradient of x-velocity:  D[:, 0, :]  (component 0, all dims)
+    >>> dvx_dx = D[:, 0, 0]
+    >>> dvx_dy = D[:, 0, 1]
+    >>>
+    >>> # Divergence (trace of gradient tensor):
+    >>> div_v = D[:, 0, 0] + D[:, 1, 1]
+    """
+    # --- Normalise coordinates ------------------------------------------
+    if hasattr(coords, "magnitude"):
+        from underworld3.scaling import non_dimensionalise
+        coords_array = np.asarray(non_dimensionalise(coords), dtype=np.float64)
+    else:
+        coords_array = np.asarray(coords, dtype=np.float64)
+
+    if coords_array.ndim == 1:
+        coords_array = coords_array.reshape(-1, var.mesh.dim)
+
+    mesh = var.mesh
+
+    # --- Sync local vector (collective — all ranks must call) ----------
+    mesh.update_lvec()
+
+    # --- Create -> evaluate -> destroy (no cache) ------------------------
+    evaluator = DMFieldEvaluator()
+    try:
+        evaluator.create(mesh, var)
+        B, D, H = evaluator.evaluate(
+            coords_array, mesh, var,
+            gradient=gradient, hessian=hessian,
+        )
+    finally:
+        evaluator.destroy()
+
+    return B, D, H
diff --git a/src/underworld3/function/_dmfield_wrapper.pyx b/src/underworld3/function/_dmfield_wrapper.pyx
new file mode 100644
index 00000000..087bee32
--- /dev/null
+++ b/src/underworld3/function/_dmfield_wrapper.pyx
@@ -0,0 +1,165 @@
+# cython: language_level=3
+"""
+Cython wrapper for PETSc DMField — FE-exact field/gradient/Hessian evaluation.
+
+Uses Cython extern declarations (via ``petsc_extras.pxi``) — no ctypes,
+no hardcoded ``.dylib`` paths. Portable across macOS, Linux, and CI.
+
+Parallel contract
+-----------------
+``DMFieldEvaluate`` internally calls ``DMLocatePoints``, which is
+**COLLECTIVE on the mesh DM's communicator**.  All ranks must call,
+even with zero local points.  Results are rank-local; each rank
+evaluates its own coordinate set and receives its own output arrays.
+
+Unlocated points
+----------------
+Unlocatable points may raise ``PETSc error 62`` or return ``NaN``,
+depending on batching, partition ownership, and the internal
+``DMLocatePoints`` path.  Points that fail the cell-interior test
+(error path: ``dmfieldds.c:261``) propagate as a loud error; points
+absent from the point-location SF graph return ``NaN``.
+
+Outputs are initialised to ``NaN``, so slots that DMFieldEvaluate
+leaves unwritten (the NaN path above) are detectable with ``np.isnan()``.
+"""
+
+import numpy as np
+cimport numpy as np
+
+from petsc4py import PETSc as PyPETSc
+from petsc4py.PETSc cimport DM, PetscDM
+from petsc4py.PETSc cimport Vec, PetscVec
+
+include "../cython/petsc_extras.pxi"
+
+
+cdef class DMFieldEvaluator:
+    """Thin wrapper for a PETSc DMField — create, evaluate, destroy.
+
+    No persistent state beyond a single ``evaluate()`` call.  Create a
+    fresh instance per evaluation; the overhead (~0.5 us) is negligible
+    compared to the evaluate cost (~50-200 us).
+    """
+
+    cdef _p_DMField* _dmf
+
+    def __cinit__(self):
+        self._dmf = NULL
+
+    def create(self, mesh, source_var):
+        """Create DMField from a MeshVariable via ``DMFieldCreateDS``.
+
+        ``mesh.update_lvec()`` must be called (collectively) before this.
+        """
+        cdef _p_DMField* dmf = NULL
+        cdef PetscDM c_dm = (<DM>mesh.dm).dm
+        cdef PetscVec c_lvec = (<Vec>mesh.lvec).vec
+        cdef PetscErrorCode ierr
+
+        ierr = DMFieldCreateDS(c_dm, source_var.field_id, c_lvec, &dmf)
+        CHKERRQ(ierr)
+        self._dmf = dmf
+
+    def evaluate(self, coords_array, mesh, source_var,
+                 gradient=True, hessian=False):
+        """Evaluate at *coords_array*, return ``(B, D, H)``.
+
+        Parameters
+        ----------
+        coords_array : ndarray (n_points, dim)
+            Non-dimensional [0-1] coordinates.
+        mesh : Mesh
+            The mesh (needed for DM/lvec access).
+        source_var : MeshVariable
+            The source variable (metadata only — used for nc/dim).
+        gradient : bool
+            Return first derivatives D (default True).
+        hessian : bool
+            Return second derivatives H (default False).
+
+        Returns
+        -------
+        B : ndarray (n_points, nc) or None
+            Field values.  Unlocated points are ``NaN``.
+        D : ndarray (n_points, nc, dim) or None
+            First derivatives.  ``D[k, j, i]`` = partial(component *j*) /
+            partial x_i.  ``None`` when ``gradient=False``.
+        H : ndarray (n_points, nc, dim, dim) or None
+            Second derivatives.  ``H[k, j, i, l]`` = partial^2(component *j*) /
+            partial x_i partial x_l.  ``None`` when ``hessian=False``.
+        """
+        cdef int n_points = coords_array.shape[0]
+        cdef int dim = mesh.dim
+        cdef int nc = source_var.num_components
+
+        # DMLocatePoints is collective on the mesh DM — all ranks must
+        # participate, even those with zero local points.
+        mesh.dm.localizeCoordinates()
+
+        # Build points Vec on COMM_SELF — each rank owns its own local
+        # coordinates.  DMFieldEvaluate (which calls DMLocatePoints
+        # internally) is collective on the DM; the Vec data is rank-local.
+        pts_np = np.ascontiguousarray(coords_array.ravel(), dtype=np.float64)
+        pts_py = PyPETSc.Vec().createWithArray(pts_np, comm=PyPETSc.COMM_SELF)
+        pts_py.setBlockSize(dim)
+        cdef PetscVec c_pts = (<Vec>pts_py).vec
+
+        # Allocate outputs — initialised to NaN.
+        # If PETSc's evaluator skips unlocated slots (the standard
+        # behaviour of DMFieldEvaluate_DS), the NaN survives and
+        # callers can detect unlocated points with np.isnan().
+        B_arr = (np.full(n_points * nc, np.nan, dtype=np.float64)
+                 if nc > 0 else None)
+        D_arr = (np.full(n_points * nc * dim, np.nan, dtype=np.float64)
+                 if gradient and nc > 0 else None)
+        H_arr = (np.full(n_points * nc * dim * dim, np.nan, dtype=np.float64)
+                 if hessian and nc > 0 else None)
+
+        cdef PetscScalar* B_cptr = \
+            <PetscScalar*>np.PyArray_DATA(B_arr) if B_arr is not None and B_arr.size > 0 else NULL
+        cdef PetscScalar* D_cptr = \
+            <PetscScalar*>np.PyArray_DATA(D_arr) if D_arr is not None and D_arr.size > 0 else NULL
+        cdef PetscScalar* H_cptr = \
+            <PetscScalar*>np.PyArray_DATA(H_arr) if H_arr is not None and H_arr.size > 0 else NULL
+
+        cdef PetscErrorCode ierr
+        ierr = DMFieldEvaluate(
+            self._dmf, c_pts, PETSC_REAL, B_cptr, D_cptr, H_cptr,
+        )
+        pts_py.destroy()
+        CHKERRQ(ierr)
+
+        # ── Reshape to corrected layout ────────────────────────────────
+        # PETSc stores point-major, then component, then direction:
+        #   rD[(i * nc + j) * dimC + l]
+        # So C-order reshape: (n_points, nc, dim) gives D[k, j, i].
+        B_out = B_arr.reshape(n_points, nc) if B_arr is not None else None
+        D_out = (D_arr.reshape(n_points, nc, dim)
+                 if D_arr is not None else None)
+        H_out = (H_arr.reshape(n_points, nc, dim, dim)
+                 if H_arr is not None else None)
+
+        return B_out, D_out, H_out
+
+    def destroy(self):
+        """Release the DMField.  Idempotent — safe to call multiple times."""
+        cdef PetscErrorCode ierr
+        if self._dmf == NULL:
+            return
+        try:
+            ierr = DMFieldDestroy(&self._dmf)
+            CHKERRQ(ierr)
+        except Exception:
+            # Sanctioned failure: during interpreter shutdown PETSc may
+            # have already been finalised; the OS reclaims memory.
+            pass
+        finally:
+            self._dmf = NULL
+
+    def __dealloc__(self):
+        self.destroy()
+
+    def __repr__(self):
+        status = "valid" if self._dmf != NULL else "destroyed"
+        return f"DMFieldEvaluator({status})"
diff --git a/tests/test_0580_dmfield_evaluate.py b/tests/test_0580_dmfield_evaluate.py
new file mode 100644
index 00000000..357ce650
--- /dev/null
+++ b/tests/test_0580_dmfield_evaluate.py
@@ -0,0 +1,384 @@
+"""
+Test DMFieldEvaluate — FE-exact field and gradient evaluation at arbitrary points.
+
+Verifies that ``underworld3.function._dmfield_evaluate.dmfield_evaluate`` returns machine-precision-exact
+values and gradients for FE-representable fields.  Tests the corrected
+``(n, nc, dim)`` layout for vector-gradient tensors, NaN semantics for
+unlocated points, and parallel-safe collective operation.
+
+.. note::
+   DMFieldEvaluate uses PETSc's ``DMFieldEvaluate`` which computes the FE
+   interpolant at each query point.  For fields that are in the FE space
+   (e.g., a quadratic velocity field on P2 elements), the gradient is exact
+   to machine precision at any interior point.
+"""
+import numpy as np
+import pytest
+import underworld3 as uw
+
+from underworld3.function._dmfield_evaluate import dmfield_evaluate
+
+
+# ---------------------------------------------------------------------------
+# Fixtures
+# ---------------------------------------------------------------------------
+
+@pytest.fixture(scope="module")
+def quad_mesh():
+    """Quadratic velocity field (exact on P2) and scalar (exact on P1)."""
+    mesh = uw.meshing.StructuredQuadBox(elementRes=(6, 6))
+    v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2)
+    s = uw.discretisation.MeshVariable("s", mesh, 1, degree=1)
+    # v = (x^2, y^2)  -- gradient: dvx/dx=2x, dvy/dy=2y, cross terms=0
+    v.data[:, 0] = v.coords[:, 0] ** 2
+    v.data[:, 1] = v.coords[:, 1] ** 2
+    # s = x + y     -- gradient: ds/dx=1, ds/dy=1
+    s.data[:, 0] = s.coords[:, 0] + s.coords[:, 1]
+    return mesh, v, s
+
+
+@pytest.fixture(scope="module")
+def tri_mesh():
+    """Unstructured simplex box -- non-affine elements."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.15
+    )
+    v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2)
+    v.data[:, 0] = v.coords[:, 0] ** 2
+    v.data[:, 1] = v.coords[:, 1] ** 2
+    return mesh, v
+
+
+# ---------------------------------------------------------------------------
+# Tests
+# ---------------------------------------------------------------------------
+
+class TestDMFieldEvaluate:
+    """FE-exact field and gradient evaluation."""
+
+    # ---- Basic field values ------------------------------------------------
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_field_values_tri(self):
+        """B (field values) should equal the variable's nodal data.
+
+        Uses an unstructured simplex mesh — DMLocatePoints locates all
+        P2 nodes reliably (0/433 failures on tri vs ~7.6% on quads).
+        """
+        mesh = uw.meshing.UnstructuredSimplexBox(
+            minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.15
+        )
+        v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2)
+        s = uw.discretisation.MeshVariable("s", mesh, 1, degree=1)
+        v.data[:, 0] = v.coords[:, 0] ** 2
+        v.data[:, 1] = v.coords[:, 1] ** 2
+        s.data[:, 0] = s.coords[:, 0] + s.coords[:, 1]
+        # Evaluate at each variable's own nodes
+        B_v, _, _ = dmfield_evaluate(v, v.coords, gradient=False)
+        assert np.allclose(B_v[:, 0], v.coords[:, 0] ** 2, atol=1e-13)
+        assert np.allclose(B_v[:, 1], v.coords[:, 1] ** 2, atol=1e-13)
+        B_s, _, _ = dmfield_evaluate(s, s.coords, gradient=False)
+        assert np.allclose(B_s[:, 0], s.coords[:, 0] + s.coords[:, 1], atol=1e-13)
+
+    # ---- Gradient accuracy -------------------------------------------------
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_gradient_scalar(self, quad_mesh):
+        """Scalar gradient: s = x + y -> ds/dx=1, ds/dy=1."""
+        mesh, _, s = quad_mesh
+        _, D, _ = dmfield_evaluate(s, s.coords, gradient=True)
+        # D[k, j, i]: scalar nc=1, so D[k, 0, 0] = ds/dx, D[k, 0, 1] = ds/dy
+        assert np.allclose(D[:, 0, 0], 1.0, atol=1e-13), "ds/dx"
+        assert np.allclose(D[:, 0, 1], 1.0, atol=1e-13), "ds/dy"
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_gradient_vector_identity(self, quad_mesh):
+        """Vector gradient: v = (x^2, y^2) -> dvx/dx=2x, dvy/dy=2y.
+        
+        Evaluates the P2 velocity gradient at a degree-1 variable's nodes
+        (the natural target for derivative quantities like strain rate).
+        """
+        mesh, v, s = quad_mesh
+        _, D, _ = dmfield_evaluate(v, s.coords, gradient=True)
+        # D[k, j, i] -> D[:, 0, 0] = dvx/dx, D[:, 1, 1] = dvy/dy
+        # Cross terms: D[:, 0, 1] = dvx/dy, D[:, 1, 0] = dvy/dx
+        cx, cy = s.coords[:, 0], s.coords[:, 1]
+        assert np.allclose(D[:, 0, 0], 2 * cx, atol=1e-12), "dvx/dx"
+        assert np.allclose(D[:, 1, 1], 2 * cy, atol=1e-12), "dvy/dy"
+        assert np.max(np.abs(D[:, 0, 1])) < 1e-12, "dvx/dy"
+        assert np.max(np.abs(D[:, 1, 0])) < 1e-12, "dvy/dx"
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_gradient_vector_nontrivial(self):
+        """Non-trivial gradients on unstructured simplex mesh.
+        
+        Evaluates the P2 velocity gradient at a degree-1 variable's nodes.
+        """
+        mesh = uw.meshing.UnstructuredSimplexBox(
+            minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.15
+        )
+        v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2)
+        v.data[:, 0] = v.coords[:, 0] ** 2
+        v.data[:, 1] = v.coords[:, 1] ** 2
+        # Degree-1 variable for evaluating the gradient
+        sr = uw.discretisation.MeshVariable("sr", mesh, 1, degree=1)
+        _, D, _ = dmfield_evaluate(v, sr.coords, gradient=True)
+        cx, cy = sr.coords[:, 0], sr.coords[:, 1]
+        assert np.allclose(D[:, 0, 0], 2 * cx, atol=1e-12), "dvx/dx"
+        assert np.allclose(D[:, 1, 1], 2 * cy, atol=1e-12), "dvy/dy"
+
+    # ---- Asymmetric layout (CRITICAL-2 guard) -----------------------------
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_gradient_asymmetric_layout(self, quad_mesh):
+        """v = (y, 0) — dvx/dy=1, all others 0.  Pins D[k, j, i] layout."""
+        mesh, _, s = quad_mesh
+        v = uw.discretisation.MeshVariable("v_asym", mesh, mesh.dim, degree=1)
+        v.data[:, 0] = v.coords[:, 1]   # vx = y
+        v.data[:, 1] = 0.0              # vy = 0
+        _, D, _ = dmfield_evaluate(v, v.coords, gradient=True)
+        # D[k, j, i]: j=component, i=spatial dim
+        assert np.allclose(D[:, 0, 1], 1.0, atol=1e-13),  "dvx/dy should be 1"
+        assert np.max(np.abs(D[:, 1, 0])) < 1e-13, "dvy/dx should be 0"
+        assert np.max(np.abs(D[:, 0, 0])) < 1e-13, "dvx/dx should be 0"
+        assert np.max(np.abs(D[:, 1, 1])) < 1e-13, "dvy/dy should be 0"
+
+    # ---- Hessian ----------------------------------------------------------
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_hessian_asymmetric_layout(self):
+        """w = (x^2, xy, 0) — in P2, pins H[k, j, i, l] layout.
+        
+        Uses a simplex mesh (DMLocatePoints can locate all P2 nodes) and
+        evaluates the Hessian at a degree-1 variable's nodes.
+        Field is quadratic, so second derivatives are exact for P2.
+        """
+        mesh = uw.meshing.UnstructuredSimplexBox(
+            minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.15
+        )
+        w = uw.discretisation.MeshVariable("w", mesh, 2, degree=2)
+        wc_x, wc_y = w.coords[:, 0], w.coords[:, 1]
+        w.data[:, 0] = wc_x * wc_x               # w0 = x^2
+        w.data[:, 1] = wc_x * wc_y               # w1 = xy
+        # Degree-1 variable for evaluating the Hessian
+        sr = uw.discretisation.MeshVariable("sr", mesh, 1, degree=1)
+        _, _, H = dmfield_evaluate(w, sr.coords, gradient=True, hessian=True)
+        # H[k, j, i, l]: partial^2(component j)/partial x_i partial x_l
+        # w0 = x^2: H[:, 0, 0, 0] = 2, all others 0
+        assert np.allclose(H[:, 0, 0, 0], 2.0, atol=1e-12), "partial^2 w0/partial x^2"
+        assert np.max(np.abs(H[:, 0, 0, 1])) < 1e-12, "partial^2 w0/partial x partial y"
+        assert np.max(np.abs(H[:, 0, 1, 0])) < 1e-12, "partial^2 w0/partial y partial x"
+        assert np.max(np.abs(H[:, 0, 1, 1])) < 1e-12, "partial^2 w0/partial y^2"
+        # w1 = xy: H[:, 1, 0, 1] = H[:, 1, 1, 0] = 1, diagonal = 0
+        assert np.allclose(H[:, 1, 0, 1], 1.0, atol=1e-12), "partial^2 w1/partial x partial y"
+        assert np.allclose(H[:, 1, 1, 0], 1.0, atol=1e-12), "partial^2 w1/partial y partial x"
+        assert np.max(np.abs(H[:, 1, 0, 0])) < 1e-12, "partial^2 w1/partial x^2"
+        assert np.max(np.abs(H[:, 1, 1, 1])) < 1e-12, "partial^2 w1/partial y^2"
+
+    @pytest.mark.level_2
+    @pytest.mark.tier_b
+    def test_hessian_scalar(self, quad_mesh):
+        """Hessian of a linear field is zero."""
+        mesh, _, s = quad_mesh
+        _, _, H = dmfield_evaluate(s, s.coords, gradient=True, hessian=True)
+        # H[k, j, i, l]: scalar nc=1 -> H.shape == (n, 1, dim, dim)
+        assert H.shape == (s.coords.shape[0], 1, mesh.dim, mesh.dim)
+        assert np.max(np.abs(H)) < 1e-10
+
+    # ---- Strain rate from gradient -----------------------------------------
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_strain_rate(self, quad_mesh):
+        """Strain rate computed from D at a degree-1 variable's nodes."""
+        mesh, v, _ = quad_mesh
+        # Create a degree-1 variable to represent where we want strain rate
+        sr = uw.discretisation.MeshVariable("sr", mesh, 1, degree=1)
+        # Evaluate P2 velocity at the degree-1 variable's nodes
+        _, D, _ = dmfield_evaluate(v, sr.coords, gradient=True)
+        # D[k, j, i] -> strain rate symmetric part
+        exx = D[:, 0, 0]
+        eyy = D[:, 1, 1]
+        exy = 0.5 * (D[:, 0, 1] + D[:, 1, 0])
+        # For v = (x^2, y^2): eps_xx = 2x, eps_yy = 2y, eps_xy = 0
+        cx, cy = sr.coords[:, 0], sr.coords[:, 1]
+        assert np.allclose(exx, 2 * cx, atol=1e-12)
+        assert np.allclose(eyy, 2 * cy, atol=1e-12)
+        assert np.max(np.abs(exy)) < 1e-12
+        # Second invariant: sqrt((eps_xx^2 + eps_yy^2) / 2)
+        eII = np.sqrt((exx ** 2 + eyy ** 2) / 2)
+        expected = np.sqrt((4 * cx ** 2 + 4 * cy ** 2) / 2)
+        assert np.allclose(eII, expected, atol=1e-12)
+
+    # ---- Fill mesh variable directly (no Projection) -----------------------
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_fill_mesh_variable(self):
+        """Fill a degree-1 mesh variable from D without uw.systems.Projection.
+        
+        Uses its own mesh to avoid fixture contamination from tests that
+        register extra variables on the shared fixture mesh.
+        """
+        mesh = uw.meshing.StructuredQuadBox(elementRes=(6, 6))
+        v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2)
+        eII_target = uw.discretisation.MeshVariable("eII", mesh, 1, degree=1)
+        v.data[:, 0] = v.coords[:, 0] ** 2
+        v.data[:, 1] = v.coords[:, 1] ** 2
+        # Evaluate v at the target variable's nodes (degree-1 = mesh vertices)
+        _, D, _ = dmfield_evaluate(v, eII_target.coords, gradient=True)
+        eII_target.data[:, 0] = D[:, 0, 0]  # just eps_xx
+        B_check, _, _ = dmfield_evaluate(eII_target, eII_target.coords, gradient=False)
+        assert np.allclose(B_check[:, 0], 2 * eII_target.coords[:, 0], atol=1e-12)
+
+    # ---- Unlocated points (CRITICAL-1 guard) ------------------------------
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_unlocated_points_nan(self, quad_mesh):
+        """Points outside domain return NaN -- no silent garbage."""
+        mesh, v, _ = quad_mesh
+        outside = np.array([[2.0, 2.0], [-1.0, -1.0], [0.5, 2.0]])
+        B, D, H = dmfield_evaluate(v, outside, gradient=True, hessian=True)
+        assert np.all(np.isnan(B)), "B should be NaN for outside points"
+        assert np.all(np.isnan(D)), "D should be NaN for outside points"
+        assert np.all(np.isnan(H)), "H should be NaN for outside points"
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_upper_face_no_crash(self, quad_mesh):
+        """Points near the closed-domain upper face -- no crash, no garbage.
+        
+        Uses y=1.0-1e-10 instead of exact y=1.0 because PETSc 3.25.3
+        ``DMFieldEvaluate_DS`` raises error 62 when ``DMLocatePoints``
+        fails to locate points exactly on the boundary.
+        """
+        mesh, v, _ = quad_mesh
+        eps = 1e-10
+        upper = np.column_stack([np.linspace(0.0, 1.0, 20), np.full(20, 1.0 - eps)])
+        B, _, _ = dmfield_evaluate(v, upper, gradient=False)
+        assert np.all(np.isfinite(B)), "Values near upper face should be finite"
+
+    # ---- No-cache determinism ----------------------------------------------
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_repeated_calls_consistent(self, quad_mesh):
+        """Without cache, repeated calls return identical results.
+        
+        Evaluates at a degree-1 variable's nodes (gradient target).
+        """
+        mesh, v, s = quad_mesh
+        B1, D1, _ = dmfield_evaluate(v, s.coords, gradient=True)
+        B2, D2, _ = dmfield_evaluate(v, s.coords, gradient=True)
+        assert np.allclose(B1, B2)
+        assert np.allclose(D1, D2)
+
+    # ---- Empty coordinates --------------------------------------------------
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_empty_coords(self, quad_mesh):
+        """Zero-point evaluation returns correct-shape empty arrays."""
+        mesh, v, _ = quad_mesh
+        empty = np.empty((0, mesh.dim))
+        B, D, H = dmfield_evaluate(v, empty, gradient=True, hessian=True)
+        assert B.shape == (0, v.num_components)
+        assert D.shape == (0, v.num_components, mesh.dim)
+        assert H.shape == (0, v.num_components, mesh.dim, mesh.dim)
+
+    # ---- Non-square shape (nc != dim) --------------------------------------
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_b
+    def test_non_square_nc_dim(self):
+        """nc=3, dim=2 -> D shape (n, 3, 2) -- catches reshape regression.
+        
+        Uses ``VarType.MATRIX`` with a ``(3, 1)`` shape tuple so that
+        the variable has 3 components in 2D (which would otherwise
+        fail vtype auto-inference).
+        """
+        mesh = uw.meshing.StructuredQuadBox(elementRes=(6, 6))
+        v3 = uw.discretisation.MeshVariable("v3", mesh, (3, 1), degree=1,
+                                             vtype=uw.VarType.MATRIX)
+        cx, cy = v3.coords[:, 0], v3.coords[:, 1]
+        v3.data[:, 0] = cx
+        v3.data[:, 1] = cy
+        v3.data[:, 2] = 0.0
+        _, D, _ = dmfield_evaluate(v3, v3.coords, gradient=True)
+        assert D.shape == (v3.coords.shape[0], 3, 2), \
+            f"Expected (n, 3, 2), got {D.shape}"
+
+    # ---- NULL-pointer paths -------------------------------------------------
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_hessian_only_no_gradient(self, quad_mesh):
+        """gradient=False, hessian=True -- D=NULL path through PETSc."""
+        mesh, v, s = quad_mesh
+        B, D, H = dmfield_evaluate(v, s.coords, gradient=False, hessian=True)
+        assert D is None
+        assert B is not None
+        assert H is not None
+        assert H.shape == (s.coords.shape[0], v.num_components, mesh.dim, mesh.dim)
+
+    @pytest.mark.level_1
+    @pytest.mark.tier_a
+    def test_values_only(self, quad_mesh):
+        """gradient=False, hessian=False -- D=H=NULL path through PETSc."""
+        mesh, v, s = quad_mesh
+        B, D, H = dmfield_evaluate(v, s.coords, gradient=False, hessian=False)
+        assert D is None
+        assert H is None
+        assert B is not None
+
+    # ---- Create/destroy stress test -----------------------------------------
+
+    @pytest.mark.level_2
+    @pytest.mark.tier_b
+    def test_repeated_create_destroy(self):
+        """10 create/destroy cycles -- no PETSc leak or crash."""
+        for i in range(10):
+            mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4))
+            v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2)
+            sr = uw.discretisation.MeshVariable("sr", mesh, 1, degree=1)
+            v.data[:, 0] = v.coords[:, 0]
+            v.data[:, 1] = v.coords[:, 1]
+            # Evaluate gradient at the degree-1 variable's nodes
+            B, D, _ = dmfield_evaluate(v, sr.coords, gradient=True)
+            assert np.allclose(D[:, 0, 0], 1.0, atol=1e-13)
+
+    # ---- MPI parallel test (CRITICAL-3 guard) ------------------------------
+
+    @pytest.mark.level_2
+    @pytest.mark.tier_b
+    @pytest.mark.mpi(min_size=2)
+    def test_parallel_rank_local(self):
+        """Each rank evaluates its own points; zero-point rank doesn't deadlock."""
+        mesh = uw.meshing.StructuredQuadBox(elementRes=(6, 6))
+        v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2)
+        v.data[:, 0] = v.coords[:, 0] ** 2
+        v.data[:, 1] = v.coords[:, 1] ** 2
+
+        # Create degree-1 variable on ALL ranks (DM is shared across ranks)
+        sr = uw.discretisation.MeshVariable("sr", mesh, 1, degree=1)
+
+        rank = uw.mpi.rank
+        if rank == 0:
+            # Evaluate gradient at the degree-1 variable's nodes
+            B, D, _ = dmfield_evaluate(v, sr.coords, gradient=True)
+            assert B is not None
+            assert D is not None
+            cx = sr.coords[:, 0]
+            assert np.allclose(D[:, 0, 0], 2 * cx, atol=1e-12)
+        else:
+            # Rank with zero local points -- must still call (collective)
+            empty = np.empty((0, mesh.dim))
+            B, D, _ = dmfield_evaluate(v, empty, gradient=True)
+            assert B.shape == (0, v.num_components)