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.