Solving laplace equation using lagrange multipliers for dirichlet bc

Question

Solve

Laplace(u) = 0 in [0,1]x[0,1]

with

u = x^2 - y^2 on the boundary

We use Lagrange multipliers to implement dirichlet bc. The stress on the boundary is computed by a post-processing step in the end. The scheme is explained on page 400-401 of this paper

Gunzburger and Hou,
TREATING INHOMOGENEOUS ESSENTIAL BOUNDARY CONDITIONS IN FINITE ELEMENT METHODS AND THE CALCULATION OF BOUNDARY STRESSES
SIAM J. Num. Anal., vol. 20, no. 2, pp. 390-424, April 1992

The code is given below.

import scipy.sparse as sps
from dolfin import *

parameters.linear_algebra_backend = "uBLAS"

class boundary(SubDomain):
   def inside(self,x,on_boundary):
      return on_boundary

uexact = Expression("x[0]*x[0] - x[1]*x[1]")

n = 50
mesh = UnitSquareMesh(n,n)
V = FunctionSpace(mesh, 'CG', 1)

u = TestFunction(V)
v = TrialFunction(V)

ubc = Function(V)

# Project boundary condition
M = assemble(u*v*ds)
b = assemble(uexact*v*ds)
bd = boundary()
bc = DirichletBC(V, ubc, bd)
binds = bc.get_boundary_values().keys()
rows, cols, values = M.data()
M = sps.csc_matrix((values, cols, rows))
M = M[binds,:][:,binds]
rhs = b[binds]

x,info = sps.linalg.cg(M, rhs)
ubc.vector().zero()
ubc.vector()[binds] = x

# Solve
a = inner(grad(u), grad(v))*dx
L = Constant(0)*v*dx
u = Function(V)
solve(a==L, u, bc)
File("sol.pvd") << u

# compute boundary stress du/dn
a = inner(grad(u), grad(v))*dx
b = assemble(a-L)
b = b[binds]
tau,info = sps.linalg.cg(M, b)
stress = Function(V)
stress.vector().zero()
stress.vector()[binds] = tau
File("stress.pvd") << stress

So this method works. The stress is not accurate near the corners but we cannot expect anything better.

My question is more to do with the implementation details. I have to use scipy to construct the small mass matrix associated to boundary dofs. Because of this I had to use a CG solver from scipy to do the projection onto the boundary. Is it possible to use other solvers available in fenics ?

Answer 1

Hi, the smaller system can also be extracted with petsc4py. Consider

from petsc4py import PETSc
from dolfin import *

uexact = Expression('x[0]*x[0] - x[1]*x[1]')

n = 40
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, 'CG', 1)

u = TrialFunction(V)
v = TestFunction(V)
ubc = Function(V)

# Project boundary condition
M = PETScMatrix()
b = PETScVector()
assemble(u*v*ds, tensor=M)
assemble(uexact*v*ds, tensor=b)
bc = DirichletBC(V, ubc, DomainBoundary())
binds = bc.get_boundary_values().keys()
binds.sort()

# Extract submatrices
bdr_dofs = PETSc.IS()
bdr_dofs.createGeneral(binds)
M_petsc = PETSc.Mat()
b_petsc = PETSc.Vec()
x_petsc = PETSc.Vec()
M.mat().getSubMatrix(bdr_dofs, bdr_dofs, M_petsc)
b.vec().getSubVector(bdr_dofs, b_petsc)
b.vec().getSubVector(bdr_dofs, x_petsc)
M = PETScMatrix(M_petsc)
rhs = PETScVector(b_petsc)
x = PETScVector(x_petsc)
x.zero()
print M.size(0), M.size(1)

# Solve the projection system
solve(M, x, rhs)

# Assign x to boundary vector
ubc.vector()[binds] = x

# Solve
a = inner(grad(u), grad(v))*dx
L = Constant(0)*v*dx
u = Function(V)
solve(a==L, u, bc)

# compute boundary stress du/dn
a = inner(grad(u), grad(v))*dx
assemble(a-L, tensor=b)
b.vec().getSubVector(bdr_dofs, b_petsc)
b = PETScVector(b_petsc)
solve(M, x, b)
stress = Function(V)
stress.vector()[binds] = x

plot(u)
plot(stress)
interactive()

This approach has an advantage that it will let you convert the extracted petc4py.PETSc.Mat to DOLFIN's PETScMatrix, similarly for vectors, and then you can use all the solvers available in solve.

Comments

Thanks MiroK. I use mac binary which doesnt seem to have petsc4py. I will compile dolfin on my own and try your suggestion.