MeshLine.refine

[gdmcbain]

While working on the point source #632, trouble was encountered with MeshLine().refined() and this was traced to element_finder:

>>> m = MeshLine().refined()
>>> m.element_finder()(np.array([.7]))
array([0])
>>> MeshLine(m.p[0]).element_finder()(np.array([.7]))
array([1])

I'm not sure which of these is right. I mean the first is wrong since 0.7 is definitely in the second element, but I'm not sure about the second. Are the points passed to the constructed assumed increasing?

[gdmcbain]

There is Mesh.sort_t flag which defaults to False. MeshTri1 has sort_t=True which is explained in the docstring of Mesh.

But given m = MeshLine().refined(), should m be valid? It has

 >>> m.p
array([[0. , 1. , 0.5]])
>>> m.t
array([[0, 2],
       [2, 1]])

which looks O. K., not sorted, but O. K., but then m.element_finder() is wrong.

>>> m.element_finder()(np.array([.2]))
array([1])
>>> m.element_finder()(np.array([.7]))
array([0])

[kinnala]

So element_finder seems to do argsort so that's not the issue.

[kinnala]

I think there is some sort of assumption inside finder regarding the order of self.t.

[kinnala]

It does assume that the start point of interval is self.t[0].

[kinnala]

There is Mesh.sort_t flag which defaults to False. MeshTri1 has sort_t=True which is explained in the docstring of Mesh.

But given m = MeshLine().refined(), should m be valid? It has

 >>> m.p
array([[0. , 1. , 0.5]])
>>> m.t
array([[0, 2],
       [2, 1]])

which looks O. K., not sorted, but O. K., but then m.element_finder() is wrong.

>>> m.element_finder()(np.array([.2]))
array([1])
>>> m.element_finder()(np.array([.7]))
array([0])

The original version before #630 gave the correct result here. I'll have to check #630 carefully.

[kinnala]

Alright, so it was a change from self.t[0] to self.t[1]. Maybe we shouldn't assume anything.

[kinnala]

Does #634 fix the issue in the point source?

[gdmcbain]

No, it gives different answers for

mesh = skfem.MeshLine(np.linspace(0, 1, 2 ** 1 + 1))

and

mesh = skfem.MeshLine().refined(1)

and indeed

>>> MeshLine(np.array([0, .5, 1])).element_finder()(np.array([.7]))
array([0])
>>> MeshLine().refined(1).element_finder()(np.array([.7]))
array([1])

The first of these is wrong, since

>>> m = MeshLine(np.array([0, .5, 1]))
>>> m.p[0, m.t]
array([[0. , 0.5],
       [0.5, 1. ]])

(Both meshes look the same in terms of m.p[0, m.t].)