InteriorBasis.point_source

[gdmcbain]

Fix #629 and exemplify.

This won't work for MeshLine until #628 is resolved, but is fine for MeshTri, &c.

Should probably compare the finite element solution with Green's function, as discussed in the docstring. For a test, I'd suggest one dimension, pending #628.

[kinnala]

Is this docstring true for other than Lagrange elements?

[gdmcbain]

Hmm, I think so, in the Galerkin sense. Initially I reread Jan's point-source function in #561 and was reminded of what we'd done for the probe for Nico. Then i thought that at the continous variational level, interpolation is the inner product of a finite element function with a Dirac delta and a point source is the inner product of a Dirac delta with a test function. Those are adjoints and on discretization they're transposes. That line of thoight doesn't make any reference to the txpe of basis.

Similarly for skfem.models.poisson.unit_load and the functional that integrates, replacing the Dirac delta with the constant 1.

[gdmcbain]

I don't have a reference for this. I couldn't find one in a quick look. I remember that for years users were asking about point sources on other finite element msiling lists. Somtimes the discussion suggested mollification: convolving the point source with a Gaussian kernel, but I didn't like that so much as it still didn't seem well suited to the quadrature. Then FreeFEM added one recently; I refrained from looking at their implementation as it's GPL.

[gdmcbain]

Possible hole in the argument: the argument shows that the transpose property holds but it only follows from that that the function returns the Galerkin point source vector if the probes function works correctly. I didn't look into whether that works for all kinds of bases. Does it? I assumed that it did.

[gdmcbain]

A couple of brief references to point sources in some of the older FEM books on my self:

• Segerlind, L. J. (1976). Applied Finite Element Analysis. New York: Wiley, §8.6 ‘Point sources’
• Hughes, T. J. R. (2000). The Finite Element Method. Mineola, New York: Dover, exercise 3.11.2(b)

Unfortunately Hughes doesn't give the solution to the exercise (delta function loading on a one-dimensional, quadratic, three-node element), but I think it's looking the similar to my pair of linear two-node elements.

Segerlind notes:

The integral of a function multiplied by an impulse function, however, is simply the function evaluated at X₀ and Y₀.

which is elementary and the subsequent example is only for ElementTriP1.

I'm yet to turn up anything on nonlagrangian elements.

[gdmcbain]

from pathlib import Path

import numpy as np

import skfem
from skfem.models import laplace
from skfem.visuals.matplotlib import plot, savefig

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

lap = skfem.asm(laplace, basis)
print(lap.toarray())

a = 3 / 4
rhs = basis.point_source(np.array([a]))
print("RHS:", rhs)
condensed = skfem.condense(lap, rhs, D=basis.find_dofs())
print(condensed[0].toarray())
print(condensed[1])

solution = skfem.solve(*condensed)
print(solution)

ax = plot(basis, solution, Nrefs=3)
ax.set_title(f"{type(basis.elem)}")
x = np.linspace(mesh.p[0].min(), mesh.p[0].max())
ax.plot(x, np.stack([(1 - a) * x, (1 - x) * a]).min(0), linestyle="--", label="exact")
ax.plot(
    mesh.p[0],
    solution[basis.nodal_dofs[0]],
    marker="o",
    linestyle="None",
    label="nodal DoFs",
)
ax.legend()
savefig(Path(__file__).with_suffix(".png"))

I think that that's doing the right thing. The plot looks a bit funny though. I copied from

scikit-fem/docs/examples/ex34.py

Line 47 in 13b5d69

plot(basis, x, Nrefs=3)

ElementLineMini looks O. K. too.

[gdmcbain]

I mean Hermite elements are not a good choice for this problem since they're C¹ and the solution isn't, but they're trying:

Actually

condensed = skfem.condense(lap, rhs, D=basis.find_dofs())

isn't the right boundary condition for Hermite either, that's for a built-in Euler–Bernoulli beam, isn't it.

[gdmcbain]

Perhaps I'll look at the elastic plate with ElementTriMorley too, adapting ex18.

from skfem import *
from skfem.helpers import ddot, dd

import numpy as np


mesh = MeshTri.init_circle(4)
element = ElementTriMorley()
mapping = MappingAffine(mesh)
ib = InteriorBasis(mesh, element, mapping, 2)


@BilinearForm
def biharmonic(u, v, w):
    return ddot(dd(u), dd(v))


stokes = asm(biharmonic, ib)
psi = solve(*condense(stokes, ib.point_source(np.array([0.2, 0.3])), D=ib.find_dofs()))


if __name__ == "__main__":
    from os.path import splitext
    from sys import argv
    from skfem.visuals.matplotlib import draw, show
    from matplotlib.tri import Triangulation

    M, Psi = ib.refinterp(psi, 3)

    ax = draw(mesh)
    ax.tricontour(Triangulation(*M.p, M.t.T), Psi)
    name = splitext(argv[0])[0]
    show()

[kinnala]

Looks good. Ready for merge or still a draft?

[kinnala]

I've actually used this kind of load in some of the articles on Kirchhoff plates and always hacked it together in nonreusable way. This would've really simplified that. Thanks!

[gdmcbain]

Looks good. Ready for merge or still a draft?

I should add a test.

I'll add a comparison with the Green's function to the example later.

I assume the Green's function for the Kirchhoff plate is known too. Have you got a reference to hand? Actually that might make a nicer example than the membrane, showing right away that it works for other than Lagrange elements.

[kinnala]

I assume the Green's function for the Kirchhoff plate is known too. Have you got a reference to hand?

I think there is but I have no reference at hand.

On the other hand, there is an infinite series solution available for a square plate and a point load. The series is not converging very quickly but it's something if you truncate it to a sum. See Section 5.1 (eq. 5.13) in this preprint: https://arxiv.org/pdf/1707.08396.pdf

[gdmcbain]

Hang on, something goes wrong if I change from MeshLine(np.linspace(0, 1, ...)) to MeshLine().refined(...).

[gdmcbain]

>>> m = fe.MeshLine().refined()
>>> fe.InteriorBasis(m, fe.ElementLineP1()).probes(np.array([[3/4]])).toarray()
array([[-0.5,  0. ,  1.5]])
>>> fe.InteriorBasis(fe.MeshLine(m.p[0]), fe.ElementLineP1()).probes(np.array([[3/4]])).toarray()
array([[0. , 0.5, 0.5]])

[gdmcbain]

It looks like MeshLine.element_finder #417 #628 again.

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

[kinnala]

It looks like MeshLine.element_finder #417 #628 again.

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

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

[gdmcbain]

Let's take that up in #633.

[kinnala]

Let me add, sort_t=False is something I'd prefer to have everywhere but it wasn't easy to drop it from MeshTri1 because other code assumes it. So maybe sorting in element_finder is sufficient?

[gdmcbain]

It remains to include the example in tests/test_examples and before that to put something testable in the example, e.g. comparison of the solution with the analytic Green's function (once I track down a reference for that—basically it's the usual logarithmic source minus an image sink on the other side of the circumference).

[gdmcbain]

I still haven't been able to find a better reference for the Green's function (Sadybekov et al just quote the result without reference or derivation), but this is it anyway.

[gdmcbain]

The PNG above will do for the Gallery.

[kinnala]

The PNG above will do for the Gallery.

Sure! Can you add it to docs/listofexamples.rst?