Interpolation nonmatching mesh in parallel

Issue #162 new
Mikael Mortensen created an issue

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)

  1. Garth N. Wells

    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?

  2. Mikael Mortensen reporter

    New version using bounding boxes. Algorithm now like this for interpolation from Function u0 to Function u1

    1. 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.
    2. Create a map from dof to component number in Mixed Space.
    3. Create bounding boxes for the partitioned mesh of u0 and distribute to all processors.
    4. Using bounding boxes compute the processes that may own the interpolation points of u1.
    5. Distribute to potential owners and subsequently try to evaluate u0. If successful, return value of u0 to owner.
  3. Mikael Mortensen reporter

    Yes, resolved for Lagrange spaces at least. Should we close it or reformulate that not all elements are handled yet?

  4. Log in to comment