Hello everybody,
I want to use dolfin to solve a problem in a porous media. For now, I consider the following problem:
Let \Omega=\Omega_1 \cup \Omega_2\subset\mathbb{R}^2 be a domain consisting of two parts. Indicate with \Gamma the interface between the two parts. Then I want to solve the PDE:
\begin{align}
-\Delta u_i &= f_i \text{, in } \Omega_i \\
\vec{n}\nabla u_i &= r(u_1 - u_2) \text{, on } \Gamma
\end{align}
Where r\in\mathbb{R} is just some parameter. Essentially a simple Poisson problem with a Robin-condition at the interface (which leads to discontinuity). Where my domain is given by a rectangle (\Omega_1) with a circle (\Omega_2) in the center.
The variational formulation is then given with a jump condition at the interface.
To solve this in FEniCS I found, until now, two different ways. Either use DG-Elements or construct a two-dimensional function space, where the first entry corresponds to the solution in \Omega_1 and the second to \Omega_2. For the second idea, one has to eliminate the DOF that are not used with .ident_zeros(). See for example here. There they have essentially the same problem but on a simpler domain.
Right now I am using the idea with the function space.
The problem I have is that the solution has not the correct behavior at the interface. The function u_2, that should only exist inside the circle, is correct and exists only \Omega_2 and all DOFs in \Omega_1 are set to 0.
But if you look at u_1 ones sees that is is not correct and also partially exists inside \Omega_2. But strangely only at the vertices right behind the interface.
Here my Code with mshr:
from dolfin import *
from mshr import *
fileU = File("test_pore_2D/.pvd")
size = 1
center = [size*0.5, size*0.5]
radius = 0.25
domain = Rectangle(Point(0, 0), Point(size, size))
circle = Circle(Point(center[0], center[1]), radius)
domain.set_subdomain(1, circle)
pore_mesh = generate_mesh(domain, 20)
cell_markers = MeshFunction("size_t", pore_mesh, 2, mesh.domains())
facet_markers = MeshFunction("size_t", pore_mesh, 1)
# Use the cell domains associated with each facet to set the boundary
for f in facets(pore_mesh):
domains = []
for c in cells(f):
domains.append(cell_markers[c])
# check that facet has adjoint with different markers
if len(domains) == 2 and domains[0] != domains[1]:
facet_markers[f] = 1
# define the measures
dS = Measure("dS", domain=pore_mesh, subdomain_data=facet_markers)
dx = Measure('dx', domain=pore_mesh, subdomain_data=cell_markers)
# define variational forms
Element1 = FiniteElement("Lagrange", pore_mesh.ufl_cell(), 1)
W_elem = MixedElement([Element1, Element1])
W = FunctionSpace(pore_mesh, W_elem)
u = TrialFunction(W)
v = TestFunction(W)
# dirichlet condition for the outer part
def dirich_boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(W.sub(0), (0.0), dirich_boundary)
f2 = Constant(150.0)
f1 = Constant(-150.0)
r = Constant(3.0)
# cell integrals
a = inner(grad(u[0]), grad(v[0])) * dx(0) + inner(grad(u[1]), grad(v[1])) * dx(1)
L = f1 * v[0] * dx(0) + f2 * v[1] * dx(1)
# interface integral
jmp_u = avg(u[0]) - avg(u[1])
jmp_v = avg(v[0]) - avg(v[1])
a += r * jmp_u * jmp_v * dS(1)
# assemble and eliminate excess dofs
A, b = assemble_system(a, L, bc)
A.ident_zeros()
# solve
w = Function(W)
solve(A, w.vector(), b)
fileU << w
If I remove the jump condition at the interface, u_1 will be correct (u_2 then not well-defined anymore, but still only exist on \Omega_2). So I think the error is connected to the interface condition. But checking assemble(1dS(1)), assemble(1dx(1)) or assemble(1dx(0))* gives the correct values.
Maybe this is connected to the restriction on the interface. Where avg(u[0]) is used, but changing this to the restriction (’+’) or (’-’) does not help.
If I use the example in the code above, this does not happen. There, mshr is not used, but the cells are marked after the creation of a simple mesh. But this is not possible for me, because of the inner circle. Creating new markers afterwards sadly does not solve the problem.
My last idea was to create the mesh with gmsh and then import it. Sadly, this leads to the same problems. Here would be my gmsh code:
SetFactory("OpenCASCADE");
//+
Point(1) = {0, 0, 0, 0.2};
//+
Point(2) = {1, 0, 0, 0.2};
//+
Point(3) = {1, 1, 0, 0.2};
//+
Point(4) = {0, 1, 0, 0.2};
//+
Line(1) = {1, 4};
//+
Line(2) = {4, 3};
//+
Line(3) = {3, 2};
//+
Line(4) = {2, 1};
//+
Point(5) = {0.75, 0.5, -0, 0.1};
//+
Point(6) = {0.25, 0.5, -0, 0.1};
//+
Point(7) = {0.5, 0.5, -0, 0.2};
//+
Circle(5) = {5, 7, 6};
//+
Circle(6) = {6, 7, 5};
//+
Curve Loop(1) = {2, 3, 4, 1};
//+
Curve Loop(2) = {6, 5};
//+
Plane Surface(1) = {1, 2};
//+
Curve Loop(3) = {6, 5};
//+
Plane Surface(2) = {3};
//+
Physical Surface("inner", 2) = {2};
//+
Physical Surface("outer", 1) = {1};
And one would have to change the first lines and the numbering of the measures to use the gmsh mesh. (dx(0) → dx(1) …) And loading it with:
pore_mesh = Mesh()
# load the whole mesh
with XDMFFile("pore_mesh/mesh.xdmf") as infile:
infile.read(pore_mesh)
# load the cell domain-markers (1=outer, 2=inner)
mvc = MeshValueCollection("size_t", pore_mesh, 2)
with XDMFFile("pore_mesh/cell_marker.xdmf") as infile:
infile.read(mvc, "name_to_read")
cell_markers = cpp.mesh.MeshFunctionSizet(pore_mesh, mvc)
Maybe one could suspect that this has no negative contribution to the solution. But to check this, I also tried to solve the problem in the link from above with mshr. There, one also will get this incorrect behavior of u_1 at the interface and no convergence.
Any help would be appreciated, thank you!