Poisson equation with multiple subdomains

This demo is implemented in a single Python file, demo_subdomains-poisson.py, which contains both the variational forms and the solver. We suggest that you familiarize yourself with the Poisson demo before studying this example, as some of the more standard steps will be described in less detail.

Implementation

This description goes through the implementation (in demo_subdomains-poisson.py) of a solver for the above described equation.

In this example, different boundary conditions are prescribed on different parts of the boundaries, and different parts of the interior have different material properties. This information must be made available to the solver. One way of doing this, is to tag the different subregions with different (integer) labels, and later integrate over the specified regions. DOLFIN provides a class MeshFunction which is useful for these types of operations: instances of this class represent functions over mesh entities (such as over cells or over facets). Mesh functions can be read from file or, if explicit formulae for the domains are known, they can be constructed by way of instances of the SubDomain class. The latter is the case here, so we begin by defining the left, right, top and bottom boundaries, and the interior obstacle domain using the SubDomain class and creating instances of these classes.

from dolfin import *

# Create classes for defining parts of the boundaries and the interior
# of the domain
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

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

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

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

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0)))

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle = Obstacle()

Note that the DOLFIN functions near and between provide robust ways of testing whether a coordinate is (to within machine precision) close to a given numerical value and in a range of values, respectively.

We next define a mesh of the domain:

mesh = UnitSquareMesh(64, 64)

The above subdomains are defined with the sole purpose of populating mesh functions. (For more complicated geometries, the mesh functions would typically be provided by other means.) We define a MeshFunction over the mesh cells to indicate which cells that correspond to the different interior subregions \(\Omega_0\) and \(\Omega_1\). Those in the interior rectangle will be tagged by 1, while the remainder is tagged by 0. We can set all the values of a MeshFunction to a given value using the set_all method. So, in order to accomplish what we want, we can set all values to 0 first, and then we can use the obstacle instance to mark the cells identified as inside the obstacle region by 1 (thus overwriting the previous value):

# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
obstacle.mark(domains, 1)

We can do the same for the boundaries using a MeshFunction defined ove the topological dimension ` mesh.topology().dim()-1`. We first tag all the edges by 0, then the edges on the left by 1, on the top by 2, on the right by 3 and on the bottom by 4:

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

Now that the geometry is defined and labeled, we can move on to defining the input source functions:

# Define input data
a0 = Constant(1.0)
a1 = Constant(0.01)
g_L = Expression("- 10*exp(- pow(x[1] - 0.5, 2))", degree=2)
g_R = Constant(1.0)
f = Constant(1.0)

Here, a0 and a1 represent the values of the coefficient \(a\) in the two regions of the domain, g_L and g_R represent the values of the Neumann boundary condition on the left and right boundaries respectively, and f represents the body source.

We may now move on to define the variational equation. As usual, we start by defining a finite element function space and basis functions on this space:

# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)

With this function space, we can define the essential (Dirichlet) boundary conditions on the top and bottom boundaries. These boundaries correspond to the facets tagged by 2 and 4, respectively, in the boundaries facet function:

# Define Dirichlet boundary conditions at top and bottom boundaries
bcs = [DirichletBC(V, 5.0, boundaries, 2),
       DirichletBC(V, 0.0, boundaries, 4)]

DOLFIN predefines the “measures” dx, ds and dS representing integration over cells, exterior facets (that is, facets on the boundary) and interior facets, respectively. These measures can take an additional integer argument. In fact, dx defaults to dx(0), ds defaults to ds(0), and dS defaults to dS(0). Integration over subregions can be specified by measures with different integer labels as arguments. However, we also need to map the geometry information stored in the mesh functions to these measures. The easiest way of accomplishing this is to define new measures with the mesh functions as additional input:

# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

We can now define the variational forms corresponding to the variational problem above using these measures and the tags for the different subregions. For simplicity, we define the full form first, and then extract the left- and right-hand sides using the UFL functions lhs() and rhs() afterwards. We can then solve as usual:

# Define variational form
F = (inner(a0*grad(u), grad(v))*dx(0) + inner(a1*grad(u), grad(v))*dx(1)
     - g_L*v*ds(1) - g_R*v*ds(3)
     - f*v*dx(0) - f*v*dx(1))

# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)

# Solve problem
u = Function(V)
solve(a == L, u, bcs)

Now, we can also evaluate various integrals of the solution or derived quantities of the solution over different regions, here are some examples:

# Evaluate integral of normal gradient over top boundary
n = FacetNormal(mesh)
m1 = dot(grad(u), n)*ds(2)
v1 = assemble(m1)
print("\int grad(u) * n ds(2) = ", v1)

# Evaluate integral of u over the obstacle
m2 = u*dx(1)
v2 = assemble(m2)
print("\int u dx(1) = ", v2)

We also plot the solution and its gradient:

# Plot solution and gradient
plot(u, title="u")
plot(grad(u), title="Projected grad(u)")
interactive()

Complete code

from dolfin import *

# Create classes for defining parts of the boundaries and the interior
# of the domain
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

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

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

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

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0)))

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle = Obstacle()

# Define mesh
mesh = UnitSquareMesh(64, 64)

# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
obstacle.mark(domains, 1)

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

# Define input data
a0 = Constant(1.0)
a1 = Constant(0.01)
g_L = Expression("- 10*exp(- pow(x[1] - 0.5, 2))", degree=2)
g_R = Constant(1.0)
f = Constant(1.0)

# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)

# Define Dirichlet boundary conditions at top and bottom boundaries
bcs = [DirichletBC(V, 5.0, boundaries, 2),
       DirichletBC(V, 0.0, boundaries, 4)]

# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Define variational form
F = (inner(a0*grad(u), grad(v))*dx(0) + inner(a1*grad(u), grad(v))*dx(1)
     - g_L*v*ds(1) - g_R*v*ds(3)
     - f*v*dx(0) - f*v*dx(1))

# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)

# Solve problem
u = Function(V)
solve(a == L, u, bcs)

# Evaluate integral of normal gradient over top boundary
n = FacetNormal(mesh)
m1 = dot(grad(u), n)*ds(2)
v1 = assemble(m1)
print("\int grad(u) * n ds(2) = ", v1)

# Evaluate integral of u over the obstacle
m2 = u*dx(1)
v2 = assemble(m2)
print("\int u dx(1) = ", v2)

# Plot solution
import matplotlib.pyplot as plt
plt.figure()
plot(u, title="Solution u")

# Plot solution and gradient
plt.figure()
plot(grad(u), title="Projected grad(u)")

# Show plots
plt.show()