20.2. Boundary Conditions

We discuss more general boundary conditions for the Poisson equation. The whole boundary is split into three non-overlapping parts: \(\partial \Omega = \Gamma_D \cup \Gamma_N \cup \Gamma_R\)

• Dirichlet boundary conditions

  \[ u(x) = u_D(x) \qquad \forall \, x \in \Gamma_D \]

  Dirichlet boundary conditions are used, for example, to specify the temperature at the boundary.

• Neumann boundary conditions

  \[ \frac{\partial u}{\partial n}(x) = g(x) \qquad \forall \, x \in \Gamma_N \]

  Neumann boundary conditions model the given heat flux at the boundary, for example \(g = 0\) in case of thermal insulation.

• Robin boundary conditions

  \[ \alpha u(x) + \frac{\partial u}{\partial n}(x) = \alpha u_0(x) \qquad \forall \, x \in \Gamma_R \]

  Robin boundary condition model the heat transmission at the boundary, the flux \(\partial_n u \sim (u_0-u)\), where \(u_0-u\) is the temperature difference between the environment and the domain.

In the derivation of the weak form we had

\[ \int_\Omega \nabla u \nabla v - \int_{\Gamma_D \cup \Gamma_N \cup \Gamma_R} \frac{\partial u }{\partial n} v = \int_\Omega f v \qquad \forall \, v \in H^1(\Omega) \]

20.2.1. Natural boundary conditions

Neumann and Robin boundary conditions are the simple ones, they can be included naturally into the variational form.

Assume we don’t have a Dirichlet boundary, i.e. \(\partial \Omega = \Gamma_N \cup \Gamma_R\). We replace \(\frac{\partial u}{\partial n}\) by the boundary condition:

\[ \int_\Omega \nabla u \nabla v - \int_{\Gamma_N} g v - \int_{\Gamma_R} \alpha (u_0-u) v = \int f v \qquad \forall \, v \in H^1(\Omega) \]

We see that \(\int g v\) and \(\int \alpha u_0 v\) are linear terms in the test-function, and thus belong to the linear-form \(f\):

\[ \text{find } u \in H^1 : \quad \int_{\Omega} \nabla u \nabla v + \int_{\Gamma_R} \alpha u v = \int_\Omega f v + \int_{\Gamma_N} g v + \int_{\Gamma_R} \alpha u_0 v \]

Thus, we have to adapt the bilinear-form and linear-form:

\[ A(u,v) = \int_{\Omega} \nabla u \nabla v \, dx + \int_{\Gamma_R} \alpha u v \, ds \]

\[ f(v) = \int_\Omega f v \, dx + \int_{\Gamma_N} g v \, ds + \int_{\Gamma_R} \alpha u_0 v \, ds \]

from ngsolve import *
from ngsolve.webgui import Draw

from netgen.occ import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

fes = H1(mesh, order=3)
u,v = fes.TnT()

alpha = 10
u0 = 293
g = 10

a = BilinearForm(grad(u)*grad(v)*dx + alpha*u*v*ds("bottom"))
f = LinearForm (g * v*ds("right") + alpha*u0*v*ds("bottom"))

gfu = GridFunction(fes)

# assemble and solve
a.Assemble()
f.Assemble()
gfu.vec.data = a.mat.Inverse() * f.vec

Draw (gfu, mesh)
print ("du/dx:")
Draw (grad(gfu)[0], mesh);

du/dx:

20.2.2. Essential boundary conditions

Essential boundary conditions are less natural: We have to set the solution field to the given Dirichlet values, and restrict the test-functions to 0 on the Dirichlet boundary:

\[ \text{find } u \in H^1, u = u_D \text{ on } \Gamma_D \text{ s.t. } \quad A(u,v) = f(v) \quad \forall \, v \in H^1, v = 0 \text{ on } \Gamma_D \]

We split the solution vector into two parts: The given coefficients on the Dirichlet boundary, and all other including internal coefficients and coefficients on the natural boundaries:

