diff --git a/src/underworld3/function/_dminterp_wrapper.pyx b/src/underworld3/function/_dminterp_wrapper.pyx index 0423c751..cc6c2d90 100644 --- a/src/underworld3/function/_dminterp_wrapper.pyx +++ b/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, &cells_view[0]) else: diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index c17126cd..ec4892b5 100644 --- a/src/underworld3/function/_function.pyx +++ b/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) diff --git a/src/underworld3/function/petsc_tools.c b/src/underworld3/function/petsc_tools.c index e8f054f4..d6b56771 100644 --- a/src/underworld3/function/petsc_tools.c +++ b/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);