Bug in LagrangeInterpolator

[matteo croci]

For high degree Lagrange elements (and under some other unclear conditions) LagrangeInterpolator introduces some error in the output function. I am attaching a script that illustrates this behaviour.

Oddly enough, if you fiddle a bit with the input parameters (i.e. lower the element degree, change the domain sizes, change the number of elements in the mesh), this behaviour might disappear (this is why the script attached is not really "minimal").

from dolfin import *
import mshr

# Create two nested meshes
sp = mshr.Rectangle(Point(-2,-2), Point(2,2))
rect = mshr.Rectangle(Point(-0.5,-0.5), Point(0.5,0.5))

sp.set_subdomain(1,rect)
mesh = mshr.generate_mesh(sp, 50)
cf = MeshFunction("size_t", mesh, 2, mesh.domains())
inner_mesh = SubMesh(mesh, cf, 1)

OUT = FunctionSpace(mesh, 'CG', 5)
IN  = FunctionSpace(inner_mesh, 'CG', 5)

# u_in1 will be interpolated using LagrangeInterpolator
# u_in2 will be interpolated using interpolate
u_in1 = Function(IN)
u_in2 = Function(IN)

expr = Expression('cos(x[0]*x[1])', degree = 6)

u_out = interpolate(expr, OUT)

interp = LagrangeInterpolator()

interp.interpolate(u_in1, u_out)

# We need to allow extrapolation here as there are two points
# which are not recognised to be inside the domain (although they
# are and the meshes are nested)
u_out.set_allow_extrapolation(True)
u_in2 = interpolate(u_out, IN)

# compare with exact value (computed numerically)
ex_sol = 9.965342812790287e-01
values = (assemble(u_in1*dx), assemble(u_in2*dx))
print "exact solution: %e, values (LagrangeInterpolator, interpolate): %s" % (ex_sol, values)
print "errors [LagrangeInterpolator, interpolate]: %s" % [abs(item - ex_sol) for item in values]

[Patrick Farrell]

Matteo's problem works with the interpolation matrix code we merged earlier this week. I suggest we just replace the LagrangeInterpolator code with code based on create_transfer_matrix (wherever it ends up). @garth-wells , what do you think?

[Prof Garth Wells]

Yes, I think LagrangeInterpolator should be replaced by code that underpins the construction of the interpolation matrix.

[Chris Richardson]

Agree, but new code needs a clean up.

[Mikael Mortensen]

This does not really look like a LagrangeInterpolator issue. Looks like the dirty hack allow_extrapolation fixes the problem for regular interpolate, and this dirty hack should also fix the problem for LagrangeInterpolator. We’re probably seeing the old issue with roundoff, as discussed here: https://bitbucket.org/fenics-project/dolfin/issues/168/evaluation-of-function-at-specific-point.

What’s the create_transfer_matrix? Does it do the same as LagrangeInterpolator, only by explicitly creating an interpolation matrix? Is it as fast, works in parallel and does it also work for non-Lagrange spaces?

[Patrick Farrell]

@chris_richardson What are the main points left to clean up? To generalise it to arbitrary backends?

@mikael_mortensen Matteo wrote some code (PETScDMCollection.cpp:131) to explicitly construct an interpolation matrix between two function spaces. It works in parallel, for non-nested meshes and for non-nested domains; we haven't tested its speed. It relies on pointwise evaluation, so it works for the same function spaces that LagrangeInterpolator does.

We wrote this for geometric multigrid (see demo/undocumented/gmg-poisson) but wonder if it could be used for interpolation more generally.

[Jan Blechta]

Remove passed milestone.

[Jan Blechta]

should be replaced by code that underpins the construction of the interpolation matrix

I don't think that it is a good reason to remove the code - on the contrary it's an advantage for specific applications. I would expect some further discussion here.

[Prof Garth Wells]

We should avoid the unnecessary creation of an assembled sparse matrix.

There should certainly be some scope to share code for the different applications.

[Jan Blechta]

@garth-wells you speak in puzzles. Do you claim the LagrangeInterpolater is too much domain specific to be included in DOLFIN?

[Prof Garth Wells]

@blechta: @pefarrell is advocating that we use an assembled matrix to compute the interpolation (the same matrix that is used for multigrid). My point is that using an assembled matrix for plain interpolation is overkill, but that the two operations have commonalities.