Impose full gradient BC via Lagrange multipliers

Hi everyone,

I have a two dimensional physical problem governed by a poisson equation

- nabla^2 u = f(x,y) in Gamma

I have measurements of the two gradient components du/dx and du/dy along some line, gamma_m.
This field I am trying to propagate into the far field.

The tangential gradient component I can just integrate and then enforce through a Dirichlet condition. In order to impose the normal component I gather I will have to go via Lagrange multipliers.

So I think the variational form is then

m = du/dn_measured

\int u v dx - \lambda \int ( du/dn - m ) v d(gamma_m) = \int f v dx

I spent a long time trying to piece it together from the information on related problems such as
https://answers.launchpad.net/fenics/+question/212434
https://answers.launchpad.net/fenics/+question/216323

but still cannot see where I am going wrong. My solution attempt so far is:

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(64, 64)
V = FunctionSpace(mesh, "Lagrange", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R

# Create mesh function over cell facets
boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)

# Mark left boundary facets as subdomain 0
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0] < DOLFIN_EPS

Gamma_Left = LeftBoundary()
Gamma_Left.mark(boundary_parts, 0)

class FarField(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ( (x[0] > 1.0-DOLFIN_EPS) \
               or (x[1]<DOLFIN_EPS) or (x[1]> 1.0-DOLFIN_EPS) )

Gamma_FF = FarField()
Gamma_FF.mark(boundary_parts, 1)

# Define boundary condition
u0 = Expression("sin(x[1]*pi)")
bcs = [DirichletBC(V, u0, Gamma_Left)]

# Define variational problem
(u, lmbd) = TrialFunctions(W)
(v, d) = TestFunctions(W)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Constant(0.0)
h = Constant(-4.0)
n = FacetNormal(mesh)

F = inner(grad(u), grad(v))*dx + d*dot(grad(u),n)*ds(0) + lmbd*dot(grad(v),n)*ds(0)-\
   (f*v*dx + g*v*ds(1) + h*d*ds(0) + lmbd*h*ds(0))

a = lhs(F)
L = rhs(F)

# Compute solution
A = assemble(a, exterior_facet_domains=boundary_parts)
b = assemble(L, exterior_facet_domains=boundary_parts)
for bc in bcs: bc.apply(A, b)

w = Function(W)
solve(A, w.vector(), b, 'lu')
(u,lmbd) = w.split()

# Plot solution
plot(u, interactive=True)

which gives a noisy result not at all resembling any solution to a poisson equation.
I would appreciate any help or pointers in the right direction - many thanks already!
Cheers
Markus