This page was generated from unit-1.3-dirichlet/dirichlet.ipynb.

1.3 Dirichlet boundary conditions

This tutorial goes in depth into the mechanisms required to solve the Dirichlet problem

\[-\Delta u = f \quad \text{ in } \Omega,\]

with a nonzero Dirichlet boundary condition

\[u|_{\Gamma_D} = g \quad \text{ on a boundary part } \Gamma_D \subset \partial\Omega.\]

The same mechanisms are used in solving boundary value problems involving operators other than the Laplacian. You will see how to perform these tasks in NGSolve: * extend Dirichlet data from boundary parts, * convert boundary data into a volume source, * reduce inhomogeneous Dirichlet case to the homogeneous case, and * perform all these tasks automatically within a utility.

If you are interested only in an automatic utility that does all these steps under the hood, then please skip to the last part of the tutorial.

Spaces with Dirichlet conditions on part of the boundary

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import unit_square

mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
mesh.GetBoundaries()

('bottom', 'right', 'top', 'left')

The unit_square has its boundaries marked as left, right, top and bottom. Suppose we want non-homogeneous Dirichlet boundary conditions on

\[\Gamma_D = \Gamma_{left} \cup \Gamma_{right}.\]

Then, we set the space as follows:

fes = H1(mesh, order=2, dirichlet="left|right")

Compare this space with the one without the dirichlet flag:

fs2 = H1(mesh, order=2)
fes.ndof, fs2.ndof    # total number of unknowns

(125, 125)

Thus, the dirichlet flag did not change ndof. In NGSolve the unknowns are split into two groups: * dirichlet dofs (or constrained dofs), * free dofs.

The facility FreeDofs gives a BitArray such that FreeDofs[dof] is True if and only if dof is a free degree of freedom.

print("free dofs of fs2 without \"dirichlet\" flag:\n",
      fs2.FreeDofs())
print("free dofs of fes:\n", fes.FreeDofs())

free dofs of fs2 without "dirichlet" flag:
 0: 11111111111111111111111111111111111111111111111111
50: 11111111111111111111111111111111111111111111111111
100: 1111111111111111111111111
free dofs of fes:
 0: 00001111000011110000111111111111111111010011011111
50: 11111110110110111111111111111010111011111111111111
100: 1111111111111111111111111

• The space fs2 without dirichlet flag has only free dofs (no dirichlet dofs).

• The other space fes has a few dofs that are marked as not free. These are the dofs that are located on the boundary regions we marked as dirichlet.

Extension of boundary data

We use the standard technique of reducing a problem with essential non-homogeneous boundary conditions to one with homogeneous boundary condition using an extension. The solution \(u\) in \(H^1\) satisfies

\[u|_{\Gamma_D} = g\]

and

\[\int_\Omega \nabla u \cdot \nabla v_0 = \int_\Omega f v_0\]

for all \(v_0\) in \(\in H_{0,D}^1 = \{ v \in H^1: v|_{\Gamma_D} = 0\}\). Split the solution

\[u = u_0 + u_D\]

where \(u_D\) is an extension of \(g\) into \(\Omega\). Then we only need to find \(u_0\) in \(H^1_{0,D}\) satisfying the homogeneous Dirichlet problem

\[\int_\Omega \nabla u_0 \cdot \nabla v_0 = \int_\Omega f v_0 - \int_\Omega \nabla u_D \cdot \nabla v_0\]

for all \(v_0\) in \(H_{0,D}^1\). These are the issues to consider in this approach:

• How to define an extension \(u_D\) in the finite element space?

• How to form and solve the system for \(u_0\)?

Let us address the first in the following example.

Suppose we are given that

\[g = \sin(y) \qquad \text{ on } \; \Gamma_D.\]

g = sin(y)

We interpolate \(g\) on the boundary of the domain and extend it to zero on the elements not having an intersection with \(\Gamma_D\).

gfu = GridFunction(fes)
gfu.Set(g, BND)
Draw(gfu);

The keyword BND tells Set that g need only be interpolated on those parts of the boundary that are marked dirichlet.

Thus, gfu now contains the extension \(u_D\). Next, we turn to set up the system for \(u_0\).

Forms and assembly

In NGSolve, bilinear and linear forms are defined independently of the dirichlet flags. Matrices and vectors are set up with respect to all unknowns (free or constrained) so they may be restricted to any group of unknowns later.

u, v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble();

