invalid interpolator of MeshTri returned by refinterp

[gdmcbain]

While developing an example of two-dimensional creeping flow #58 , I ran into trouble trying to evaluate the stream-function at a point. (The problem there has a known closed-form solution and should equal to 1/64 at the origin, which would seem to be a convenient thing to test.)

In ex02, a MeshTri M is generated from the Morley InteriorBasis ib with refinterp, but the M.interpolator doesn't seem to work.

To reproduce:

import ex02
M, X = ex02.ib.refinterp(ex02.x, 3)
M.interpolator(X)

raises

RuntimeError: In initialize: Triangulation is invalid

I didn't find any examples using MeshTri.interpolator but

ex12.m.interpolator(ex12.x)(0, 0)

returns array(0.24952565) which is close to the exact solution, ¼.

[gdmcbain]

Ah: I see that the MeshTri has coincident points

M.p[:, (0, 135)]

array([[0., 0.],
          [0., 0.]])

[kinnala]

refinterp cuts the connectivity between the neighboring elements to enable the visualisation of discontinuous bases, etc.

So in practice, the resulting mesh has duplicate nodes. This might be one reason for the error, as I believe the resulting mesh is passed directly to matplotlib.

Edit: you already noticed ;-)

[gdmcbain]

O. K. then. So that seems reasonable. Maybe there should just be something in the docstring of either InteriorBasis.refinterp or MeshTri.interpolator or both warning against doing this.

For refinterp, I guess that that might be that the returned Mesh isn't necessarily valid for all purposes, maybe isn't 'conforming'. Something like from Gmsh's §1.2 Mesh: finite element mesh generation

A finite element mesh is a tessellation of a given subset of the three-dimensional space by elementary geometrical elements of various shapes (in Gmsh’s case: lines, triangles, quadrangles, tetrahedra, prisms, hexahedra and pyramids), arranged in such a way that if two of them intersect, they do so along a face, an edge or a node, and never otherwise.

or Braess's (2001, Finite Elements, 2nd edn, §5.1) ‘If $T_i\cap T_j$ consists of exactly one point, then it is a common vertex of $T_i$ and $T_j$, …’.

I note that ex02.ib.refinterp(ex02.x, 3)[0]._validate() raises

 UserWarning: Mesh._validate(): Mesh contains duplicate vertices.

For MeshTri.interpolator, that the MeshTri must be valid in some sense. I'm not sure exactly what the conditions of validity are here, maybe need to see what happens inside matplotlib.tri.LinearTriInterpolator.

[gdmcbain]

For the underlying problem though—how to interpolate the deflection defined on Morley elements—the answer is that this isn't the way to do it but then how to. I'll raise another issue for that.