Manufactured solution method for multi-material

2

In your code, you are enforcing:
D_1\frac{\partial u}{\partial x} = D_2\frac{\partial u}{\partial x} on the interface between the two materials.
When D_1\neq D_2, this will cause a discontinuity in the gradient across the interface.
To counteract this behavior, you need to do your integration by parts thoroughly, and will end up with a term like this:

F += -2*pi*cos(2*pi*0.5)*(D_01*v("-")-D_02*v("+"))*dS(3)

where dS(3) is the interface between the two domains.

I’ve attached my adaptation of your code (where I avoid using sympy and only use dolfin imports:

from fenics import *

import matplotlib.pyplot as plt

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


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


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


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

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

left = Left()
right = Right()
left_bound = Left_bound()
right_bound = Right_bound()
interface = Interface()
mesh = UnitIntervalMesh(250)

volume_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
volume_markers.set_all(0)
left.mark(volume_markers, 1)
right.mark(volume_markers, 2)

surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
surface_markers.set_all(0)
left_bound.mark(surface_markers, 1)
right_bound.mark(surface_markers, 2)
interface.mark(surface_markers, 3)

V = FunctionSpace(mesh, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)

# Variational formulation

x= SpatialCoordinate(mesh)
u_D = 1 + sin(2*pi*x[0])  # exact solution

D_01, D_02 = Constant(2), Constant(5)  # diffusion coeffs
# source term left
f = grad(D_01*grad(u_D))[0,0]

# source term right
g = grad(D_02*grad(u_D))[0,0]


dx = Measure('dx', domain=mesh, subdomain_data=volume_markers)
dS = Measure("dS", domain=mesh, subdomain_data=surface_markers)

F = dot(D_01*grad(u), grad(v))*dx(1)
F += dot(D_02*grad(u), grad(v))*dx(2)

F += f*v*dx(1)
F += g*v*dx(2)
F += -2*pi*cos(2*pi*0.5)*(D_01*v("-")-D_02*v("+"))*dS(3)
# BCs
File("sf.pvd") << surface_markers
bcs = [
    DirichletBC(V, 1, surface_markers, 1),
    DirichletBC(V, 1, surface_markers, 2)]
# Solving
uh = Function(V)
solve(lhs(F)==rhs(F), uh, bcs=bcs)
plot(uh)
plt.savefig('out.png')
print(assemble(inner(uh-u_D,uh-u_D)*dx))