Add dmfield_evaluate from PETSc

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,
+    ),
 ]
 
 

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

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

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})"