interpolate Morley

[gdmcbain]

The problem is to interpolate a function defined on Morley elements at a specified (x, y) point which isn't necessarily a node; e.g. the deflection at the centre of the Kirchhoff plate in ex02, or the value of the stream-function at the centre of the circular-cylindrical cavity in ex18.

Attempting MeshTri.interpolator with the MeshTri returned by refinterp fails #59; the mesh is invalid, having duplicate points.

I think interpolation involves two steps:

1. Find the element (if any) that contains the point.
2. Evaluate the function based on the local basis functions of that element and the elemental degrees of freedom.

and I suspect that in general the first is the harder.

Thus far, this operation is only implemented in scikit-fem as MeshTri.interpolator (P1) and MeshTri.const_interpolator (P0) and makes use of :

• matplotlib.tri.Triangulation,
• matplotlib.tri.LinearTriInterpolator
• matplotlib.tri.Triangulation.get_trifinder

Perhaps the TriFinder might be factored out of those so that it could also be used for other elements defined on linear simplicial meshes; e.g., besides Morley, also affine P2. Looking ahead, with MeshQuad._splitquads, perhaps quadrilateral meshes too.

[gdmcbain]

Finding the triangle:

from matplotlib.tri import Triangulation
from ex02 import m

triang = Triangulation(m.p[0, :], m.p[1, :], m.t.T)
trifinder = triang.get_trifinder()

m.p[:, m.t[:, trifinder(.5, .5)]]
array([[0.5   , 0.4375, 0.5625],
           [0.5   , 0.4375, 0.4375]])

which looks like the vertices of a triangle containing the specified point. Here trifinder(.5, .5) is array(8), so it's m.t[:, 8].

Next: how to evaluate the shape function of element 8 at the point.

[kinnala]

Next, one should use Mapping.invF and Element.gbasis, quite similarly as here:

scikit-fem/skfem/assembly/global_basis/facet_basis.py

Line 99 in 167ba29

Y = self.mapping.invF(x, tind=self.tind)

Instead of leaving the different evaluated basis functions separate (as is done in the previous url), we should sum their values similarly as in GlobalBasis.interpolate.

scikit-fem/skfem/assembly/global_basis/global_basis.py

Line 182 in 167ba29

for j in range(self.Nbfun):

Yesterday I finished a big project proposal so now I hope to have more time again.

[gdmcbain]

Thanks for the tips; I haven't fully understood the machinery here but continuing to scratch around:

ex02.x[ex02.ib.element_dofs[:, 8]]

is

array([ 2.33941283e-04, 2.30040321e-04, 2.55625370e-04, 3.45343932e-04, -2.10424724e-04, -9.12993610e-05])

and

 ex02.x[ex02.ib.element_dofs[:, 8]] @ [ex02.ib.basis[0][i][8] for i in range(ex02.ib.Nbfun)]

is

array([0.00023885, 0.00023825, 0.00025078])

which might be the ordinates at the three quadrature points of the element containing (0.5, 0.5).

I think the first three Morley degrees of freedom are also nodal ordinates which looks about right too.

[gdmcbain]

ex02.mapp.invF(np.array([0.5, 0.5]).reshape(-1, 1, 1), np.array([8]))

is

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

which looks like the local coordinates of the point in the reference element.

[kinnala]

Yes. I apologize, these parts of the library are not yet fully documented as I wanted to improve the documentation of user-facing parts.

[gdmcbain]

It might be that in part but also I think that interpolation isn't intrinsic to the finite element method which I associate with linear and bilinear forms, approximated by quadrature. The modules in place so far seem to reflect that, as they should. It's not clear to me at this stage whether further extrinsic functionality needs to be added, or whether everything that's needed for interpolation is already there.

[gdmcbain]

which looks like the local coordinates of the point in the reference element.

Confirming that with a further experiment, computing the local coordinates of a global point in all of the elements and seeing for which it lies in the reference triangle ξ ≥ 0, η ≥ 0, ξ + η ≤ 1.

from functools import reduce
import numpy as np
import ex02

xieta = ex02.mapp.invF(np.array((.5, .5)).reshape(-1, 1, 1))[:, :, 0]
elements = np.nonzero(reduce(np.logical_and, [xieta[0, :] >= 0, xieta[1, :] >= 0, xieta[0, :] + xieta[1, :] <= 1]))
reduce(np.intersect1d, ex02.m.t[:, elements].T)

returns array([4]) and indeed

ex02.m.p[:, reduce(np.intersect1d, ex02.m.t[:, elements].T)]

is array([[0.5], [0.5]]).

[gdmcbain]

I got this going in 672a50c in #66 for the stream-function at the centre of the cavity in ex18. It did use Mapping.invF and ElementH2.gbasis, as suggested.

It does involve recourse to lower level methods than usual and a fair bit of shape-manipulation so might benefit from encapsulation later.

I don't think it's too specific to ElementTriMorley so might well work for ElementTriP2 too.

[kinnala]

Yes, it looks like it could work for any triangular element and also could be vectorized to work with multiple points and multiple triangles at once.

