-
Notifications
You must be signed in to change notification settings - Fork 36
Interpolation nonmatching mesh in parallel
With the function interpolate_nonmatching_mesh_any
it is possible to interpolate from a Function of any space on one mesh to a Function of any space on a different mesh. An example for a Nedelec element is shown below.
from dolfin import *
from fenicstools import interpolate_nonmatching_mesh_any
import numpy
# Create two different meshes on the same domain
mesh1 = UnitSquareMesh(2, 2)
mesh2 = UnitSquareMesh(4, 4)
V1 = FunctionSpace(mesh1, "Nedelec 1st kind H(curl)", 2)
V2 = FunctionSpace(mesh2, "Nedelec 1st kind H(curl)", 2)
f = Expression(("x[0]", "x[1]"), degree=1)
u1 = interpolate_nonmatching_mesh_any(f, V1)
# Interpolate to nonmatching mesh
u2 = interpolate_nonmatching_mesh_any(u1, V2)
# Validation
u3 = interpolate_nonmatching_mesh_any(f, V2)
assert numpy.allclose(u2.vector().get_local(), u3.vector().get_local())
For Lagrange spaces this type of interpolation is possible to perform in parallel with dolfin's LagrangeInterpolator
class. For special elements, like Nedelec, there is currently no other possible option than the one provided here. Note that for interpolate_nonmatching_mesh_any
to work, the evaluate_dofs
function must be implemented for the space interpolating from. This means that for now the function will not work for enriched elements like the bubble.
The interpolation function may also be used to interpolate between meshes of different geometrical dimension. To illustrate, this code interpolates from a function in 3D to a 2D slice through the center (at z=0.5) of the 3D domain. Note that this second example does not work in parallel because the dolfin Function SubMesh does not work in parallel. However, this is only a limitation of the SubMesh and not the interpolation routine.
mesh1 = UnitCubeMesh(12, 12, 12)
V1 = FunctionSpace(mesh1, "Nedelec 1st kind H(curl)", 2)
f = Expression(("sin(x[0]*pi)", "cos(x[1]*pi)", "x[2]"))
u1 = interpolate_nonmatching_mesh_any(f, V1)
# Create a 2D unit square mesh from the boundary of mesh1
mesh2 = BoundaryMesh(mesh1, "exterior")
side = CellFunction("size_t", mesh2, 0)
AutoSubDomain(lambda x: x[2] < 1e-8).mark(side, 1)
mesh3 = SubMesh(mesh2, side, 1)
x = mesh3.coordinates()
x[:, 2] += 0.5
# Interpolate from the 3D Function to the 2D Function
V2 = FunctionSpace(mesh3, "Nedelec 1st kind H(curl)", 2)
u2 = interpolate_nonmatching_mesh_any(u1, V2)
plot(u1)
plot(u2)