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.
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 $$
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} $$
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.
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:
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}}
specifies
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)
bcs.append(bc)
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']
integrals_N.append(g*v*ds(i))
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']
integrals_R_a.append(r*u*v*ds(i))
integrals_R_L.append(r*s*v*ds(i))
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.
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.
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' %
(str(tuple(coor[d2v[dof]].tolist()))))
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.