Interpolation nonmatching mesh in parallel
Hi,
I have some code I have been using for a while that is part of fenicstools and that has been requested by several. The code is basically a stand-alone function that interpolates a dolfin Function from one mesh to a Function with the same FunctionSpace defined on another mesh.
mesh0 = UnitSquareMesh(10, 10)
mesh1 = UnitSquareMesh(20, 20)
V0 = FunctionSpace(mesh0, "CG", 1)
V1 = FunctionSpace(mesh1, "CG", 1)
u0 = interpolate(Expression("x[0]*x[1]"), V0)
u0.update()
u1 = interpolate_nonmatching_mesh(u0, V1)
The new feature is that it works in parallel. I think this functionality should be part of dolfin. It could be added to the Function class.
I attach a cpp file that contains the implementation. It can be compiled using instant, or compile_extension_module. A complete test is
from dolfin import *
interpolation_code = open('interpolation_nonmatching.cpp', 'r').read()
compiled_fem_module = compile_extension_module(code=interpolation_code)
def interpolate_nonmatching_mesh(u0, V):
"""Interpolate from GenericFunction u0 to FunctionSpace V.
The FunctionSpace V can have a different mesh than that of u0, if u0
has a mesh.
"""
u = Function(V)
compiled_fem_module.interpolate_nonmatching_mesh(u0, u)
return u
mesh0 = UnitSquareMesh(10, 10)
mesh1 = UnitSquareMesh(20, 20)
V0 = FunctionSpace(mesh0, "CG", 1)
V1 = FunctionSpace(mesh1, "CG", 1)
u0 = interpolate(Expression("x[0]*x[1]"), V0)
u0.update()
u1 = interpolate_nonmatching_mesh(u0, V1)
plot(u0, title="mesh0")
plot(u1, title="mesh1")
Q0 = FunctionSpace(mesh0, "CR", 1)
Q1 = FunctionSpace(mesh1, "CR", 1)
Mixed0 = V0*Q0
Mixed1 = V1*Q1
m0 = interpolate(Expression(("x[0]", "x[1]")), Mixed0)
m0.update()
m1 = interpolate_nonmatching_mesh(m0, Mixed1)
plot(m0[1], title="mesh0")
plot(m1[1], title="mesh1")
interactive()
Comments (10)
-
-
reporter Ok, will do tomorrow.
-
reporter - attached interpolation_nonmatching2.cpp
New version using bounding boxes. Algorithm now like this for interpolation from Function u0 to Function u1
- Create a map from all u1's dofs' coordinates to dofs living on that coordinate. This is done such that one only needs to visit (and distribute) each interpolation point once.
- Create a map from dof to component number in Mixed Space.
- Create bounding boxes for the partitioned mesh of u0 and distribute to all processors.
- Using bounding boxes compute the processes that may own the interpolation points of u1.
- Distribute to potential owners and subsequently try to evaluate u0. If successful, return value of u0 to owner.
-
- changed milestone to 1.5
-
Isn't this resolved by pull request #112?
-
reporter Yes, resolved for Lagrange spaces at least. Should we close it or reformulate that not all elements are handled yet?
-
- changed milestone to 1.6
-
Issue
#411was marked as a duplicate of this issue. -
- changed milestone to 1.7
-
- removed milestone
Removing milestone: 1.7 (automated comment)
- Log in to comment
This should be proposed and discussed on the mailing list, and include a description of the implementation/methodology. Can you also provide some performance measures/scaling?