Laplace problem with discontinous parameter produces a big error at the discontinuity

Question

I am trying to solve the Laplace problem with a discontinous parameter. The problem is given by:
$$ \Omega_1 = { (x,y) \in [0,1]^{2} : y <= 1 - x} $$
$$ \Omega_1 = {(x,y) \in [0,1]^{2} : y > 1 - x} $$
$$- div(k*\nabla T) = 0$$
with $ k = 5 $ in $ \Omega_{1}$ and $ k = 1 $ in $ \Omega_{2}$.

By chosing appropriate DB-conditions, the problem has a analytical solution (see the code for the boundary conditions and the analytical solution).

My code is nearly the same as for the example from the Fenics tutorial, chapter "Handling domains with different materials", but when I examine the difference of the analytical solution and the numerical approximation, I have a huge error near the line $ y = 1 - x$ where the discontinuity of the parameter is.

Is there anything wrong with my code?
Thanks for helping me out.

Code:

from dolfin import *

nx = 64 
ny = 64
k1 = 5.0
k2 = 1.0
mesh = UnitSquareMesh(nx,ny,"right") 

# Define the two subdomains
class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] <= 1 - x[0] else False

class Omega2(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] > 1 - x[0] else False

omega1 = Omega1()
omega2 = Omega2()


# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
omega1.mark(domains, 0)
omega2.mark(domains, 1)


#plot(mesh, interactive = True)

class AnalyticalSolution(Expression):
    def eval(self, value, x):
        if x[1] <= 1 - x[0]:
            value[0] = ( -k1 * pow(x[0],2) + k1 * pow(x[1],2) + (1 - k2 + k1) * x[0] + ( k2 - k1 ) * x[1] )
        elif x[1] > 1 - x[0]:
            value[0] = ( -k2 * pow(x[0],2) + k2 * pow(x[1],2) + x[0] )

sol = AnalyticalSolution()

plot(sol, mesh, interactive  = True, wireframe = True)

# Boundaries
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1],0) and on_boundary

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1],1) and on_boundary

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0) and on_boundary

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],1) and on_boundary

bottom = Bottom()
top = Top()
left = Left()
right = Right()

# Functionspace for non-continous right hand side
V0 = FunctionSpace(mesh, "DG", 0)
k = Function(V0)

# Initialize coeffizient k
# values of k in the two subdomains
k_values = [k1, k2] 
for cell_no in range(len(domains.array())):
    subdomain_no = domains.array()[cell_no]
    k.vector()[cell_no] = k_values[subdomain_no]

plot(k, interactive = True, wireframe = True, axes = True)

# Functionspace for continous but not C1 solution
V1 = FunctionSpace(mesh, "CG", 1)

# Dirichlet boundary conditions
# bottom
dbbot = DirichletBC(V1, Expression("-a * pow(x[0],2) + (1 - b + a) * x[0]", a = k1, b = k2), bottom)
# top
dbtop = DirichletBC(V1, Expression("-b * pow(x[0],2) + b + x[0]", b = k2), top)
# left
dbleft = DirichletBC(V1, Expression("a * pow(x[1],2) + (b-a) * x[1]", a = k1, b = k2),  left)
# right
dbright = DirichletBC(V1, Expression("-b + b * pow(x[1],2) +1 ", b = k2), right)

db = [dbbot, dbtop, dbleft, dbright]

# Define trial- and testfunctions, homogeneous right hand side
u = TrialFunction(V1)
v = TestFunction(V1)
x = Function(V1)
f = Constant('0')

# Define variational forms
a = k * inner(nabla_grad(u),nabla_grad(v)) * dx
L = f * v * dx

problem = LinearVariationalProblem(a, L, x, db)
solver = LinearVariationalSolver(problem)
solver.solve()

plot(x, interactive = True, wireframe = True, axes = True)

Verr = FunctionSpace(mesh, "CG", 4)
isol = interpolate(sol, Verr)
ix = interpolate(x,Verr)
diff = project(isol-ix, Verr)
plot(diff, interactive = True, wireframe = True, axes = True)

Comments

As far as I can see your mesh does not align with the discontinuity of your coefficient. Have you tried

mesh = UnitSquareMesh(nx,ny,"left") 

or, alternatively,

mesh = UnitSquareMesh(nx,ny,"crossed") 

instead? Also, did you double-check that your exact solution is correct?

---

Thanks for your answer. The problem with this example is, that $ k \nabla T$ is not continuous along the line $y = 1 - x$. If I try to turn the dgl into its weak formulation, I have to consider integrals along the discontinuity line $y = 1 - x$ (which I didn't). So, even for a mesh where the discontinuity is along some edges of cells like for, there is a big error, because the variational formulation doesn't match the strong formulation.

My aim was to analzye how to implement weak discontinuities using fenics, when the discontinuity is not aligned with some edges of some cells of the mesh. If I use a correct test problem in $[0,1]^2$ with the discontinuity along $ y = 1 - x$ and use the mesh

mesh = UnitSquareMesh(nx,ny,"left")

(the discontinuity is along some edges), everything works fine and I get a solution with a small enough error. If I use

mesh = UnitSquareMesh(nx,ny,"right")

I get a big error.
I tried to use the PUM library and XFEM to solve this issue, but I guess PUM is only for solutions with strong discontinuities (only Heaviside-functions are used as the enrichment functions). So I guess PUM won't work for a problem like this.

I'm still watching for a solution for problems like this. So if you have any idea, let me know.