# 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

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 < 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 > 1.0-DOLFIN_EPS) \
or (x<DOLFIN_EPS) or (x> 1.0-DOLFIN_EPS) )

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

# Define boundary condition
u0 = Expression("sin(x*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.5, 2) + pow(x - 0.5, 2)) / 0.02)")
g = Constant(0.0)
h = Constant(-4.0)
n = FacetNormal(mesh)

(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

## Question information

Language:
English Edit question
Status:
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Johannes Ring (johannr) said on 2013-05-14: #1