Discontinuous coefficient in Poisson: Subdomain vs Expression

Question

Hello all,

I've been experimenting with the Poisson equation with a discontinuous coefficient, like in this example from the tutorial and the multiple subdomain demo.
There seem to be several ways to define the coefficients, and I wonder what the difference might be.
The method described in the demo is via subdomains.
Example:

from dolfin import *
import mshr

# Define Dirichlet boundary (x = 0 or x = 1)
class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] < DOLFIN_EPS and on_boundary

class BoundaryRight(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] > 1 - DOLFIN_EPS and on_boundary

class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] <= 0.5 else False

class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] >= 0.5 else False

domain = mshr.Rectangle(Point(0., 0.), Point(1., 1.))
mesh = mshr.generate_mesh(domain, 32)

V = FunctionSpace(mesh, "CG", 1)

domains = CellFunction("size_t", mesh)
domains.set_all(0)
subdomain0 = Omega0()
subdomain0.mark(domains, 0)
subdomain1 = Omega1()
subdomain1.mark(domains, 1)
dx = Measure('dx', domain=mesh, subdomain_data=domains)

left = BoundaryLeft()
right = BoundaryRight()
bcs = [DirichletBC(V, 0., left),
       DirichletBC(V, 1., right)]

k0 = Constant(1.5)
k1 = Constant(50.)
f = Constant(0.)

# Exact solution
u_ex = Expression("x[1] <= 0.5 ? 2*x[1]*k1/(k0+k1) : \
        ((2*x[1]-1)*k0+k1)/(k0+k1)", k0=float(k0), k1=float(k1))

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

L = f*v*dx
a = inner(k0*grad(u), grad(v))*dx(0) + inner(k1*grad(u), grad(v))*dx(1)
u1 = Function(V)
solve(a == L, u1, bcs)

print("Err = %.3e" % errornorm(u_ex, u1))

Another method I consider is by interpolating an expression onto a DG-0 space, instead of using subdomains:

dx = Measure('dx')
kf = Expression("x[1]<=0.5 ? k0 : k1", k0=k0, k1=k1)
V0 = FunctionSpace(mesh, "DG", 0)
ki = interpolate(kf, V0)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(ki*grad(u), grad(v))*dx
u2 = Function(V)
solve(a == L, u2, bcs)

print("Err = %.3e" % errornorm(u_ex, u2))

On the "unstructured" mesh, the error of the second method appears to be generally smaller than with the first one. With UnitSquareMesh(32, 32) the problem is solved exactly and in both cases the error is exactly the same and near machine precision; but this is not the case I'm interested in.

Apart from the general advantages that subdomains may have for complex externally created meshs over expressions, what is the recommended way? What is the technical difference?

Thank you for your help!
David

Comments

I wanted to ask very much the same question one of these days, so thanks for beating me to it

---

yes.. would be great to have some insights =)

Answer 1

The solutions are different because the coefficients are different.

kf = Expression("x[1]<=0.5 ? k0 : k1", k0=k0, k1=k1)
V0 = FunctionSpace(mesh, "DG", 0)
ki = interpolate(kf, V0)

This results in a coefficient that equals k0 on the cells where the midpoint x1 coordinate is lesser or equal to 0.5, and k1 on the remaining cells.

subdomain1.mark(domains, 1)
dx = Measure('dx', domain=mesh, subdomain_data=domains)

The cells marked 1 are those that have all three vertices and the midpoint inside subdomain1.

The discontinuous coefficient can also be implemented a UFL conditional:

x0, x1 = MeshCoordinates(mesh)
k = conditional(x1 <= 0.5, k0, k1)

This avoids discretising the discontinous diffusion coefficient. The coefficient will be correctly evaluated at each quadrature point (but the quadrature rule will be inexact).

The last method is arguably better than the other the two since it admits optimal convergence rates for sufficiently smooth solutions when the polynomial degree of the FEM space is larger than 1. (This does not result in better a priori estimates since diffusion problems with discontinuous diffusion coefficient are not in general sufficiently smooth.)

Comments

Great, thanks for the answer!

---

In addition to Magne's excellent answer, with a Quadrature function space ("Q") you should get the same answer as with the UFL conditional approach. This will evaluate your exact coefficient at each quadrature point.

---

Ok, thanks. So summarizing, if a formula for the 'domain partition' is known, use the UFL conditional or an Expression interpolated to the "Q"-space; if not, e.g., for an externally generated mesh, you're left with the (less optimal) subdomain.mark() method.

---

It will be inefficient computationally and programmatically, but you could make two "Q"-based functions, then use (discrete) subdomain marking to get ds(1) and ds(2) and repeat your bilinear form term twice. Then you will get your interface.

I believe restricting FunctionSpaces to subsets of meshes is in the works, then this will be both efficient and easy.

---

Oh, sounds good, I guess being able to restrict FunctionSpaces to subsets of meshes would resolve this other problem I'm having. Hehe.

---

Is this preferred 3rd method still valid? When I use 1D and put
x0 = MeshCoordinates(mesh)
k = conditional(x0 <= 0.5, k0, k1)
then I get an error message:
"Expecting scalar arguments."

---

Does it work if you index into x0, e.g. x0[0]?

---

Use SpatialCoordinate. For example,

x = SpatialCoordinate(mesh)
k = conditional(x[0] <= 0.5, k0, k1)

---

Not that way, Jack (MeshCoordinates with x0[0]). It segfaults with

#0 0x00007fffe7ef9a33 in dolfin::MeshCoordinates::MeshCoordinates(std::shared_ptr) () from /usr/lib/x86_64-linux-gnu/libdolfin.so.2016.1
#1 0x00007fffbf7b4a3a in ?? () from /usr/lib/python2.7/dist-packages/dolfin/cpp/_function.x86_64-linux-gnu.so
#2 0x0000555555654e70 in PyEval_EvalFrameEx ()
#3 0x000055555564d3c5 in PyEval_EvalCodeEx ()

---

Hmm, no nested replies, so it looks like I'm replying to myself...

Thanks Magne, your update with SpatialCoordinate works great.