# 19. 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.

The main purpose of this demo is to demonstrate how to create and integrate variational forms over distinct regions of a domain and/or its boundaries.

## 19.1. Equation and problem definition¶

For illustration purposes, we consider a weighted Poisson equation over a unit square: \(\Omega = [0,1] \times [0,1]\) with mixed boundary conditions: find \(u\) satisfying

where \(\Gamma_{T}\), \(\Gamma_{B}\), \(\Gamma_{L}\), \(\Gamma_{R}\) denote the top, bottom, left and right sides of the unit square, respectively. The coefficient \(a\) may vary over the domain: here, we let \(a = a_1 = 0.01\) for \((x, y) \in \Omega_1 = [0.5, 0.7] \times [0.2, 1.0]\) and \(a = a_0 = 1.0\) in \(\Omega_0 = \Omega \backslash \Omega_1\). We can think of \(\Omega_1\) as an obstacle with differing material properties from the rest of the domain.

## 19.2. Variational form¶

We can write the above boundary value problem in the standard linear variational form: find \(u \in V\) such that

where \(V\) and \(\hat{V}\) are suitable function spaces incorporating the Dirichlet boundary conditions on \(\Gamma_{T}\) and \(\Gamma_{B}\), and

where \(f = 1.0\), \(g_L = - 10.0 e^{-(y - 0.5)^2}\) and \(g_R = 1.0\).

## 19.3. 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.) The classes
`CellFunction`

and
`FacetFunction`

are specialized
versions of the more general `MeshFunction`

. `CellFunction`

represents a function with a value for each
cell of a mesh, while `FacetFunction`

represents a function with a value for
each facet. We define a `CellFunction`

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 = CellFunction("size_t", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)
```

We can do the same for the boundaries using a `FacetFunction`

. 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 = FacetFunction("size_t", mesh)
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()
```

## 19.4. Complete code¶

```
from __future__ import print_function
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 = CellFunction("size_t", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
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 and gradient
plot(u, title="u")
plot(grad(u), title="Projected grad(u)")
interactive()
```