Implementation of Essential Boundary Conditions

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

Motivation

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.

Implementation

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}) \;.

Then

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.