Boundary Conditions

22. 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 domain and the environment.

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) \]

22.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)
gfu = GridFunction(fes)
a = BilinearForm(fes)
f = LinearForm (fes)

u = fes.TrialFunction()
v = fes.TestFunction()

alpha = 10
u0 = 293
g = 1

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

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

Draw (gfu, mesh, "temperature")
Draw (grad(gfu), mesh, "gradient");

22.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. } 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")
gfu = GridFunction(fes)
a = BilinearForm(fes)
f = LinearForm (fes)

u = fes.TrialFunction()
v = fes.TestFunction()

uD = sin(2*pi*x)
g = 1

a += grad(u)*grad(v)*dx

# assemble and solve
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.Set (uD, definedon=mesh.Boundaries("bottom"))
Draw(gfu);

Now we compute the correction:

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

Draw(gfu);
Contents