Douglas N. Arnold 2014

Specifying boundary conditions¶

Inhomogeneous natural boundary conditions¶

The first problem we considered was the PDE

$$-\Delta u + u = f$$ together with the homogeneous Neumann boundary condition $\partial u/\partial n=0$ on $\partial\Omega$. This is the easiest boundary condition to implement with finite elements: you have to do precisely nothing! (By contrast, Neumann boundary conditions are a bit of a chore for finite differences.) FEniCS can handle many other types of boundary conditions as well, just about all the boundary conditions that make sense for such an equation.

Only slightly more complicated than homogeneous Neumann boundary conditions are the inhomogeneous Neumann boundary conditions $\partial u/\partial n=g$ on $\partial\Omega$. This adds an additional term to the right hand side of the weak formulation, so $u\in V$ satisfies

$$b(u, v) = F(v), \quad v\in V,$$ where again $V=H^1(\Omega)$ and $$b(u,v)=\int_\Omega (\operatorname{grad} u \cdot \operatorname{grad}v+uv)\,dx,$$ but now $$F(v) = \int_\Omega fv\,dx + \int_{\partial\Omega} gv\,ds.$$ (Make sure that you understand this!) Of course the finite element solution $u_h\in V_h$ is still defined by $b(u_h,v)= F(v)$, $v\in V_h$, with this new definition of $F(v)$.

To implement the new term in FEniCS is straightforward. First we have to define an Expression for $g$, just as we did for $f$. For example, if $g(x, y) = x^2 y$, we would use

g = Expression('pow(x[0], 2)*x[1]')

Next we add the extra term to the definition of $F$, which becomes

F = f*v*dx + g*v*ds

The remainder of the program need not be changed.

Remark: Note that, in addition to the "measure" dx which is used to signal integration over the domain, here we use the measure ds which means integration over the boundary of the domain. It is also possible to single out a part of the domain or a part of the boundary and integrate over that with a measure like dx(1) or ds(3). We do not need this yet, but will for more complicated boundary conditions (see below).

Essential boundary conditions¶

Instead of the Neumann boundary condition, we can impose the Dirichlet boundary condition $u=g$. This is an essential boundary condition, meaning that it is imposed as a constraint on the space $V$ rather than by changing the bilinear form $b$ or the linear functional $F$. The weak formulation is now to find $u\in V=H^1(\Omega)$ such that $u=g$ on $\partial\Omega$ and satisfying

$$b(u,v) = \int_\Omega fv\,dx \quad\text{for all $v\in V$ such that $v=0$ on $\partial\Omega$}.$$ On the discrete level, $u_h$ is restricted to elements of $V_h$ which interpolate $g$ on the boundary, i.e., any DOFs associated to vertices or edges on the boundary must give the same value when applied to $u_h$ as when applied to $g$. We say that the DOFs of $u_h$ on the boundary are constrained to $g$. In the case of Lagrange elements of degree 2, for instance, this means that $u_h$ must equal $g$ at the vertices on the boundary and the midpoints of edges on the boundary. Similarly the DOFs of the test function $v$ on the boundary are constrained to be zero.

FEniCS imposes these constraints this after assembling the stiffness matrix and load vector, by modifying entries associated to boundary DOFs. For this, we need to specify to FEniCS which DOFs are to be constrained (in this case all the DOFs on the boundary), and what they are to be constrained to (namely $g$). This is accomplished by invoking the procedure

DirichletBC(Vh, g, bdry)

where the three arguments are the function space to constrain, the function to constrain to, and the portion of the domain where DOFs are to be constrained. The function $g$ can be defined as an Expression just as for the Neumann boundary condition. When the Dirichlet boundary condition is to be imposed on the whole boundary of $\Omega$, we can simply take bdry = DomainBoundary(). The relevant code is then:

g = Expression('...')
bc = DirichletBC(Vh, g, DomainBoundary())
uh = Function(Vh)
solve(b == F, uh, bc)

Note that the solve call now takes a third argument, specifying the essential boundary condition.

Mixed boundary conditions¶

Sometimes we have different boundary conditions on different, complementary, parts of the boundary. For example, let $\Gamma_D$ be the bottom side of the unit square, and $\Gamma_N$ the union of the remaining three sides, and consider the solution of the boundary value problem

