Mixed problem pressure computation

Question

Hi,

I was trying to implement a mixed formulation solver for the Stokes equations as done in the FEniCS book on pp. 385-... Unfortunately I could not view the same results. The results for the velocity were quite good, while the ones of the pressure were not. In some cases the pressure seemed just to be displaced in "z"-direction although I had set up boundary conditions for the pressure as well (neither the "pinpoint" solution used in the book, nor setting the exact pressure solution on the whole boundary worked). Could you explain to me why the boundary conditions seem sometimes to be ignored completely?
Besides, how could I implement something like \int_\Omega p dx = 0 as restriction for my pressure function space ?

My minimal example is as follows:

from dolfin import *

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

class PinPoint(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < DOLFIN_EPS and x[1] < DOLFIN_EPS

n = 16
mesh = UnitSquareMesh(n, n, "crossed")

V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
W = V*Q

gammaD = DirichletBoundary()
u_ex = Expression(("sin(4*pi*x[0])*cos(4*pi*x[1])","-cos(4*pi*x[0])*sin(4*pi*x[1])"))
pinpoint = PinPoint()
p_ex = Expression("pi*cos(4*pi*x[0])*sin(4*pi*x[1])")

bc0 = DirichletBC(W.sub(0), u_ex, gammaD)
bc1 = DirichletBC(W.sub(1), p_ex, pinpoint, "pointwise")
bc = [bc0, bc1]

f = Expression(("28*pi*pi*sin(4*pi*x[0])*cos(4*pi*x[1])","-36*pi*pi*cos(4*pi*x[0])*sin(4*pi*x[1])"))

u_exInt = interpolate(u_ex, V)    
p_exInt = interpolate(p_ex, Q)


(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

a = inner(grad(u), grad(v))*dx-p*div(v)*dx+div(u)*q*dx
L = inner(f,v)*dx

w = Function(W)
solve(a == L, w, bc)

(u, p) = w.split()

print "e_u", errornorm(u_ex,u,norm_type="L2")
print "e_p", errornorm(p_ex,p,norm_type="L2")

print "e_u2 = %e" % assemble(inner(u_exInt-u,u_exInt-u)*dx)
print "e_p2 = %e" % assemble(inner(p_exInt-p,p_exInt-p)*dx)

plot(u, axes=True, title="uh")
plot(u_ex, mesh, axes=True, title="uex")
plot(p, axes=True, title="ph")
plot(p_ex, mesh, axes=True, title="pex")
interactive()

Comments

For giving at least an update: I decided to implement an additional Lagrange multiplier after having looked at the discrete problem but I can't give an analytical proof for the continuous case. So far it worked fine to change/add the following lines:

    ...

    V = VectorFunctionSpace(mesh, elemVeloc, degPolVeloc)
    Q = FunctionSpace(mesh, elemPress, degPolPressure)
    R = FunctionSpace(mesh, "R", 0)
    W = MixedFunctionSpace([V,Q,R])

   ...

    (u, p, c) = TrialFunctions(W)
    (v, q, d) = TestFunctions(W)

    a = inner(grad(u), grad(v))*dx-p*div(v)*dx+div(u)*q*dx+c*q*dx+p*d*dx
    L = inner(f,v)*dx

   ....

(elemVeloc, degPolVeloc, etc. are the corresponding finite elements specifications)