How to match function values on boundary

Question

Given two functions f and f_b, defined on a mesh m and a mesh mb = BoundaryMesh(m, 'exterior'), is there a way to "copy" the values on the boundary from f_b to f?

Using f.interpolate(f_b) fails with "Unable to evaluate function at point" (presumably because it attempts to evaluate f_b in the interior of mb).

When using a CG_1 function space for both f and f_b, I was able to do manually copy the affected dofs on the boundary vertices using the mb.entity_map, but I don't see a way to generalize this to different function spaces.

Does someone know a way to do this?

mikael-mortensen · Answer 1 · Mar 15, 2014

Hi,

Like oyvinev kindly suggests, fenicstools is definitely an option for Lagrange spaces. I thought I should fill in with some more details though. Consider the function spaces you mention

from dolfin import *

mesh  = UnitCubeMesh(10, 10, 10)
bmesh = BoundaryMesh(mesh, 'exterior')
V = FunctionSpace(mesh, 'CG', 1)
V_b = FunctionSpace(bmesh, 'CG', 1)
f_b = interpolate(Expression('x[0]'), V_b)

You can now use fenicstools like this

from fenicstools import interpolate_nonmatching_mesh
f = interpolate_nonmatching_mesh(f_b, V)

There is another very nice option, though, using overloading of an Expression eval

class fb(Expression):
    def eval(self, values, x):
        try:
            values[0] = f_b(x)
        except:
            values[0] = 0

The very nice thing about this approach is that you can use the small boundary Function f_b in a Dirichlet boundary condition for the large function space V even though f_b is defined on a different mesh:

bc = DirichletBC(V, fb(), ...)

This may be possible for other function spaces as well, but I haven't really looked into it.

Finally I include a complete running example to show how I usually set the velocity profile of my Navier Stokes solver on an inlet of arbitrary shape:

from dolfin import *

# Create mesh and boundary mesh
mesh = UnitCubeMesh(10, 10, 10)
bmesh = BoundaryMesh(mesh, 'exterior')

# Create SubDomain and mark left hand side
Left = AutoSubDomain(lambda x, on_bnd: near(x[0], 0))
cf = CellFunction('size_t', bmesh, 0)
Left.mark(cf, 1)
smesh = SubMesh(bmesh, cf, 1)

# Solve a Poisson equation that descibes fully developed duct flow
Vs = VectorFunctionSpace(smesh, 'CG', 1)
us = TrialFunction(Vs)
vs = TestFunction(Vs)

dpdx = Constant((10.0, 0, 0))
a = inner(grad(us), grad(vs))*dx
L = dot(dpdx, vs)*dx
bc = DirichletBC(Vs, Constant((0, 0, 0)), "on_boundary")
A = assemble(a)
b = assemble(L)
bc.apply(A, b)
us = Function(Vs)
solve(A, us.vector(), b)

plot(us, title='Inlet profile computed in 2D')

# Now use this inlet profile as DirichletBC for the 3D problem
V = VectorFunctionSpace(mesh, 'CG', 1)
u = Function(V)

# Two options - first is a memory friendly overloaded Expression
class InletVelocity(Expression):
    def eval(self, values, x):
        try:
            values[:] = u(x)
        except:
            values[:] = 0
    def value_shape(self):
        return (3,)

inletvelocity = InletVelocity()
bcb = DirichletBC(V, inletvelocity, Left) 
bcb.apply(u.vector()) # Works

# Second option - interpolate 2D Function on 3D Function. Memory
# demanding and not working with regular fenics
#us = interpolate(us, V)

# fenicstools working though
from fenicstools import interpolate_nonmatching_mesh
us = interpolate_nonmatching_mesh(us, V)
bc2 = DirichletBC(V, us, Left)
bc2.apply(u.vector())
plot(u, title="Inlet profile in 3D")

interactive()

How to match function values on boundary

2 Answers