fix(interpolation): bypass DMLocatePoints when caller supplies a verified hint

src/underworld3/function/_dminterp_wrapper.pyx

@@ -135,12 +135,23 @@ cdef class CachedDMInterpolationInfo:
                 DMInterpolationDestroy(&self._ipInfo)
                 raise RuntimeError(f"DMInterpolationAddPoints failed with error {ierr}")
 
-        # Set up — calls DMLocatePoints which is COLLECTIVE on the mesh DM.
+        # Set up — calls DMLocatePoints (or, on simplex / manifold meshes,
+        # the DMLocatePoints-bypass path) which is COLLECTIVE on the mesh DM.
         # All ranks must call this, even with zero local points.
         # ignoreOutsideDomain=1: PETSc silently skips points it cannot locate
         # rather than crashing. This is essential — points_in_domain() uses a
         # kd-tree heuristic that can misclassify distant points as interior.
-        if n_points > 0:
+        #
+        # Hint policy: pass the UW3 kdtree-nearest-cell hint to PETSc only on
+        # meshes where it can be trusted as authoritative — i.e. simplex
+        # cells (planar faces, affine map, exact face-containment test) and
+        # manifold meshes (dim != cdim, where PETSc's own in-cell test is
+        # unreliable near 2-manifold simplex edges). On non-simplex meshes
+        # (quads, hexes), deformed cells can have non-planar faces and the
+        # kdtree-nearest cell can be wrong; pass NULL so PETSc's
+        # DMLocatePoints runs its own (more rigorous) cell-location search.
+        cdef bint use_hint = bool(mesh.dm.isSimplex()) or (mesh.dim != mesh.cdim)
+        if n_points > 0 and use_hint:
             cells_view = np.ascontiguousarray(self.cells)
             ierr = DMInterpolationSetUp_UW(self._ipInfo, dm, 0, 1, <size_t*> &cells_view[0])
         else:

src/underworld3/function/_function.pyx

@@ -1151,8 +1151,13 @@ def petsc_interpolate(   expr,
             # CACHE MISS - Create structure and cache it
             cached_info = CachedDMInterpolationInfo()
 
-            # Get cell hints
-            # coords is already np.ndarray type (function signature ensures this)
+            # Cell hints from a kdtree-nearest-cell lookup. The downstream
+            # `create_structure` decides — based on the mesh's cell shape —
+            # whether to pass these to PETSc as authoritative (bypass path,
+            # simplex / manifold meshes whose face containment is exact for
+            # affine transforms) or as a guess for `DMLocatePoints` to
+            # verify (non-simplex meshes whose deformed faces can be
+            # non-planar and require PETSc's pseudo-inverse projection).
             cells = mesh.get_closest_cells(coords)
 
             # Create and set up DMInterpolation structure (EXPENSIVE)

src/underworld3/function/petsc_tools.c

@@ -22,18 +22,17 @@ PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool
   PetscScalar       *a;
   PetscInt           p, q, i;
   PetscMPIInt        rank, size;
-  Vec                pointVec;
+  Vec                pointVec = NULL;
   PetscSF            cellSF;
   PetscLayout        layout;
   PetscReal         *globalPoints;
-  PetscScalar       *globalPointsScalar;
+  PetscScalar       *globalPointsScalar = NULL;
   const PetscInt    *ranges;
   PetscMPIInt       *counts, *displs;
   const PetscSFNode *foundCells;
   const PetscInt    *foundPoints;
   PetscMPIInt       *foundProcs, *globalProcs;
   PetscInt           n, N, numFound;
-  PetscErrorCode     ierr;
 
   PetscFunctionBegin;
   PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
@@ -66,41 +65,45 @@ PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool
   PetscCall(PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs));
   /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
 #else
-  #if defined(PETSC_USE_COMPLEX)
-  PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar));
-  for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
-  #else
-  globalPointsScalar = globalPoints;
-  #endif
-  PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec));
   PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs));
   for (p = 0; p < N; ++p) foundProcs[p] = size;
   cellSF = NULL;
-  /* the Underworld code is used to find good guesses for the owning cells */
-  if (owning_cell)
-  {
-    PetscSFNode *sf_cells;
-    ierr = PetscMalloc1(N, &sf_cells);
-    CHKERRQ(ierr);
-    size_t range = 0;
-    for (size_t p = 0; p < (size_t)N; p++)
-    {
-      sf_cells[p].rank = 0;
-      sf_cells[p].index = owning_cell[p];
-      if (owning_cell[p] > range)
-      {
-        range = owning_cell[p];
-      }
+  if (owning_cell) {
+    /*
+      Bypass DMLocatePoints when the caller supplies a hint.
+
+      The caller is expected to call this path only on meshes where the
+      hint is authoritative for cell containment — simplex cells (planar
+      faces, affine reference map) and 2-manifold meshes (where PETSc's
+      DMPlexLocatePoint_Simplex_2D_Internal pseudo-inverse projection is
+      itself unreliable near triangle edges in 3-D). On those meshes the
+      kdtree-nearest cell from `Mesh.get_closest_cells` correctly
+      contains every query coord, so PETSc's re-verification adds no
+      information and only introduces wrong-cell fallback opportunities.
+
+      Each rank claims every point whose hint is a valid cell id (we
+      treat the size_t value (size_t)-1 as a non-claim sentinel). The
+      MPI_Allreduce(MIN, foundProcs) below then assigns each global
+      point to the lowest-rank claimant.
+    */
+    numFound    = 0;
+    foundPoints = NULL;
+    foundCells  = NULL;
+    for (p = 0; p < N; ++p) {
+      if (owning_cell[p] != (size_t)-1) foundProcs[p] = rank;
     }
-    ierr = PetscSFCreate(PETSC_COMM_SELF, &cellSF);
-    CHKERRQ(ierr);
-    // PETSC_OWN_POINTER => sf_cells memory control goes to cellSF
-    // nroots must be > max(iremote.index), so use range + 1
-    ierr = PetscSFSetGraph(cellSF, range + 1, N, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER);
-    CHKERRQ(ierr);
+  } else {
+    /* Build pointVec lazily — only the DMLocatePoints path needs it. */
+    #if defined(PETSC_USE_COMPLEX)
+    PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar));
+    for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
+    #else
+    globalPointsScalar = globalPoints;
+    #endif
+    PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec));
+    PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF));
+    PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells));
   }
-  PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF));
-  PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells));
 #endif
 
   /*
@@ -176,9 +179,12 @@ PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool
 #else
   PetscCall(PetscFree2(foundProcs, globalProcs));
   PetscCall(PetscSFDestroy(&cellSF));
-  PetscCall(VecDestroy(&pointVec));
+  /* pointVec / globalPointsScalar are only constructed in the !owning_cell branch above. */
+  if (!owning_cell) {
+    PetscCall(VecDestroy(&pointVec));
+    if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar));
+  }
 #endif
-  if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar));
   if (!redundantPoints) PetscCall(PetscFree3(globalPoints, counts, displs));
   PetscCall(PetscLayoutDestroy(&layout));
   PetscFunctionReturn(PETSC_SUCCESS);