If \(A=\) a.mat is the matrix just assembled, then we want to solve for

\[A (u_0 + u_D) = f \quad \Rightarrow \quad A u_0 = f - A u_D\]

or

\[\begin{split}\left( \begin{array}{cc} A_{FF} & A_{FD} \\ A_{DF} & A_{DD} \end{array} \right) \left( \begin{array}{c} u_{0,F} \\ 0 \end{array} \right) = \left( \begin{array}{c} {f}_F \\ {f}_D \end{array} \right) - \left( \begin{array}{cc} A_{FF} & A_{FD} \\ A_{DF} & A_{DD} \end{array} \right) \left( \begin{array}{c} u_{D,F} \\ u_{D,D} \end{array} \right)\end{split}\]

where we have block partitioned using free dofs (\(F\)) and dirichlet dofs (\(D\)) as if they were numbered consecutively (which may not be the case in practice) for ease of presentation. The first row gives

\[A_{FF} u_{0,F} = f_F - [A u_D]_F.\]

Since we have already constructed \(u_D\), we need to perform these next steps:

• Set up the right hand side from \(f\) and \(u_D\).

• Solve a linear system which involves only \(A_{FF}\).

• Add solution: \(u = u_0 + u_D\).

Solve for the free dofs

We need to assemble the right hand side of \(A_{FF} u_{0,F} = f_F - [A u_D]_F\), namely

\[r = f - A u_D.\]

f = LinearForm(1*v*dx).Assemble()
r = f.vec - a.mat * gfu.vec

The implementation of

\[\begin{split}u = u_D + \left( \begin{array}{cc} A_{FF}^{-1} & 0 \\ 0 & 0 \end{array} \right) r\end{split}\]

by sparse solvers is achieved by the following:

gfu.vec.data += a.mat.Inverse(freedofs=fes.FreeDofs()) * r
Draw(gfu);

The automatic utility BVP

NGSolve also provides a BVP facility in the solvers submodule, within which the above steps are performed automatically. You provide \(A\), \(f\), a grid function gfu with your boundary condition \(g\), and a preconditioner. Then BVP solves the problem with non-homogeneous Dirichlet boundary condition and overwrites gfu with the solution.

gfu.Set(g, BND)
c = Preconditioner(a,"local")   #<- Jacobi preconditioner
#c = Preconditioner(a,"direct") #<- sparse direct solver
c.Update()
solvers.BVP(bf=a, lf=f, gf=gfu, pre=c)
Draw(gfu)

CG iteration 1, residual = 1.3219112283104724
CG iteration 2, residual = 0.7394997024575486
CG iteration 3, residual = 0.7004532102522516
CG iteration 4, residual = 0.5522264752073291
CG iteration 5, residual = 0.4534567012521492
CG iteration 6, residual = 0.2734129900941939
CG iteration 7, residual = 0.1911773899084286
CG iteration 8, residual = 0.10796683450156781
CG iteration 9, residual = 0.048718669465369416
CG iteration 10, residual = 0.03259935779207689
CG iteration 11, residual = 0.02069778320711081
CG iteration 12, residual = 0.014717573799955764
CG iteration 13, residual = 0.008800511509260013
CG iteration 14, residual = 0.005523832428735901
CG iteration 15, residual = 0.002872336024428405
CG iteration 16, residual = 0.001362590231106108
CG iteration 17, residual = 0.0006670722719336465
CG iteration 18, residual = 0.0003559270578789513
CG iteration 19, residual = 0.0001814449062072086
CG iteration 20, residual = 9.463997169105238e-05
CG iteration 21, residual = 5.447680492062797e-05
CG iteration 22, residual = 2.594180421781381e-05
CG iteration 23, residual = 1.6156303938285905e-05
CG iteration 24, residual = 9.86686340931647e-06
CG iteration 25, residual = 5.006687407986134e-06
CG iteration 26, residual = 2.4373113551839774e-06
CG iteration 27, residual = 1.067564188042133e-06
CG iteration 28, residual = 6.638270930774624e-07
CG iteration 29, residual = 3.3285126210921684e-07
CG iteration 30, residual = 1.4040502809762196e-07
CG iteration 31, residual = 7.915425165707552e-08
CG iteration 32, residual = 4.008252875972726e-08
CG iteration 33, residual = 2.2817766749220224e-08
CG iteration 34, residual = 1.1067817979380727e-08

BaseWebGuiScene