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:
u is reconstructed from u_r;
local element contributions are evaluated using u;
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.