Table Of Contents

Previous topic

Notes on solving PDEs by the Finite Element Method

Next topic

Developer Guide

This Page

Implementation of Essential Boundary Conditions

The essential boundary conditions can be applied in several ways. Here we describe the implementation used in SfePy.


Let us solve a linear system A x = b with n \times n matrix A with n_f values in the x vector known. The known values can be for example EBC values on a boundary, if A comes from a PDE discretization. If we put the known fixed values into a vector x_f, that has the same size as x, and has zeros in positions that are not fixed, we can easily construct a n \times n_r matrix T that maps the reduced vector x_r of size n_r = n - n_f, where the fixed values are removed, to the full vector x:

x = T x_r + x_f \;.

With that the reduced linear system with a n_r \times n_r can be formed:

T^T A T x_r = T^T (b - A x_f)

that can be solved by a linear solver. We can see, that the (non-zero) known values are now on the right-hand side of the linear system. When the known values are all zero, we have simply

T^T A T x_r = T^T b \;,

which is convenient, as it allows simply throwing away the A and b entries corresponding to the known values already during the finite element assembling.


All PDEs in SfePy are solved in a uniform way as a system of non-linear equations

f(u) = 0 \;,

where f is the nonlinear function and u the vector of unknown DOFs. This system is solved iteratively by the Newton method

u^{new} = u^{old} - (\tdiff{f}{u^{old}})^{-1} f(u^{old})

until a convergence criterion is met. Each iteration involves solution of the system of linear equations

K \Delta u = r \;,

where the tangent matrix K and the residual r are

K \equiv \tdiff{f}{u^{old}} \;,

r \equiv f(u^{old}) \;.


u^{new} = u^{old} - \Delta u \;.

If the initial (old) vector u^{old} contains the values of EBCs at correct positions, the increment \Delta u is zero at those positions. This allows us to assemble directly the reduced matrix T^T K
T, the right-hand side T^T r, and ignore the values of EBCs during assembling. The EBCs are satisfied automatically by applying them to the initial guess u^{0}, that is given to the Newton solver.

Linear Problems

For linear problems we have

f(u) \equiv A u - b = 0 \;,

\tdiff{f}{u} = A \;,

and so the Newton method converges in a single iteration:

u^{new} = u^{old} - A^{-1} (A u^{old} - b) = A^{-1} b \;.

Evaluation of Residual and Tangent Matrix

The evaluation of the residual f as well as the tangent matrix K within the Newton solver proceeds in the following steps:

  • The EBCs are applied to the full DOF vector u.
  • The reduced vector u_r is passed to the Newton solver.
  • Newton iteration loop:
    • Evaluation of f_r or K_r:
      1. u is reconstructed from u_r;
      2. local element contributions are evaluated using u;
      3. local element contributions are assembled into f_r or K_r - values corresponding to fixed DOF positions are thrown away.
    • The reduced system K_r \Delta u_r = r_r is solved.
    • Solution is updated: u_r \leftarrow u_r - \Delta u_r.
    • The loop is terminated if a stopping condition is satisfied, the solver returns the final u_r.
  • The final u is reconstructed from u_r.