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

.. math::

   x = T x_r + x_f \;.

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

.. math::

   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

.. math::

   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

.. math::

   f(u) = 0 \;,

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

.. math::

   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

.. math::

   K \Delta u = r \;,

where the tangent matrix :math:`K` and the residual :math:`r` are

.. math::

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

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

Then

.. math::

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

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

Linear Problems
^^^^^^^^^^^^^^^

For linear problems we have

.. math::

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

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

and so the Newton method converges in a single iteration:

.. math::

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

Evaluation of Residual and Tangent Matrix
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

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

- The EBCs are applied to the full DOF vector :math:`u`.
- The reduced vector :math:`u_r` is passed to the Newton solver.
- Newton iteration loop:

  - Evaluation of :math:`f_r` or :math:`K_r`:

    #. :math:`u` is reconstructed from :math:`u_r`;
    #. local element contributions are evaluated using :math:`u`;
    #. local element contributions are assembled into :math:`f_r` or
       :math:`K_r` - values corresponding to fixed DOF positions are thrown
       away.

  - The reduced system :math:`K_r \Delta u_r = r_r` is solved.
  - Solution is updated: :math:`u_r \leftarrow u_r - \Delta u_r`.
  - The loop is terminated if a stopping condition is satisfied, the solver
    returns the final :math:`u_r`.

- The final :math:`u` is reconstructed from :math:`u_r`.