$$-\Delta u = 0 \quad \text{in $\Omega$}, \quad u = g \quad \text{ on $\Gamma_D$}, \quad \frac{\partial u}{\partial n} =0 \quad \text{ on $\Gamma_N$}.$$ The weak formulation then seeks $u\in V=H^1(\Omega)$ constrained to equal $g$ on $\Gamma_D$ such that $$\int_\Omega\operatorname{grad} u\cdot\operatorname{grad}v\,dx = 0$$ for all $v\in V$ constrained to zero on $\Gamma_D$. In this case we cannot use DomainBoundary(), but must define the subset $\Gamma_D$. This is done by defining a python function and passing it as the third argument of DirichletBC. We might choose to call the function GammaD. It takes two arguments, point x, and a Boolean on_boundary. When FEniCS needs to identify whether some point of the domain belongs to $\Gamma_D$ it will call GammaD(x, on_boundary) with x set to the point (the array of its coordinates x[0] and x[1]), and on_boundary set to True or False according to whether the point $x$ belongs to $\partial\Omega$ or not. Our function must return the value True or False according to whether the point $x$ belongs to $\Gamma_D$ or not. Since $\Gamma_D$ is the bottom edge of the square, the test is simply whether the second coordinate is equal to zero. However one should never test floating point numbers for exact equality, since they may suffer from small round off errors. FEniCS provides near(x[1], 0.) for this. It returns True if x[1] is within round off error of zero. Thus we can specify $\Gamma_D$ as follows:

def GammaD(x, on_boundary):
    return near(x[1], 0)

From this point on the code is similar to that discussed above.

g = Expression('...')
bc = DirichletBC(Vh, g, GammaD)
uh = Function(Vh)
solve(b == F, uh, bc)

Note that for this problem, with $f=0$ and no inhomogeneous Neumann data, the linear functional $F=0$. This can be entered as F = Constant(0.)*v*dx.

Remark. People are often confused about what to do at the interface where the boundary conditions switch from Neumann to Dirichlet. In our example, this would consist of the two points $(0,0)$ and $(1,0)$, the end points of the bottom side $\Gamma_D$. Should the DOFs associated with these points be constrained, as for a Dirichlet BC or left unconstrained as for a Neumann BC? The answer is that they should be constrained. The constraint $u=g$ is a continuous condition: if it holds on a subset of $\partial\Omega$, it should hold on the closure of that set.

Multiple Dirichlet boundary conditions¶

It is possible to have multiple Dirichlet boundary conditions, say $u=g_1$ on the bottom side and $u=g_2$ on the top side. We would then define two different DirichletBC, each just as above, say bc1 and bc2, and then call solve as

solve(b == F, uh, [bc1, bc2])

Neumann BC on part of the boundary¶

When inhomogeneous Neumann conditions are imposed on part of the boundary, we may need to include an integral like $\int_{\Gamma_N} gv\,ds$ in the linear functional $F$. If we can define the expression $g$ on the whole boundary, but so that it is zero except on $\Gamma_N$ (extension by zero), we can simply write this integral as $\int_{\partial\Omega} gv\,ds$ and nothing new is needed. Otherwise, we need to define a "measure" ds(1) which integrates only over $\Gamma_N$. The same approach can be used for multiple Neumann boundary conditions. The definition of such special measures is confusing because it involves several complex datastructures. One must first define a MeshFunction which assigns an integer to each edge, then define the relevant portion of the boundary using a subclass of the FEniCS class SubDomain, then use this subclass to set the values of the mesh function, and then use the mesh function to make the new measure. Here is some typical code for this, which redefines ds so that ds(1) integrates over the left half of the boundary (where the first coordinate is $<1/2$), ds(0) integrates over the right half of the boundary, and ds() integrates over the whole boundary.

# create a mesh function which assigns an unsigned integer (size_t) to each edge
mf = MeshFunction("size_t", mesh, 1) # 3rd argument is dimension of an edge
mf.set_all(0) # initialize the function to zero
# create a SubDomain subclass, specifying the portion of the boundary with x[0] < 1/2
class Lefthalf(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < 0.5 and on_boundary
lefthalf = Lefthalf() # instantiate it
# use this lefthalf object to set values of the mesh function to 1 in the subdomain
lefthalf.mark(mf, 1)
# define a new measure ds based on this mesh function
ds = Measure("ds")[mf]

Robin boundary conditions¶

Robin boundary conditions take the form $\partial u/\partial n + \alpha u = g$, i.e., a linear combination of Dirichlet and Neumann boundary conditions. These can be treated much as for Neumann boundary conditions, in that they are natural, not essential. In the case of Robin boundary conditions, not only is the functional $F$ modified, but also the bilinear form $b$.

Periodic boundary conditions¶

In the case of periodic boundary conditions, the boundary condition does not constrain individual points on the boundary, but relates the values at two different points. We will not discuss this further except to say that this possibility is available in FEniCS.

# this cell is present to set the notebook style
from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/dna1.css", "r").read()
    return HTML(styles)
css_styling()