I'll consider whether to replace MeshTri.interpolator or do something else.

[gdmcbain]

Yes, on that vectorization to multiple elements, I had wanted to pass the elements containing the specified points as the tind keyword instead of omitting that, letting it default to all the elements, and then extracting the ones (one) of interest by indexing:

672a50c#diff-a4c8fc413d1defaa1a6464dbfbce0f52R64

    psi0 = (np.hstack([element.gbasis(mapping, local_point, k)[0][triangle]
                                for k in range(ib.Nbfun)])

but I couldn't work out yet the shape wanted for tind.

[kinnala]

This is work in progress:

points = np.array([[0.45, 0.78], [0.23, 0.55]]).T
triangulation = Triangulation(mesh.p[0, :], mesh.p[1, :], mesh.t.T)
triangles = triangulation.get_trifinder()(points[:, 0], points[:, 1])
local_points = mapping.invF(points[:, :, np.newaxis], tind=triangles)

Next step is to evaluate the basis at local_points.

[kinnala]

Maybe this would work:

def interpolator(mesh, basis, dofs):
    from matplotlib.tri import Triangulation
    
    triang = Triangulation(mesh.p[0, :], mesh.p[1, :], mesh.t.T)
    finder = triang.get_trifinder()
    
    def interpfun(x):
        tris = finder(*x)
        pts = basis.mapping.invF(x[:, :, np.newaxis], tind=tris)
        w = np.zeros(x.shape[1])
        for k in range(basis.Nbfun):
            phi = basis.elem.gbasis(basis.mapping, pts, k, tind=tris)
            w += dofs[basis.element_dofs[k, tris]]*phi[0].flatten()
        return w
    
    return interpfun

# For example,
# f = interpolator(mesh, ib, psi)
# f(points)

[gdmcbain]

Thanks for looking into this.

I wonder if the use of scipy.spatial.Delaunay instead of matplotlib.tri.Triangulation would generalize this from MeshTri to any simplicial mesh (i.e. really only MeshTet but maybe also MeshLine). Would this function then be defined as a method of a mixin class MeshSimplicial?

[gdmcbain]

points = np.array([[0.45, 0.78], [0.23, 0.55]]).T
triangulation = Triangulation(mesh.p[0, :], mesh.p[1, :], mesh.t.T)
triangles = triangulation.get_trifinder()(points[:, 0], points[:, 1])
local_points = mapping.invF(points[:, :, np.newaxis], tind=triangles)

isn't quite right since if given just a single global point, it returns two local points

points = np.array([[0.45], [0.23]]).T
triangles = triangulation.get_trifinder()(points[:, 0], points[:, 1])
mapping.invF(points[:, :, np.newaxis], tind=triangles)

returns [[[-1.59661657], [ 0.94083269]], [[ 2.9832518 ], [ 1.92616944]]] with shape (2, 2, 1).

[gdmcbain]

function interpolator doesn't work either. I'm not sure what shape to pass points to pass to it, but it has to be (2, number_of_points), no?

But then (in ex18)

interpolator(mesh, ib, psi)(np.array([[0.], [0.]]))

or

interpolator(mesh, ib, psi)(mesh.p[:, :3])

raises

ValueError: operand has more dimensions than subscripts given in einstein sum, but no '...' ellipsis provided to broadcast the extra dimensions.

This is from the line

phi = basis.elem.gbasis(basis.mapping, pts, k, tind=tris)

[kinnala]

Oops. I tested this only with subclasses of ElementH1. The implementation of gbasis is different in ElementH2 so it should be taken into account.

[gdmcbain]

Ah, and ElementH1.gbasis doesn't memoize the inverse Vandermonde. #68

By the way: how do you mark up those tidy quotes from sources as in #62 (comment) ?

[kinnala]

When you read the source code in Github and mouse over the line number, three dots apper next to it. Pressing that gives you a URL.

[kinnala]

Thanks for looking into this.

I wonder if the use of scipy.spatial.Delaunay instead of matplotlib.tri.Triangulation would generalize this from MeshTri to any simplicial mesh (i.e. really only MeshTet but maybe also MeshLine). Would this function then be defined as a method of a mixin class MeshSimplicial?

In particular, which function do you mean? Delaunay.find_simplex?

[gdmcbain]

Sorry, scipy.spatial was a false lead. Yes, I was thinking of that method but it only seems to work with the triangulation generated by its own Delaunay class, I couldn't see how to pass in the one in MeshTri.

[kinnala]

For triangular meshes, you can now try InteriorBasis.interpolator. I'm gonna write some tests for it now.

[kinnala]

Ok, there are still some issues. I'm trying to fix them asap.

[kinnala]

cc99ad9 should fix interpolating ElementH2 bases. Support for multiple elements was missing in MappingAffine.F.

[gdmcbain]

Excellent! This works now, as demonstrated in #66 for the stream-function in the centre of the cavity in ex18. Thank you!

[gdmcbain]

I note that the new InteriorBasis.interpolator works for ElementTriP2 too.