\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[1])*sin(pi*x[0])*sin(pi*x[0])*sin(pi*x[1])","-cos(pi*x[0])*sin(pi*x[1])*sin(pi*x[1])*sin(pi*x[0])"),domain=mesh, degree = 2)
p_e = Expression("sin(2*pi*x[0])*sin(2*pi*x[1])",domain=mesh,degree=2)
f = -div(grad(u))+grad(p)
bc = DirichletBC(W.sub(0), u_e , boundary)
a= inner(grad(u),grad(v))*dx-div(v)*p*dx+c0*q*dx+ div(u)*q*dx+p*c1*dx
L= dot(f,v)*dx
w = Function(W)
solve(a == L, w, bc)
(u, c) = w.split(deepcopy=True)

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).

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[1]*x[0]*x[0]*(1-x[0])*(1-x[0])*(1-x[1])*(2*x[1]-1)","2*x[0]*x[1]*x[1]*(1-x[1])*(1-x[1])*(1-x[0])*(1-2*x[0])"),domain=mesh, degree = 2)
p_e = Expression("x[0]*(1-x[0])*(1-x[1])-1/12",domain=mesh,degree=2)
#f = Expression(("4*x[1]*(x[1]-1)*(2*x[1]-1)*(6*x[0]*x[0]-6*x[0]+1)-x[0]*x[0]*(12-24*x[1])*(x[0]-1)*(x[0]-1)+(2*x[0]-1)*(1-x[1])","-x[1]*x[1]*(24*x[0]-12)*(x[1]-1)*(x[1]-1)-4*x[0]*(x[0]-1)*(2*x[0]-1)*(6*x[1]*x[1]-6*x[1]+1)-x[0]*(1-x[0])"),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
L= dot(-div(grad(u_e))+grad(p_e),v)*dx
w = Function(W)
solve(a == L, w, bc)
(u, p, c) = w.split()
error = dot(p-p_e,p-p_e) *dx
E = sqrt(abs(assemble(error)))
print(E)

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.