Subdomains and boundary conditions

So far, we have only looked briefly at how to specify boundary conditions. In this chapter, we look more closely at how to specify boundary conditions on specific parts (subdomains) of the boundary and how to combine multiple boundary conditions. We will also look at how to generate meshes with subdomains and how to define coefficients with different values in different subdomains.

Combining Dirichlet and Neumann conditions

Let's return to the Poisson problem from the chapter Fundamentals: Solving the Poisson equation and see how to extend the mathematics and the implementation to handle a Dirichlet condition in combination with a Neumann condition. The domain is still the unit square, but now we set the Dirichlet condition $ u=\ub $ at the left and right sides, $ x=0 $ and $ x=1 $, while the Neumann condition $$ \begin{equation*} -{\partial u\over\partial n}=g \end{equation*} $$ is applied to the remaining sides $ y=0 $ and $ y=1 $.

PDE problem

Let $ \GD $ and $ \GN $ denote the parts of the boundary $ \partial\Omega $ where the Dirichlet and Neumann conditions apply, respectively. The complete boundary-value problem can be written as $$ \begin{alignat}{2} - \nabla^2 u &= f \quad&&\mbox{in } \Omega, \\ u &= \ub &&\mbox{on } \GD, \\ - {\partial u\over\partial n} &= g &&\mbox{on } \GN \tp \end{alignat} $$ Again, we choose $ u=1+x^2 + 2y^2 $ as the exact solution and adjust $ f $, $ g $, and $ \ub $ accordingly: $$ \begin{align*} f(x, y) &= -6,\\ g(x, y) &= \left\lbrace\begin{array}{ll} 0, & \quad y=0,\\ -4, & \quad y=1, \end{array}\right.\\ \ub(x, y) &= 1 + x^2 + 2y^2\tp \end{align*} $$

For ease of programming, we define $ g $ as a function over the whole domain $ \Omega $ such that $ g $ takes on the correct values at $ y=0 $ and $ y=1 $. One possible extension is $$ \begin{equation*} g(x,y) = -4y\tp \end{equation*} $$

Variational formulation

The first task is to derive the variational formulation. This time we cannot omit the boundary term arising from the integration by parts, because $ v $ is only zero on $ \GD $. We have $$ \begin{equation*} -\int_\Omega (\nabla^2 u)v \dx = \int_\Omega\nabla u\cdot\nabla v \dx - \int_{\partial\Omega}{\partial u\over \partial n}v \ds, \end{equation*} $$ and since $ v=0 $ on $ \GD $, $$ \begin{equation*} - \int_{\partial\Omega}{\partial u\over \partial n}v \ds = - \int_{\GN}{\partial u\over \partial n}v \ds = \int_{\GN}gv \ds, \end{equation*} $$ by applying the boundary condition on $ \GN $. The resulting weak form reads $$ \begin{equation} \int_{\Omega} \nabla u \cdot \nabla v \dx = \int_{\Omega} fv \dx - \int_{\GN} gv \ds\tp \tag{4.1} \end{equation} $$ Expressing this equation in the standard notation $ a(u,v)=L(v) $ is straightforward with $$ \begin{align} a(u, v) &= \int_{\Omega} \nabla u \cdot \nabla v \dx, \tag{4.2}\\ L(v) &= \int_{\Omega} fv \dx - \int_{\GN} gv \ds\tp \tag{4.3} \end{align} $$

FEniCS implementation

How does the Neumann condition impact the implementation? Let us revisit our previous implementation ft01_poisson.py from the section FEniCS implementation and examine which changes we need to make to incorporate the Neumann condition. It turns out that only two changes are necessary:

The function boundary defining the Dirichlet boundary must be modified.
The new boundary term must be added to the expression for L.

The first adjustment can be coded as

tol = 1E-14

def boundary_D(x, on_boundary):
    if on_boundary:
        if near(x[0], 0, tol) or near(x[0], 1, tol):
            return True
        else:
            return False
    else:
        return False

A more compact implementation reads

def boundary_D(x, on_boundary):
    return on_boundary and (near(x[0], 0, tol) or near(x[0], 1, tol))

The second adjustment of our program concerns the definition of L, which needs to include the Neumann condition:

g = Expression('-4*x[1]', degree=1)
L = f*v*dx - g*v*ds

The ds variable implies a boundary integral, while dx implies an integral over the domain $ \Omega $. No other modifications are necessary.

Note that the integration *ds is carried out over the entire boundary, including the Dirichlet boundary. However, since the test function v vanishes on the Dirichlet boundary (as a result specifying a DirichletBC), the integral will only include the contribution from the Neumann boundary.