\[\begin{split} u = \left( \begin{array}{c} u_D \\ u_f \end{array} \right) \end{split}\]

(f like free). Accordingly, the matrix and the right hand side are split as

\[\begin{split} A = \left( \begin{array}{cc} A_{DD} & A_{Df} \\ A_{fD} & A_{ff} \end{array} \right) \qquad f = \left( \begin{array}{c} f_D \\ f_f \end{array} \right) \end{split}\]

The test functions are reduced to the free nodes. Thus, the equations to solve are

\[ A_{fD} u_D + A_{ff} u_f = f_f \]

The given \(u_D\) is moved to the right hand side, and thus

\[ u_f = {A_{ff}}^{-1} (f_f - A_{fD} u_D). \]

It is convenient to assemble the complete matrix, invert a part of the matrix, end extend by zero:

\[\begin{split} A^{-1,ff} := \left( \begin{array}{cc} 0 & 0 \\ 0 & A_{ff}^{-1} \end{array} \right) \end{split}\]

The case of homogeneous Dirichlet boundary condtions \(u_D = 0\) is simple: We just say

\[ u = A^{-1,ff} f \]

In NGSolve, the finite element space maintains a BitArray marking the free degrees of freedom. It is used to invert the sub-matrix:

# fes = H1(mesh, dirichlet="left|bottom")
# u.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec

Non-homogeneous Dirichlet b.c. are reduced to homogeneous once as follows: Choose some function \(\widetilde u\) such that

\[ \widetilde u = u_D \quad \text{on } \Gamma_D \]

and set \(u = \widetilde u + w\). So we have to find a correction \(w\) with \(w = 0\) on \(\Gamma_D\) and

\[ A(w,v) = f(v) - A(\widetilde u , v) \qquad \forall \, v \text{ with } v = 0 \text{ on } \Gamma_D \]

In matrix notation, the correction \(w\) is

\[ w = A^{-1,ff} (f - A \tilde u) \]

from netgen.occ import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

fes = H1(mesh, order=3, dirichlet="bottom|left|right")
u,v = fes.TnT()

a = BilinearForm(grad(u)*grad(v)*dx)
f = LinearForm (fes)

a.Assemble()
f.Assemble();

Now we set the Dirichlet boundary values and don’t worry about the rest. The Set - function of a GridFunction does some kind of interpolation.

gfu = GridFunction(fes)

uD = sin(2*pi*x)
gfu.Set (uD, definedon=mesh.Boundaries("bottom"))

ea = { "euler_angles" : [-60,0,105] }
Draw(gfu, **ea, deformation=True, scale=0.3);

Now we compute the correction:

r = f.vec - a.mat * gfu.vec
gfu.vec.data += a.mat.Inverse(fes.FreeDofs()) * r

Draw(gfu, **ea, deformation=True, scale=0.3);

print (fes.FreeDofs())

0: 00000000000000000000001111111110000000001111111111
50: 11111111111111111111111111111111111111111111111111
100: 11111111111111111111111111111111111100000000110011
150: 11001100111100111100111100111100111100111100111100
200: 11111111001111001111001111001111001111001111001111
250: 00111111111111111111111111111111111111111111111111
300: 11111111111100111100111100111100111100111100111100
350: 11110011111111111111111111111111111111111111111111
400: 11111111111111111111111111111111111111111111111111
450: 11111111111111111111111111111111111111111111111111
500: 11111111111111111111111111111111111111111111111111
550: 11111111111111111111111111111111111111111111111111
600: 11111111111111111111111111111111111111111111111111
650: 11111111111111111111111111111111111111111111111111
700: 11111111111111111111111111111111111111111111111111
750: 11111111111111111111111111111111111111111111111111
800: 11111111111111111111111111111111111111111111111111
850: 11111111111111111111111111111111111111111111111111
900: 11111111111111111111111111111111111111111111111111
950: 11111111111111111111111111111111111111111111111111
1000: 11111111111111111111111111111111111111111111111111
1050: 1111111111111111111111111111111111111111111111