element_finder and external points

[gdmcbain]

The behaviour of MeshLine.element_finder is anomalous when some of the points passed to the function returned are outside. For example:

MeshLine().element_finder()(np.array([[-0.3, 0.3, 0.9]]))

returns array([[0, 0]]). The first point is external and silently dropped. The two returned elements are those containing the remaining two points, but how is the caller to know?

MeshLine().element_finder()(np.array([[0.3, 0.9, 1.5]]))

raises IndexError: index 2 is out of bounds for axis 0 with size 2.

More examples:

MeshTri.init_circle().element_finder()(*[np.array([2])]*2)

returns array([36]) while

MeshTri.init_circle().element_finder()(*[np.array([-2])]*2)

returns array([39]) but

 MeshTri().element_finder()(*[np.array([2])]*2)

raises IndexError: index 2 is out of bounds for axis 2 with size 2. Actually MeshTri().element_finder() seems to raise this for any point, even for one that should be internal.

 MeshTri().element_finder()(*[np.array([0.5])]*2)

This arose while looking into the method of characteristics #600, for the feet of characteristics through points near an inlet.

What's a sensible thing to do here? Return a sentinel value like −1 for points for which a containing element can't be found? Raise a specific exception? Return the closest element? (Is that what's happening for MeshTri.init_circle()?)

[gdmcbain]

Return a sentinel value like −1

Not a good sentinel value since it's a valid index into the elements.

[kinnala]

I guess the part with MeshTri().element_finder() is fixed in #667 . There was a bug affecting meshes with less than 5 elements.

[kinnala]

However, there are still no checks if the correct element was found or not.

[kinnala]

Now there is a check in #667 that the correct element was found. It raises ValueError if the point is outside of the mesh. What do you think about that?

[kinnala]

(Note: Work remains for MeshLine)

[kinnala]

Alright, fixed that

[gdmcbain]

In the method of characteristics #600 , it's not an error for an evaluation point to fall outside the domain. One follows a material particle back from a quadrature point at the end of a time step to find where it was at the start and evaluates some field (say temperature) there that the particle will carry with it. A very normal occurrence is that a fluid particle will flow into the domain through an inlet in which case the 'foot' of the characteristic is contained in no element and the evaluation by interpolation is replaced typically by an assignment of the property as in a Dirichlet inlet boundary condition. If each point were handled separately, raising an exception for one upstream of the inlet wiuld be reasonable ; it could be caught and specially handled but if the 'feet' of all characteristics ending at quadrature points were sought at once with a single call to the element finder as usual in scikit-fem, wouldn't a ValueError thwart the general functioning for the interior points?

If instead a sentinel value for the relevant entries in the array were returned that could easily be detected and treated specially and appropriately.

However, i'm getting the feeling that the method of characteristics might not be a good choice for strong advection in scikit-fem. I hope to investigate other techniques (dG ?), in which case this objection would be less relevant.

[gdmcbain]

Besides the method of characteristics though, there might be other applications where a subset of points falling outside the domain are special but not erroneous; e.g. tracking particles, some of which eventually hit a wall or leave via an outlet. Also moving or deforming meshes?

[kinnala]

I'd simply for-loop over the points with a try-catch. You need to provide some example snippet so I know what is expected. Sentinel integer is not so simple because it is always a valid index.

[gdmcbain]

Sentinel integer is not so simple because it is always a valid index.

Doesn't an int greater than or equal to the length of the dimension of the array raise an IndexError?

You need to provide some example snippet so I know what is expected.

O. K. Basically I think there's one point for each quadrature point in the domain, being the estimate of where the fluid that is at the quadrature point at the end of the time-step was at the beginning of the time-step, given the advecting velocity field.

Since this morning, I've thought of an example problem for the method of characteristics #600 that might not be too arduous: a closed cavity (maybe like ex18) (so that all the feet should be within the domain and the present issue won't arise) with steady flow (so that the feet are constant from time-step to time-step and the probes matrix only needs to be assembled once) and unsteady advection or advection–diffusion with a given nonuniform initial scalar field.

I'd simply for-loop over the points with a try-catch

Yeah, O. K., I could do that, loop over the quadrature points; there are a few of them and I usually avoid big for-loops in Python, but it should work. I had been thinking of a single call to probes (or a call per time-step in case of an unsteady advecting velocity field).

[kinnala]

O. K. Basically I think there's one point for each quadrature point in the domain, being the estimate of where the fluid that is at the quadrature point at the end of the time-step was at the beginning of the time-step, given the advecting velocity field.

Where does element_finder come into the picture? Is the time step so large that the point can move over several elements?

[kinnala]

What is being done about a point passing the boundary of the domain?

[gdmcbain]

Where does element_finder come into the picture? Is the time step so large that the point can move over several elements?

I see, right. No, i wouldn't think one would choose a time-step as large as that. If i recall correctly the method is unconditionally stable, but still for practical accuracy that's not going to be much good. So yes, i had been thinking of using Basis.probes (because it's to hand) but now that you put it like that, perhaps something more specialized might be better : first try the same element, then the one adjacent, and failing that quite possibly abandon, raising an exception, or adaptively reducing the time-step.

Thanks.

This discussion should be in #600 rather than here.

[gdmcbain]

What is being done about a point passing the boundary of the domain?

I imagine the usual high Peclet number approximation is that an incoming point gets the value of the advected quantity assigned.