7. Dirichlet boundary conditions

The Helmholtz problem we solved in the previous part was chosen to have homogeneous Neumann or natural boundary conditions, which can be implemented simply by cancelling the zero surface integral. We can now instead consider the case of Dirichlet, or essential boundary conditions. Instead of the Helmholtz problem we solved before, let us now specify a Poisson problem with homogeneous Dirichlet conditions, find u in some finite element space V such that:

(7.1)2u=fu=0 on Γ

In order to implement the Dirichlet conditions, we need to decompose V into two parts:

(7.2)V=V0VΓ

where VΓ is the space spanned by those functions in the basis of V which are non-zero on Γ, and V0 is the space spanned by the remaining basis functions (i.e. those basis functions which vanish on Γ). It is a direct consequence of the nodal nature of the basis that the basis functions for VΓ are those corresponding to the nodes on Γ while the basis for V0 is composed of all the other functions.

We now write the weak form of (7.1), find u=u0+uΓ with u0V0 and uΓVΓ such that:

(7.3)Ωv0(u0+uΓ)dxΓv0(u0+uΓ)nds=0=Ωv0fdxv0V0uΓ=0 on Γ

There are a number of features of this equation which require some explanation:

  1. We only test with functions from V0. This is because it is only necessary that the differential equation is satisfied on the interior of the domain: on the boundary of the domain we need only satisfy the boundary conditions.

  2. The surface integral now cancels because v0 is guaranteed to be zero everywhere on the boundary.

  3. The uΓ definition actually implies that uΓ=0 everywhere, since all of the nodes in VΓ lie on the boundary.

This means that the weak form is actually:

(7.4)Ωv0udx=Ωv0fdxv0V0uΓ=0

7.1. An algorithm for homogeneous Dirichlet conditions

The implementation of homogeneous Dirichlet conditions is actually rather straightforward.

  1. The system is assembled completely ignoring the Dirichlet conditions. This results in a global matrix and vector which are correct on the rows corresponding to test functions in V0, but incorrect on the VΓ rows.

  2. The global vector rows corresponding to boundary nodes are set to 0.

  3. The global matrix rows corresponding to boundary nodes are set to 0.

  4. The diagonal entry on each matrix row corresponding to a boundary node is set to 1.

This has the effect of replacing the incorrect boundary rows of the system with the equation ui=0 for all boundary node numbers i.

Hint

This algorithm has the unfortunate side effect of making the global matrix non-symmetric. If a symmetric matrix is required (for example in order to use a symmetric solver), then forward substition can be used to zero the boundary columns in the matrix, but that is beyond the scope of this module.

7.2. Implementing boundary conditions

Let:

f=(16π2(x11)2x212(x11)28(x11)x12x21)sin(4πx0)

With this definition, (7.4) has solution:

u=sin(4πx0)(x11)2x21
Exercise

fe_utils/solvers/poisson.py contains a partial implementation of this problem. You need to implement the assemble() function. You should base your implementation on your fe_utils/solvers/helmholtz.py but take into account the difference in the equation, and the boundary conditions. The fe_utils.solvers.poisson.boundary_nodes() function in fe_utils/solvers/poisson.py is likely to be helpful in implementing the boundary conditions. As before, run:

python fe_utils/solvers/poisson.py --help

for instructions (they are the same as for fe_utils/solvers/helmholtz.py). Similarly, test/test_12_poisson_convergence.py contains convergence tests for this problem.

7.3. Inhomogeneous Dirichlet conditions

The algorithm described here can be extended to inhomogeneous systems by setting the entries in the global vector to the value of the boundary condition at the corresponding boundary node. This additional step is required for the mastery exercise, but will be explained in more detail in the next section.