Hello!
My question is close related with the topic
http://scicomp.stackexchange.com/questions/7204/poisson-equation-impose-full-gradient-as-boundary-condition-via-lagrange-multip
which was discussed in another forum. Shortly, I need to impose a general non-linear boundary condition (given as a restriction R(u) = 0) on certain boundary and the best way is to define a Lagrange multiplier field (lm) on the boundary. As a test, I supose my restriction is a simple Dirichlet bc "u = u_bot" => "R(u) = u - u_bot" on a specific booundary. My implementation of the problem is the following:
from dolfin import *
import numpy
set_log_level(ERROR)
#-------------- Preprocessing step -----------------
# Create mesh and define function space
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'Lagrange', 1)
#print V.print_dofmap()
# Define boundary segments for the boundary condition via Lagrange multiplier
# Create mesh function over cell facets
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
# Mark upper boundary facets as subdomain 0
class UpperBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[1] - 1) < tol
Gamma_Up = UpperBoundary()
Gamma_Up.mark(boundary_parts, 0)
# Mark lower boundary facets as subdomain 1
class LowerBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[1]) < tol
def map(self, x, y):
y = x
Gamma_Low = LowerBoundary()
Gamma_Low.mark(boundary_parts, 1)
# Create function space over Lower boundary edges???
V_lm = FunctionSpace(mesh, 'Lagrange', 1, constrained_domain = Gamma_Low)
#print V_lm.print_dofmap()
#Mixed finite element space
W = V * V_lm
print W
#-------------- Solution and problem definition step -----------------
# given mesh and boundary_parts
u_up = Constant(373.15)
bcs = [DirichletBC(W.sub(0), u_up, Gamma_Up)]
# Define variational problem
f_e = Expression("1000")
u_bot = Expression("273.15")
#expression = -inner(grad(u), grad(v))*dx - v*f_e*dx - v*f_e*ds
(u , lm) = TrialFunctions(W)
(v, v_lm) = TestFunctions(W)
expression = -inner(grad(u), grad(v))*dx + (lm*v + v_lm*(u-u_bot))*ds(1) - v*f_e*dx
a = lhs(expression)
L = rhs(expression)
# 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)
U = Function(W)
solve(A, U.vector(), b, 'lu')
u, lm = U.split()
#Plotting
plot(u,interactive=True)
plot(lm,interactive=True)
file = File("Temperature.pvd")
file << u
For the formulation works, "lm" and its test functions "v_lm" must be defined only on the boundary ds(1). That its, "lm" must have a DofMap restricted only to the mesh nodes lying on the boundary "1". This problem run cleanly for coarse meshes (7x7) but for finer the solution is "nan" (see the file after run). I suspect that the system is singular, because the DofMap for "lm" is actually not restricted to the specified boundary; as result, the formulation is giving zero rows in the system matrix because there is nothing to assembly in the interior nodes for the variable "lm".
So, my questions are:
- Am I doing a right definition of the "constrained_domain" space for the Lagrange Multipliers field? I use a identity mapping "x = y" in the class LowerBoundary, because I do not wish to do any special mapping with the boundary.
- There is another way to define a Lagrange multiplier field restricted only to a specific boundary of the domain?
Thanks!
Diego
