Hi everyone, I am solving a Nonlinear Darcy with porosity depending of the pressure. The boundary conditions are p=h in Gamma_D and u \cdot n = g in Gamma_N. The boundary condition for the pressure is natural and the velocity is essential. My problem is with the imposition of the Neumann condition in the space. I have tried to define the normal trace of the velocity like https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/mixed-poisson/python/documentation.html . On other hand, I have tried imposse the normal trace like a Dirichet condition using W.sub(0).sub(1) and doesn’t work either. By last, I have impossed the normal trace doing zero the term which is annulled in the exact velocity with the normal. In every case I haven’t got the correct solution. Is there any form to implement that?
from dolfin import *
import math
#Define Natural Boundary
class GammaD(SubDomain):
def inside(self, x, on_boundary):
return near(x[1],1.0) or near(x[0],1.0)
# Define essential boundary
class GammaNT(SubDomain):
def inside(self, x, on_boundary):
return near(x[1],0.0) #or near(x[1],1.0)
class GammaNU(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],0.0)
num = 4; h = []; errorp = []; erroru = []; errorpu = [];
ratep = []; rateu = []; ratepu = [];
rateu.append(0.0); ratep.append(0.0); ratepu.append(0.0);
Lx = 1.0; Ly = 1.0;
for dk in range(num):
d=0.2/pow(2,dk)
# Create mesh
mesh = RectangleMesh(Point(0.0, 0.0), Point(Lx, Ly), int(Lx/d), int(Ly/d),"crossed")
#Boundaries References: 1=Natural; 2=Essential
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(0)
gammaD = GammaD()
gammaD.mark(sub_domains, 1)
gammant = GammaNT()
gammant.mark(sub_domains, 2)
gammanu = GammaNU()
gammanu.mark(sub_domains, 3)
#Defining Measures
ds = Measure('ds', domain = mesh, subdomain_data = sub_domains)
n = FacetNormal(mesh)
hh = mesh.hmax()
h.append(hh)
# Define function spaces and mixed (product) space
V = VectorElement("P", mesh.ufl_cell(), 1)
Q = FiniteElement("P", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([V,Q]))
# Define trial and test functions
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
# Define source function
u1 = "pow(x[1],3"
u2 = "pow(x[0],3)"
u_ex=Expression(("-pow(x[1],3)","pow(x[0],3)"), domain=mesh, degree=3)
unt =Expression(("-pow(x[1],3)","0"), domain=mesh, degree=3)
unu =Expression(("0","-pow(x[0],3)"), domain=mesh, degree=3)
p_ex = Expression("pow(x[0],2)-pow(x[1],2)", domain=mesh, degree=2)
u3 = "-pow(x[1],3)"
u4 = "-pow(x[0],3)"
ut=Expression(u3, domain=mesh, degree=7)
uu=Expression(u4, domain=mesh, degree=7)
# Define variational form
f = u_ex-grad(p_ex)
B = (inner(u, v) + div(v)*p - div(u)*q -0.5*inner(u-grad(p),v+grad(q)))*dx +inner(div(u),div(v))*dx
F = inner(f,v)*dx + inner(dot(v,n), p_ex)*ds(1)-0.5*inner(f,v+grad(q))*dx
# Define function of normal trace like a dirichlet condition
class BoundarySource(UserExpression):
def __init__(self, mesh, **kwargs):
self.mesh = mesh
super().__init__(**kwargs)
def eval_cell(self, values, x, ufc_cell):
cell = Cell(self.mesh, ufc_cell.index)
n = cell.normal(ufc_cell.local_facet)
g = pow(x[0],3)
values[0] = g*n[0]
values[1] = g*n[1]
def value_shape(self):
return (2,)
G = BoundarySource(mesh, degree=2)
#bc1 = DirichletBC(W.sub(0), unt, sub_domains, 2)
#bc2 = DirichletBC(W.sub(0), unt, sub_domains, 3)
bc1 = DirichletBC(W.sub(0).sub(1), ut, sub_domains, 2)
bc2 = DirichletBC(W.sub(0).sub(1), uu, sub_domains, 3)
bc = [bc1,bc2]
# Compute solution
w = Function(W)
solve(B == F, w, bc)
(u, p) = w.split()
u.rename("velocity","velocity"); p.rename("pressure","pressure");
#Averages
pex_avg = Constant(assemble(p_ex*dx)/assemble(1.0*dx(domain=mesh)))
p_avg = Constant(assemble(p*dx)/assemble(1.0*dx(domain=mesh)))
#Error
ep = (grad(p_ex - p))**2*dx
ep = sqrt(abs(assemble(ep)))
errorp.append(ep)
e0 = (u_ex - u)**2*dx
ed = (div(u_ex - u))**2*dx
eu = sqrt(abs(assemble(e0+ed)))
erroru.append(eu)
error= ep+eu
errorpu.append(error)
#Logarithmic rates
if(dk>0):
ratep.append(ln(errorp[dk]/errorp[dk-1])/(ln(float(h[dk])/h[dk-1])))
if(dk>0):
rateu.append(ln(erroru[dk]/erroru[dk-1])/(ln(float(h[dk])/h[dk-1])))
if(dk>0):
ratepu.append(ln(errorpu[dk]/errorpu[dk-1])/(ln(float(h[dk])/h[dk-1])))
# ******** Generating table of convergences **** #
print( 'h & ||p-p_h||_1 & || u - u_h ||_0 & rate p & rate u & rate p-u')
print('=======================================================================')
for i in range(num):
print('%4.4g & %4.4g & %4.4g & %4.4g & %4.4g & %4.4g '% (h[i], errorp[i], erroru[i], ratep[i], rateu[i], ratepu[i]))
#Saving the solution
ufile_pvd = File("ejem2/velocity.pvd")
ufile_pvd << u
pfile_pvd = File("ejem2/pressure.pvd")
pfile_pvd << p
p_in=interpolate(p_ex, W.sub(1).collapse())
p_in.rename("exact-pressure","exact-pressure");
pinfile_pvd = File("ejem2/p_in.pvd")
pinfile_pvd << p_in
u_in=interpolate(u_ex, W.sub(0).collapse())
Greetings