Working with BoundaryMesh

Question

I would like to have a form like

bf * v * ds(1)

where v is the TestFunction over a 2D domain, ds(1) refers to a particular boundary of that domain, and bf is a 1D function defined only on that boundary. Here is an example code

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(20, 20)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
leftboundary = LeftBoundary()

boundaries = FacetFunction('size_t', mesh)
boundaries.set_all(0)
leftboundary.mark(boundaries, 1)
ds = Measure('ds')[boundaries]

boundarymesh = BoundaryMesh(mesh, 'exterior')
left_mesh = SubMesh(boundarymesh, leftboundary)
L = FunctionSpace(left_mesh, 'CG', 1)
bf = Function(L)
# here goes some more code to find the correct bf
# but here a dummy function will do
bf.vector()[:] = np.ones(21, dtype=float)

V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
L = bf * v * ds(1)
assemble(L)

The error message I get for the last line (version 1.4.0) is:

AssertionError: Expecting a completed form with domains at this point.

I feel like I should somehow map the function defined on the boundary bf to the 2D mesh, but I can't find an example of this somewhere.

Comments

I noticed a similar (unanswered) question, so for reference:
http://fenicsproject.org/qa/6810/vertex-mapping-from-submesh-boundarymesh-back-actual-mesh

Answer 1

The following solution was pointed out to me:

from dolfin import *
parameters["allow_extrapolation"] = True
...
L = interpolate(bf, V) * v * ds(1)
assemble(L)