$$ \newcommand{\dt}{\Delta t} \newcommand{\tp}{\thinspace .} \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\x}{\boldsymbol{x}} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\uI}{u_{_0}} \newcommand{\ub}{u_{_\mathrm{D}}} \newcommand{\GD}{\Gamma_{_\mathrm{D}}} \newcommand{\GN}{\Gamma_{_\mathrm{N}}} \newcommand{\GR}{\Gamma_{_\mathrm{R}}} \newcommand{\inner}[2]{\langle #1, #2 \rangle} $$




Setting multiple Dirichlet, Neumann, and Robin conditions

Consider again the variable-coefficient Poisson problem from the section Defining subdomains for different materials. We will now discuss how to implement general combinations of boundary conditions of Dirichlet, Neumann, and Robin type for this model problem.

Three types of boundary conditions

We extend our repertoire of boundary conditions to three types: Dirichlet, Neumann, and Robin. Dirichlet conditions apply to some parts \( \GD^0, \GD^1, \ldots \), of the boundary: $$ u = \ub^0\hbox{ on }\GD^0,\quad u = \ub^1\hbox{ on }\GD^1, \quad \ldots,$$ where \( \ub^i \) are prescribed functions, \( i=0,1,\ldots \). On other parts, \( \GN^0, \GN^1, \ldots \), we have Neumann conditions: $$ -\kappa{\partial u\over\partial n} = g_{0}\hbox{ on }\GN^0,\quad -\kappa{\partial u\over\partial n} = g_{1}\hbox{ on }\GN^1,\quad \ldots, $$ Finally, we have Robin conditions: $$ \begin{equation*} -\kappa{\partial u\over\partial n} = r(u-s), \tag{4.5} \end{equation*} $$ where \( r \) and \( s \) are specified functions. The Robin condition is most often used to model heat transfer to the surroundings and arise naturally from Newton's cooling law. In that case, \( r \) is a heat transfer coefficient, and \( s \) is the temperature of the surroundings. Both can be space and time-dependent. The Robin conditions apply at some parts \( \GR^0, \GR^1, \ldots \), of the boundary: $$ -\kappa{\partial u\over\partial n} = r_0(u-s_0)\hbox{ on }\GR^0,\quad -\kappa{\partial u\over\partial n} = r_1(u-s_1)\hbox{ on }\GR^1,\quad \ldots $$

PDE problem

With the notation above, the model problem to be solved with multiple Dirichlet, Neumann, and Robin conditions can be formulated as follows: $$ \begin{alignat}{2} -\nabla\cdot(\kappa\nabla u) &= f \quad&&\mbox{in } \Omega, \tag{4.6}\\ u &= \ub^i &&\mbox{on } \GD^i,\quad i=0,1,\ldots \tag{4.7}\\ -\kappa{\partial u\over\partial n} &= g_i &&\mbox{on } \GN^i,\quad i=0,1,\ldots \tag{4.8}\\ -\kappa{\partial u\over\partial n} &= r_i(u-s_i) \quad&&\mbox{on } \GR^i,\quad i=0,1,\ldots \tag{4.9} \end{alignat} $$

Variational formulation

As usual, we multiply by a test function \( v \) and integrate by parts: $$ \begin{equation*} -\int_\Omega \nabla\cdot(\kappa\nabla u) v \dx = \int_\Omega \kappa\nabla u\cdot \nabla v \dx - \int_{\partial\Omega}\kappa\frac{\partial u}{\partial n}v \ds\tp \end{equation*} $$ On the Dirichlet part of the boundary (\( \GD^i \)), the boundary integral vanishes since \( v = 0 \). On the remaining part of the boundary, we split the boundary integral into contributions from the Neumann parts (\( \GN^i \)) and Robin parts (\( \GR^i \)). Inserting the boundary conditions, we obtain $$ \begin{align*} -\int_{\partial\Omega} \kappa\frac{\partial u}{\partial n}v \ds &= -\sum_i\int_{\GN^i} \kappa\frac{\partial u}{\partial n} \ds -\sum_i\int_{\GR^i} \kappa\frac{\partial u}{\partial n} \ds\\ &= \sum_i\int_{\GN^i}g_i \ds + \sum_i\int_{\GR^i}r_i(u-s_i) \ds\tp \end{align*} $$ We thus obtain the following variational problem: $$ \begin{equation} F = \int_{\Omega} \kappa\nabla u\cdot \nabla v \dx + \sum_i\int_{\GN^i} g_iv \ds + \sum_i\int_{\GR^i}r_i(u-s_i)v \ds - \int_{\Omega} fv \dx =0\tp \tag{4.10} \end{equation} $$

We have been used to writing this variational formulation in the standard notation \( a(u,v)=L(v) \), which requires that we identify all integral depending on the trial function \( u \), and collect these in \( a(u,v) \), while the remaining integrals go into \( L(v) \). The integrals from the Robin condition must for this reason be split into two parts: $$ \begin{equation*} \int_{\GR^i}r_i(u-s_i)v \ds = \int_{\GR^i} r_iuv \ds - \int_{\GR^i}r_is_iv \ds\tp \end{equation*} $$ We then have $$ \begin{align} a(u, v) &= \int_{\Omega} \kappa\nabla u\cdot \nabla v \dx + \sum_i\int_{\GR^i}r_iuv \ds, \tag{4.11}\\ L(v) &= \int_{\Omega} fv \dx - \sum_i\int_{\GN^i} g_i v \ds + \sum_i\int_{\GR^i}r_is_iv \ds\tp \tag{4.12} \end{align} $$ Alternatively, we may keep the formulation (4.10) and either solve the variational problem as a nonlinear problem (F == 0) in FEniCS or use the FEniCS functions lhs and rhs to extract the bilinear and linear parts of F:

a = lhs(F)
L = rhs(F)

Note that if we choose to solve this linear problem as a nonlinear problem, the Newton iteration will converge in a single iteration.

FEniCS implementation

Let us examine how to extend our Poisson solver to handle general combinations of Dirichlet, Neumann, and Robin boundary conditions. Compared to our previous code, we must consider the following extensions:

A general approach to the first task is to mark each of the desired boundary parts with markers 0, 1, 2, and so forth. Here we aim at the four sides of the unit square, marked with 0 (\( x=0 \)), 1 (\( x=1 \)), 2 (\( y=0 \)), and 3 (\( y=1 \)). The markers will be defined using a MeshFunction, but contrary to the section Defining subdomains for different materials, this is not a function over cells, but a function over the facets of the mesh. We use a FacetFunction for this purpose:

boundary_markers = FacetFunction('size_t', mesh)

As in the section Defining subdomains for different materials we use a subclass of SubDomain to identify the various parts of the mesh function. Problems with domains of more complicated geometries may set the mesh function for marking boundaries as part of the mesh generation. In our case, the boundary \( x = 0 \) can be marked as follows:

class BoundaryX0(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0, tol)

bx0 = BoundaryX0()
bx0.mark(boundary_markers, 0)

Similarly, we create the classes BoundaryX1 (\( x=1 \)), BoundaryY0 (\( y=0 \)), and BoundaryY1 (\( y=1 \)) boundary, and mark these as subdomains 1, 2, and 3, respectively.

For generality of the implementation, we let the user specify what kind of boundary condition that applies to each of the four boundaries. We set up a Python dictionary for this purpose, with the key as subdomain number and the value as a dictionary specifying the kind of condition as key and a function as its value. For example,

boundary_conditions = {0: {'Dirichlet': u_D},
                       1: {'Robin':     (r, s)},
                       2: {'Neumann':   g},
                       3: {'Neumann',   0}}


As explained in the section Setting multiple Dirichlet conditions, multiple Dirichlet conditions must be collected in a list of DirichletBC objects. Based on the boundary_conditions data structure above, we can construct this list by the following code snippet:

bcs = []
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
                         boundary_markers, i)

