interpolate non-Lagrange finite element function from coarse mesh to refined mesh in parallel

Question

If we have a finite element function in some FunctionSpace and then refine the mesh and create a new FunctionSpace with the refined mesh, we can use the interpolate function to interpolate the function from the first space into the second or vice versa. Unfortunately, this works in serial, but not in parallel. (In parallel, if the "allow_extrapolation" parameter is not set, the code crashes, while if it is set, the code completes but gives the wrong answer.) This is discussed, for example, here.

If the FunctionSpaces use Lagrange elements, then instead of the usual interpolate function one can use LagrangeInterpolator().interpolate, and this gives correct results, also in parallel.

My question is: for non-Lagrange spaces, e.g., Nedelec spaces, is it possible to interpolate from a coarse mesh space to a refined mesh space in parallel or is this just not available in FEniCS?

Answer 1

I think you are right (that there is no way of doing it at present).

Since there is a defined parent <-> child relation between refined meshes, there should be a function like "interpolate_parent()" which works cell-wise and would work in parallel.

Answer 2

Normally, yes, this is not possible for non-Lagrange space. However, there is a (perhaps ugly) hack that I can think of. Use LagrangeInterpolator to interpolate first from Nedelec on original mesh onto a Lagrange space on the new mesh. Then interpolate with regular interpolate function from the intermediate Lagrange space onto Nedelec on the new mesh. A script may explain better:

from dolfin import *

mesh1 = UnitSquareMesh(40, 40)   # Fine mesh
mesh2 = UnitSquareMesh(30, 30)   # Some other mesh
N1 = FunctionSpace(mesh1, "N1curl", 1)
N2 = FunctionSpace(mesh2, "N1curl", 1)
CG2 = VectorFunctionSpace(mesh2, "CG", 2)  # Intermediate space on mesh2

# Just create some simple solution on N1
e = Expression(("x[0]", "x[1]"))
u_n1 = interpolate(e, N1)

# Now we want u_n1 interpolated on N2. 
u_cg = Function(CG2)
u_n2 = Function(N2)

# Make first intermediate interpolation
lp = LagrangeInterpolator()
lp.interpolate(u_cg, u_n1)

# Move from u_cg to final N2
u_n2 = interpolate(u_cg, N2)
plot(u_n2)
interactive()

Comments

Actually this does not work. The problem is that the N1curl space is not contained in the Lagrange CG2 space, because N1curl vector fields are, in general, not continuous. Therefore if we interpolate an N1curl function to a CG2 space, and then interpolate that to an N1curl space on a different mesh, we will usually not get the same result as if we interpolated the first N1curl function directly onto the second N1curl space.

Here is a simple code that demonstrates that the problem. We compute the interpolant from coarse mesh to fine mesh N1curl space first in the way you suggest (passing through CG2), and then directly. The output shows that the two results are nowhere near equal.

from dolfin import *
import numpy as np
mesh1 = UnitSquareMesh(1, 1)
mesh2 = UnitSquareMesh(2, 2)
N1 = FunctionSpace(mesh1, "N1curl", 1)
N2 = FunctionSpace(mesh2, "N1curl", 2)
CG2 = VectorFunctionSpace(mesh2, "CG", 2)
u_n1 = Function(N1)
u_n1.vector()[:] = np.array([0., 1., 2., 3., 4.])
u_cg = Function(CG2)
lp = LagrangeInterpolator()
lp.interpolate(u_cg, u_n1)
u_n2 = interpolate(u_cg, N2)
u_n2a = interpolate(u_n1, N2)
print np.linalg.norm(u_n2.vector().array() - u_n2a.vector().array())

Output is 0.768..., not near zero.

---

Ah, ok, an ugly hack that didn't work then:-( It seemed to work for my linear Expression and your testcase above at least gives me output 1e-14 if using N1curl of order one on mesh2 as well. Does it generally work for that particular case?

---

Yes, I overstated the case. I believe that your approach does work for the case of interpolating
N1curl on a mesh to N1curl on a refined mesh. Thanks.