Interpolation onto nonmatching mesh should set function values to zero at nodes outside the original mesh

Fenicstools has a function "interpolate_nonmatching_mesh" (see here) which can interpolate a function defined on one FunctionSpace onto another (unrelated) FunctionSpace, even if they are defined on different meshes. From what I gather, some (most?) of this functionality has been incorporated into dolfin. The dolfin demo demo_nonmatching-interpolation.py illustrates this. However, the dolfin interpolation cannot deal with the case that the second mesh (onto which the function is interpolated) has vertices outside the domain of the first mesh. On the other hand, Fenicstools simply sets the values to zero at these vertices, which is very handy in certain situations. Could this be adopted into dolfin?

Below is a modified version of the demo illustrating this issue (works with dolfin 1.5.0 and Fenicstools 1.5.0). Change the value of "allow_extrapolation" and switch between the two different interpolation methods to see the difference between the interpolation methods in dolfin and Fenicstools.

from __future__ import print_function
from dolfin import *
from fenicstools.Interpolation import interpolate_nonmatching_mesh

# Change 'allow_extrapolation' to 'True' to avoid crash below
# (but it then gives wrong results).
parameters.allow_extrapolation = False

# Create mesh and define function spaces
mesh0 = UnitSquareMesh(16, 16)
mesh1 = RectangleMesh(0, 0, 2, 1, 64, 64)

V0 = FunctionSpace(mesh0, "CG", 3)
V1 = FunctionSpace(mesh1, "CG", 1)

# Define function
v0 = interpolate(Expression("sin(10.0*x[0])*sin(10.0*x[1])",
                            element=FiniteElement('CG', triangle, 3)),
                 V0)
v1 = Function(V1)

# (i) Interpolate using dolfin. This crashes if 'allow_extrapolation'
# is False and gives wrong results otherwise.
v1.interpolate(v0)

# (ii) Interpolate using Fenicstools. This gives sensible results, since
# it sets the function values at vertices in mesh1 which are outside of
# mesh0 to zero.
#
#v1 = interpolate_nonmatching_mesh(v0, V1)

# Plot functions
plot(v0, title="v0")
plot(v1, title="v1")
interactive()

Comments (1)