# Variational form of stokes problem with homogeneous dirichlet boundary condition

Sir/Mam,
The Model problem is:

• \Delta(u)+\grad(p)=f in \Omega
div(u)=0 in \Omega
u=0 on \del_\Omega
\int_\Omega p dx=0

I impose the average zero condition of p with the help of lagrange multiplier. Weak formulation of my code is showing error.

the code is given below:

from dolfin import *

#Create mesh and define function space

Nx=64
mesh = UnitSquareMesh(Nx,Nx)
V = VectorElement("CG", mesh.ufl_cell(), 1)
F = FiniteElement("CG", mesh.ufl_cell(), 1)
R = VectorElement("Real",mesh.ufl_cell(), 0)
# Make a mixed space
me = MixedElement([V,F,R])
W = FunctionSpace(mesh,me)

# Define boundary condition

def boundary(x, on_boundary):
return on_boundary

# Define test and trial function
(u,p,c0) = TrialFunctions(W)
(v,q,c1) = TestFunctions(W)

u_e = Expression(("cos(pi*x)*sin(pi*x)*sin(pi*x)*sin(pi*x)","-cos(pi*x)*sin(pi*x)*sin(pi*x)*sin(pi*x)"),domain=mesh, degree = 2)
p_e = Expression("sin(2*pi*x)*sin(2*pi*x)",domain=mesh,degree=2)
bc = DirichletBC(W.sub(0), u_e , boundary)
L= dot(f,v)*dx
w = Function(W)
solve(a == L, w, bc)
(u, c) = w.split(deepcopy=True)


Here are a few notes:

• R should be a scalar FiniteElement, not a VectorElement.
• To use the method of manufactured solutions, you want to plug the exact solution (\mathbf{u}_\text{e},p_\text{e}) into the strong form to get \mathbf{f}, not the TrialFunctions used to define the left-hand-side form.
• The call to the split method at the end needs a three-tuple to assign its output to.
• The choices of equal-order pressure and velocity spaces are not inf-sup stable, so you would either want to use a different combination (e.g., the Taylor–Hood element, where V would has one higher polynomial degree than F) or a stabilized formulation (e.g., the PSPG method).
1 Like

can you explain your second point more clearly by writing its weak formulation?

In terms of your code, the second point would amount to using

f = -div(grad(u_e))+grad(p_e)


as the source term instead of

f = -div(grad(u))+grad(p)


as i already apply the mean zero condition of p in weak formulation but it provide almost same error as the mesh is more and more finer. I dont understand what is the exact error. please help me!


from fenics import *

Nx=16
mesh = UnitSquareMesh(Nx,Nx)
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
F = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
R = FiniteElement("Real",mesh.ufl_cell(), 0)
me = MixedElement([V,F,R])
W = FunctionSpace(mesh,me)

def boundary(x, on_boundary):
return on_boundary

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

u_e = Expression(("2*x*x*x*(1-x)*(1-x)*(1-x)*(2*x-1)","2*x*x*x*(1-x)*(1-x)*(1-x)*(1-2*x)"),domain=mesh, degree = 2)
p_e = Expression("x*(1-x)*(1-x)-1/12",domain=mesh,degree=2)
#f = Expression(("4*x*(x-1)*(2*x-1)*(6*x*x-6*x+1)-x*x*(12-24*x)*(x-1)*(x-1)+(2*x-1)*(1-x)","-x*x*(24*x-12)*(x-1)*(x-1)-4*x*(x-1)*(2*x-1)*(6*x*x-6*x+1)-x*(1-x)"),degree=2)

bc = DirichletBC(W.sub(0), u_e , boundary)

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

The problem here is that the string passed to Expression in legacy FEniCS is actually a fragment of C++ code, so 1/12 is being evaluated using integer division, and truncates to zero. (This was what would happen in Python 2 as well, but in Python 3, / is floating-point division, while // is integer division.) This potential pitfall is avoided in standard usage of the new FEniCSx, which is now recommended for general use in new projects, but you can fix your current code by replacing 1/12 with 1.0/12.0 in the C++ code.