Hello everyone,
I have a problem in the use of a mixed formulation of Poisson with the couple (BDM, DG) and with derivation operator.
So, I start with a mixed formulation of the Poisson equation but with “CG” spaces as experiment. To be more precise, I drive a flow from left to right with the Cls
p=1, left
p=0, right
u.n=0, top and bottom
I would now like to use a stable inf-sup pair and I saw that the tutorial recommended the pair (BDM, DG). So I adapted my scripts knowing that the condition u.n=0 is directly done by a Dirichlet condition (by the nature of the BDM space).
However, when compiling, I get a very strange behavior of the pressure and its derivatives.
After a plot, it seems to be linear along the x direction as it should be (i.e. y=-x) but by calculating its partial derivative along x (with the numerical operator “.dx(0)”) and integrating on the domain, i.e \int_\Omega \partial_x p~dx, I find a null result!
I should however find something close to -1.
For the first formulation, with the couple (CG,CG), I find 0.99999.
Note that by changing the DG to CG (with order 1), getting the (BDM,CG) couple, I get the expected result.
Do you have any idea to solve this problem?
Thank you
Here is a minimal running example:
from __future__ import print_function
from dolfin import *
from fenics import *
# Create classes for defining parts of the boundaries and the interior
# of the domain
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
class Obstacle(SubDomain):
def inside(self, x, on_boundary):
return (between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0)))
class Obstacle2(SubDomain):
def inside(self, x, on_boundary):
return (between(x[1], (0.1, 0.3)) and between(x[0], (0.7, 0.8)))
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
# Define mesh
mesh = UnitSquareMesh(70, 70)
# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
dx = Measure("dx", domain=mesh,subdomain_data=domains)
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
# Configuration with BDM
V = FiniteElement("BDM", mesh.ufl_cell(), 1)
Q = FiniteElement("DG", mesh.ufl_cell(), 0)
TH = V * Q
W = FunctionSpace(mesh, TH)
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
n = FacetNormal(mesh)
mu = 0.001
f = Constant(0.0)
K1 = as_matrix([[1, 0], [0,1]])
# Formulation avec BDM
bc_symx1 = DirichletBC(W.sub(0), (0.0,0.0), boundaries, 2)
bc_symx2 = DirichletBC(W.sub(0), (0.0,0.0), boundaries, 4)
bc_outflowS = DirichletBC(W.sub(1), 0.0, boundaries, 3)
bc1S = [bc_symx1, bc_symx2, bc_outflowS]
#Formulation du pb varia
a1S = mu*inner(inv(K1)*u,v)*dx\
- div(v)*p*dx - q*div(u)*dx + dot(v,n)*p*(ds(2)+ds(4))
L1S = -dot(v,n)*1*ds(1)
w1S = Function(W)
solve(a1S == L1S, w1S, bc1S)
(u1S, p1S) = w1S.split()
P11S = assemble((p1S.dx(0))*dx)
print(P11S)