A new aspect of the variational problem is the two distinct boundary integrals over \( \GN^i \) and \( \GR^i \). Having a mesh function over exterior cell facets (our boundary_markers object), where subdomains (boundary parts) are numbered as \( 0,1,2,\ldots \), the special symbol ds(0) implies integration over subdomain (part) 0, ds(1) denotes integration over subdomain (part) 1, and so on. The idea of multiple ds-type objects generalizes to volume integrals too: dx(0), dx(1), etc., are used to integrate over subdomain 0, 1, etc., inside \( \Omega \).

To express integrals over the boundary parts using ds(i), we must first redefine the measure ds in terms of our boundary markers:

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

Similarly, if we want integration over different parts of the domain, we redefine dx as

dx = Measure('dx', domain=mesh, subdomain_data=domain_markers)

where domain_markers is a CellFunction defining subdomains in \( \Omega \).

Suppose we have a Robin condition with values r and s on subdomain R, and a Neumann condition with value g on subdomain N. The variational form can then be written

a = kappa*dot(grad(u), grad(v))*dx + r*u*v*ds(R)
L = f*v*dx - g*v*ds(N) + r*s*v*ds(R)

In our case, things get a bit more complicated since the information about integrals in Neumann and Robin conditions are in the boundary_conditions data structure. We can collect all Neumann conditions by the following code snippet:

integrals_N = []
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        if boundary_conditions[i]['Neumann'] != 0:
            g = boundary_conditions[i]['Neumann']

Applying sum(integrals_N) will apply the + operator to the variational forms in the integrals_N list and result in the integrals we need for the right-hand side L of the variational form.

The integrals in the Robin condition can similarly be collected in lists:

integrals_R_a = []
integrals_R_L = []
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        r, s = boundary_conditions[i]['Robin']

We are now in a position to define the a and L expressions in the variational formulation:

a = kappa*dot(grad(u), grad(v))*dx + sum(integrals_R_a)
L = f*v*dx - sum(integrals_N) + sum(integrals_R_L)

Alternatively, we may use the FEniCS functions lhs and rhs as mentioned above to simplify the extraction of terms for the Robin integrals:

