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

2.1. Natural boundary conditions

In the derivation of the weak form we had

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

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\) on the right hand side:

\[ \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

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"))

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

The temperature:

Draw (gfu, mesh); 

The heat flux:

Draw (-grad(gfu), mesh, vectors=True)
Draw (grad(gfu)[0], mesh);

Exercise:

• Look at the components of the gradient (grad(gfu)[0] gives the x-component). Where do you find the given values of \(g\) for the Neumann bc?

• How does the solution change by changing the heat transmission coefficient \(\alpha\) ?

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(\Omega), u = u_D \text{ on } \Gamma_D \; \text{such that}\)

\[ 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, and padding by 0:

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

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

fes = H1(mesh, order=3, dirichlet="bottom|left|right")
u,v = fes.TnT()
gfu = GridFunction(fes)
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
f = LinearForm (fes)  # nothing here

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.

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

Now we compute the correction:

inv = a.mat.Inverse(fes.FreeDofs())

residuum = f.vec - a.mat * gfu.vec
gfu.vec.data += inv * residuum

Draw (gfu);

Exercise:

• Change the boundary conditions to homogeneous Neumann bc at the left and right boundaries

• Replace sin(2*pi*x) by higher frequencies

• To set different functions at different parts, use mesh.BoundaryCF( { "bottom" : sin(x), "top" : 1 }, default=0)

Draw (grad(gfu)[0], mesh)

BaseWebGuiScene