[petsc-users] Finding which cell an arbitrary point belongs to in DMPlex

[Matthew Knepley]

On Fri, Sep 25, 2020 at 7:44 AM Jeremy Theler <jeremy at seamplex.com> wrote:

> On Thu, 2020-09-24 at 18:30 -0400, Matthew Knepley wrote:
> > > There is also DMPlexFindVertices() which finds the nearest vertex
> > > to
> > > > the given coords in the given radius.
> > >
> > > At first I had understood that this function performed a nearest-
> > > neighbor search but after a closer look it sweeps the local DM and
> > > marks whether the sought points coincide with a mesh node within
> > > eps or
> > > not. Neat.
>
> This DMPlexFindVertices() sweeps over DMGetCoordinatesLocal() which
> returns both the local and ghost coordinates, so at the end of the day
> I might get more than one process claiming to have found the same node.
>
> How can I ignore ghost points so each vertex actually belongs to the
> process that found it?
>

This is just a debugging function from Vaclav. It is not intended to be used
for production code. If you wanted to change it, you can just have it search
DMGetCoordinates(). It is consistent with Plex since mesh points exist in
the
local space, and you only have a notion of global with the SF. I think the
right
modification would be to optionally pass the pointSF and not return points
contained in it.


> > > > I wrote it some time ago mainly for debug purposes. It uses just
> > > > brute force. I'm not sure it deserves to exist :-) Maybe we
> > > should
> > > > somehow merge these functionalities.
> > >
> > > It works, although a kd-tree-based search would be far more
> > > efficient
> > > than a full sweep over the DM.
> >
> > We should not need to do that. LocatePoints() does not sweep the
> > mesh.
> > It just does grid hashing. kd is a little better with really
> > irregular distributions,
> > but hashing should be fine.
>
> Yes, it seems to be pretty efficent (although there is no support for
> 3D so far).
>
> Two more things about plexgeometry.c:
>
>  1. shouldn't line 224 in DMPlexLocatePoint_Simplex_3D_Internal()
> compare against -eps instead of against zero as donde in line 145 in
> DMPlexLocatePoint_Simplex_2D_Internal()?
>

I hate to do it, but I made the change. It is practical. The Right Thing to
do is use the robust predicates
from Shewchuk and Bailey. That is on my list for someday.


>  2. wouldn't it be better to replace DMGetDimension() by
> DMGetCoordinateDim() in line 45 inside DMPlexFindVertices? I have a 2D
> mesh with 3D coordinates and PetscUnlikely() is triggered.
>

You are right. I fixed it. I will put this branch in today hopefully.

  Thanks,

     Matt


> Thanks
> --
> jeremy
>
-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200925/569c0a9f/attachment.html>