integrals_R = []
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        r, s = boundary_conditions[i]['Robin']
        integrals_R.append(r*(u - s)*v*ds(i))

F = kappa*dot(grad(u), grad(v))*dx + \ 
    sum(integrals_R) - f*v*dx + sum(integrals_N)
a, L = lhs(F), rhs(F)

This time we can more naturally define the integrals from the Robin condition as r*(u - s)*v*ds(i).

The complete code can be found in the function solver_bcs in the program ft10_poisson_extended.py.

Test problem

We will use the same exact solution \( \uex=1+x^2+2y^2 \) as in the chapter Fundamentals: Solving the Poisson equation, and thus take \( \kappa=1 \) and \( f=-6 \). Our domain is the unit square, and we assign Dirichlet conditions at \( x=0 \) and \( x=1 \), a Robin condition at \( y=0 \), and a Neumann condition at \( y=1 \). With the given exact solution \( \uex \), we realize that the Neumann condition at \( y=1 \) is \( -\partial u / \partial n = - \partial u / \partial y = 4y = 4 \), while the Robin condition at \( y=0 \) can be selected in many ways. Since \( -\partial u/\partial n=\partial u/\partial y=0 \) at \( y=0 \), we can select \( s=\uex \) and specify \( r \neq 0 \) arbitrarily in the Robin condition. We will set \( r = 1000 \) and \( s = \uex \).

The boundary parts are thus \( \GD^0 \): \( x=0 \), \( \GD^1 \): \( x=1 \), \( \GR^0 \): \( y=0 \), and \( \GN^0 \): \( y=1 \).

When implementing this test problem, and especially other test problems with more complicated expressions, it is advantageous to use symbolic computing. Below we define the exact solution as a sympy expression and derive other functions from their mathematical definitions. Then we turn these expressions into C/C++ code, which can then be used to define Expression objects.

# Define manufactured solution in sympy and derive f, g, etc.
import sympy as sym
x, y = sym.symbols('x[0], x[1]')            # needed by UFL
u = 1 + x**2 + 2*y**2                       # exact solution
u_e = u                                     # exact solution
u_00 = u.subs(x, 0)                         # restrict to x = 0
u_01 = u.subs(x, 1)                         # restrict to x = 1
f = -sym.diff(u, x, 2) - sym.diff(u, y, 2)  # -Laplace(u)
f = sym.simplify(f)                         # simplify f
g = -sym.diff(u, y).subs(y, 1)              # compute g = -du/dn
r = 1000                                    # Robin data, arbitrary
s = u                                       # Robin data, u = s

# Collect variables
variables = [u_e, u_00, u_01, f, g, r, s]

# Turn into C/C++ code strings
variables = [sym.printing.ccode(var) for var in variables]

# Turn into FEniCS Expressions
variables = [Expression(var, degree=2) for var in variables]

# Extract variables
u_e, u_00, u_01, f, g, r, s = variables

# Define boundary conditions
boundary_conditions = {0: {'Dirichlet': u_00},   # x = 0
                       1: {'Dirichlet': u_01},   # x = 1
                       2: {'Robin':     (r, s)}, # y = 0
                       3: {'Neumann':   g}}      # y = 1

The complete code can be found in the function demo_bcs in the program ft10_poisson_extended.py.

Debugging boundary conditions

It is easy to make mistakes when implementing a problem with many different types of boundary conditions, as in the present case. One method to debug boundary conditions is to run through all vertex coordinates and check if the SubDomain.inside method marks the vertex as on the boundary. Another useful method is to list which degrees of freedom that are subject to Dirichlet conditions, and for first-order Lagrange (\( \mathsf{P}_1 \)) elements, print the corresponding vertex coordinates as illustrated by the following code snippet:

if debug1:

    # Print all vertices that belong to the boundary parts
    for x in mesh.coordinates():
        if bx0.inside(x, True): print('%s is on x = 0' % x)
        if bx1.inside(x, True): print('%s is on x = 1' % x)
        if by0.inside(x, True): print('%s is on y = 0' % x)
        if by1.inside(x, True): print('%s is on y = 1' % x)

    # Print the Dirichlet conditions
    print('Number of Dirichlet conditions:', len(bcs))
    if V.ufl_element().degree() == 1:  # P1 elements
        d2v = dof_to_vertex_map(V)
        coor = mesh.coordinates()
        for i, bc in enumerate(bcs):
            print('Dirichlet condition %d' % i)
            boundary_values = bc.get_boundary_values()
            for dof in boundary_values:
                print('   dof %2d: u = %g' % (dof, boundary_values[dof]))
                if V.ufl_element().degree() == 1:
                    print('    at point %s' %

Calls to the inside method. In the code snippet above, we call the inside method for each coordinate of the mesh. We could also place a printout inside the inside method. Then it will be surprising to see that this method is called not only for the points assoicated with degrees of freedom. For \( \mathsf{P}_1 \) elements the method is also called for each midpoint on each facet of the cells. This is because a Dirichlet condition is by default set only if the entire facet can be said to be subject to the condition defining the boundary.