PointLocator parallel logic issues

[permcody]

The framework contains several bad parallel constructs of this form:

          // First find the element the hit lands in
          const Elem * elem = (*pl)(some_point);

          if (elem && elem->processor_id() == some_proc.processor_id())

This turns out to be unreliable when points happen to fall on edges or nodes between processor boundaries. Reference this ticket when fixing these issues:

Related: #4407 #3431 #2899 #4439

[permcody]

Possible candidates:

permcj@thearchitect { ~/projects/moose(devel) }$ ack -Q "elem && elem->processor_id()"
framework/src/postprocessors/PointValue.C
66:  if (elem && elem->processor_id() == processor_id())

framework/src/transfers/MultiAppVariableValueSamplePostprocessorTransfer.C
70:          if (elem && elem->processor_id() == from_mesh.processor_id())

framework/src/transfers/MultiAppVariableValueSampleTransfer.C
76:          if (elem && elem->processor_id() == from_mesh.processor_id())

framework/src/vectorpostprocessors/PointSamplerBase.C
100:  if (elem && elem->processor_id() == processor_id())

[YaqiWang]

Seems to me locator should return a vector of elements. The size could be zero (outside the domain), one (inside an element), two (inside an element side), more than two (on lines or vertices). How this vector will be used is up to the users.

[jwpeterson]

On Fri, Dec 19, 2014 at 2:10 PM, Yaqi notifications@github.com wrote:

Seems to me locator should return a vector of elements. The size could be zero (outside the domain), one (inside an element), two (inside an element side), more than two (on lines or vertices). How this vector will be used is up to the users.

The PointLocators are based on a binary/quad/oct tree search which returns as soon as it finds an element or fails to find any element...

So, what you suggest could be implemented, but would require rewriting the internals of the PointLocator class significantly.

John

[YaqiWang]

The users can check if the neighboring elements own the point or not. Here the neighboring elements are all elements connected to the elements through not only sides but lines, vertices. I bet it will be really painful to do so...

[permcody]

We just have to be careful not to assume that you always receive the same element on all processors. That may not always be necessary depending on what you want or need to do but when your parallel communication relies on a specific id for broadcasting - you can get in trouble. It's up to the user

[YaqiWang]

If we can check if the point is a vertex of an element or in a side or line, then we can go through the connectivity to find all the other elements owning the point. This seems quite doable.

[friedmud]

Like Cody mentions... This is implantation specific. Some places in the code might not care if we get different elements back on different processors. However, in the places where that is important we need to be more careful.

Let me do some thinking on this... it's nobtrivial to get this right and make it fast...

[friedmud]

Lol - should be implementation... damn iPhone autocorrect!

[permcody]

"nobtrivial" wow - a 2-for-1 in that comment!

[permcody]

We'll close this when we are sure that the bug has been removed.

[friedmud]

Well - this ticket is about more than what my PR was for. My PR only fixed up this type of stuff in one specific place. There are more in the framework...

[roystgnr]

"This seems quite doable." Doable and mostly done, @YaqiWang - see https://github.com/libMesh/libmesh/blob/master/include/geom/elem.h#L323 for a method you can use given a point and its PointLocator result to see what other elements contain it as well.

What I mean by "mostly": By default we use a relative tolerance of TOLERANCE in the contains_point() tests there; let me know if you want to make that adjustable.

[YaqiWang]

Thanks @roystgnr, glad to know libmesh has such a function already.

[roystgnr]

Actually, come to think of it, the libMesh function isn't always what you need here. If your mesh is a simple manifold containing the point you're looking for, then find_point_neighbors will be able to find all other elements containing that point by searching along that manifold.

If however you're searching for a point along a "slit" where two manifold boundaries overlap, or a point that's overlapped by elements of more than one dimension, etc., then you have to use the allowed_subdomains argument to PointLocator to select which set of elements you might hit, and then find_point_neighbors will only give you other elements from the same local manifold as the first one you found.