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

with a **nonzero** Dirichlet boundary condition

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 just interested in the *automatic utility,* then please skip to the last part of the tutorial.

## Spaces with Dirichlet conditions on part of the boundary¶

```
[1]:
```

```
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
mesh.GetBoundaries()
```

```
[1]:
```

```
('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

Then, we set the space as follows:

```
[2]:
```

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

Compare this space with the one without the `dirichlet`

flag:

```
[3]:
```

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

```
[3]:
```

```
(129, 129)
```

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.

```
[4]:
```

```
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: 11111111111111111111111111111
free dofs of fes:
0: 00001111000011110000111111111111111111101010110111
50: 11111111101101101111111111111110101110111111111111
100: 11111111111111111111111111111
```

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

and

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

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

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

```
[5]:
```

```
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\).

```
[6]:
```

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

```
[6]:
```

```
BaseWebGuiScene
```

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.

```
[7]:
```

```
u, v = fes.TnT()
a = BilinearForm(fes, symmetric=True)
a += grad(u)*grad(v)*dx
a.Assemble()
```

```
[7]:
```

```
<ngsolve.comp.BilinearForm at 0x7f77ceac3730>
```

If \(A=\) `a.mat`

is the matrix just assembled, then we want to solve for

or

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

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

```
[8]:
```

```
f = LinearForm(fes)
f += 1*v*dx
f.Assemble()
r = f.vec.CreateVector()
r.data = f.vec - a.mat * gfu.vec
```

The implementation of

by sparse solvers is achieved by the following:

```
[9]:
```

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

```
[9]:
```

```
BaseWebGuiScene
```

## 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.

```
[10]:
```

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

```
WARNING: maxsteps is deprecated, use maxiter instead!
CG iteration 1, residual = 1.3567996411869832
CG iteration 2, residual = 0.7483785609652776
CG iteration 3, residual = 0.6889833652395589
CG iteration 4, residual = 0.5654260806499144
CG iteration 5, residual = 0.47533706455917296
CG iteration 6, residual = 0.2727362779053307
CG iteration 7, residual = 0.17822059813787386
CG iteration 8, residual = 0.10777829194512963
CG iteration 9, residual = 0.06145108121302209
CG iteration 10, residual = 0.04209335389817511
CG iteration 11, residual = 0.02688622849286256
CG iteration 12, residual = 0.017003894315589064
CG iteration 13, residual = 0.0113162393824333
CG iteration 14, residual = 0.007706696746039635
CG iteration 15, residual = 0.003970044884338698
CG iteration 16, residual = 0.0018386562110491644
CG iteration 17, residual = 0.001007830954001793
CG iteration 18, residual = 0.0005324563317848768
CG iteration 19, residual = 0.00029545746894748284
CG iteration 20, residual = 0.00012801576145707164
CG iteration 21, residual = 6.086297488979755e-05
CG iteration 22, residual = 3.078769310428112e-05
CG iteration 23, residual = 1.8984767787489565e-05
CG iteration 24, residual = 1.0360038364735514e-05
CG iteration 25, residual = 5.364802003113384e-06
CG iteration 26, residual = 2.5598220926214547e-06
CG iteration 27, residual = 1.5328957516271298e-06
CG iteration 28, residual = 8.968897598370481e-07
CG iteration 29, residual = 4.181374210203346e-07
CG iteration 30, residual = 1.8694088830583122e-07
CG iteration 31, residual = 1.0408737597479772e-07
CG iteration 32, residual = 4.6394901186507276e-08
CG iteration 33, residual = 2.3245623086990255e-08
CG iteration 34, residual = 1.0824134397758478e-08
```

```
[10]:
```

```
BaseWebGuiScene
```

```
[ ]:
```

```
```