I have a physical problem governed by the Poisson equation in two dimensions $$ -\nabla^2 u = f(x,y), \; in \; \Omega $$ I have measurements of the two gradient components $\partial{u}/\partial{x}$ and $\partial{u}/\partial{y}$ along some part of the boundary, $\Gamma_m$, so would like to impose $$ \frac{\partial{u}}{\partial{x_i}}_0 = g_m, \; on \; \Gamma_m $$ and propagate into the far field.
The tangential gradient component, $\frac{\partial{u}}{\partial{x}}_{(t,0)}$, I can just integrate and then enforce through a Dirichlet condition, such that $$ \int_{\Gamma_m}\frac{\partial{u}}{\partial{x}}_{(t,0)} \, ds = u_0 $$ In order to simultaneously impose the normal component, $\frac{\partial{u}}{\partial{x}}_{(n,0)}$, I gathered I would have to go via Lagrange multipliers.
So I think the variational form is then $$ \int_\Omega \nabla{u} \cdot\nabla{v}\, dx - \lambda \int_{\Gamma_m} ( \frac{\partial{u}}{\partial{x}}_{(n,0)}-g_m ) v \,ds = \int_\Omega 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 runs but gives a noisy result not at all resembling a solution to a Poisson equation.
It seems to have something to do with the combined function spaces, but I cannot find the mistake.
I would appreciate any help or pointers in the right direction - many thanks already!
Cheers
Markus