Boundary-condition elimination ============================== Once :math:`\mathbf{K}` and :math:`\mathbf{f}` are assembled (:doc:`assembly`), the prescribed-displacement DOFs supplied through :meth:`Model.fix() ` are **eliminated by partitioning** before the linear solver sees the system. This chapter covers the partition algebra, why femorph-solver doesn't use Lagrange multipliers for simple Dirichlet constraints, and how the constrained-DOF reactions fall out for free. Partition into free / constrained DOFs -------------------------------------- Order the global DOF vector so the free (unconstrained) DOFs :math:`\mathbf{u}_f \in \mathbb{R}^{N_f}` precede the constrained ones :math:`\mathbf{u}_c \in \mathbb{R}^{N_c}`, :math:`N = N_f + N_c`. The partitioned linear system reads .. math:: :label: partitioned \begin{bmatrix} \mathbf{K}_{ff} & \mathbf{K}_{fc} \\ \mathbf{K}_{cf} & \mathbf{K}_{cc} \end{bmatrix} \begin{bmatrix} \mathbf{u}_f \\ \mathbf{u}_c \end{bmatrix} \;=\; \begin{bmatrix} \mathbf{f}_f \\ \mathbf{f}_c \end{bmatrix}, where :math:`\mathbf{u}_c` is **fully prescribed** (the user- supplied values from each ``Model.fix(node, dof, value=)`` call) and :math:`\mathbf{f}_c` collects the applied loads on the constrained DOFs (typically zero for pure-Dirichlet rows; non-zero only when the user has applied a force on a node that's also pinned in another DOF). The first block of :eq:`partitioned` rearranges to .. math:: :label: reduced \mathbf{K}_{ff}\, \mathbf{u}_f \;=\; \mathbf{f}_f \,-\, \mathbf{K}_{fc}\, \mathbf{u}_c, which is a square symmetric-positive-definite system in :math:`\mathbf{u}_f` only. This is what the linear backend (Pardiso / CHOLMOD / MUMPS / SuperLU) actually factors; the constrained rows / columns are dropped before the factor. Reactions follow from the second block of :eq:`partitioned`: .. math:: :label: reactions \mathbf{r}_c \;=\; \mathbf{K}_{cf}\, \mathbf{u}_f \,+\, \mathbf{K}_{cc}\, \mathbf{u}_c \,-\, \mathbf{f}_c. Two implementation notes follow from this: * :math:`\mathbf{K}_{cf}` is the transpose of :math:`\mathbf{K}_{fc}` (symmetry), so we don't store it separately; the reaction product re-uses the assembled upper-triangle structure. * :math:`\mathbf{K}_{cc}` blocks are typically tiny (a handful of constrained DOFs per node × the number of constrained nodes), so the dense :math:`\mathbf{K}_{cc}\, \mathbf{u}_c` product is a free post-step. Why partition rather than Lagrange multipliers? ----------------------------------------------- Two standard approaches handle :math:`\mathbf{u}_c = \bar{\mathbf{u}}_c`: 1. **Partition + reduce.** Drop the constrained rows / columns and solve :eq:`reduced` directly. Pros: the reduced :math:`\mathbf{K}_{ff}` is symmetric positive- definite — every direct backend's optimal path applies. The full :math:`(N + N_c) \times (N + N_c)` linear system never materialises. Cons: rebuilding :math:`\mathbf{K}_{ff}` on every step requires bookkeeping on the free / constrained partition. 2. **Lagrange multipliers.** Augment with the constraint :math:`\mathbf{C}\, \mathbf{u} = \bar{\mathbf{u}}_c` as an auxiliary block: .. math:: \begin{bmatrix} \mathbf{K} & \mathbf{C}^{\!\top} \\ \mathbf{C} & 0 \end{bmatrix} \begin{bmatrix} \mathbf{u} \\ \boldsymbol{\lambda} \end{bmatrix} \;=\; \begin{bmatrix} \mathbf{f} \\ \bar{\mathbf{u}}_c \end{bmatrix}. Pros: the user-facing matrix never changes pattern when constraints are added / removed. Cons: the augmented system is **symmetric indefinite** (the zero block on the diagonal kills positive-definiteness), so direct backends need a Bunch–Kaufman LDL factorisation rather than the faster Cholesky. Iterative solvers like CG can't be used directly without pre-conditioning the saddle-point structure. For **simple Dirichlet** BCs (the only kind femorph-solver exposes today through ``Model.fix``), partition + reduce wins because: * Linear-elastic problems are symmetric positive-definite by construction; throwing that away to use Lagrange would cost a 2× factorisation slow-down for no gain. * The free / constrained partition is cheap to compute once the user finishes calling ``Model.fix`` — DOFs are permuted in-place and the partitioning of :math:`\mathbf{K}_{ff}` re-uses the CSR pattern from :doc:`assembly`. Lagrange-multiplier paths *are* needed for **multi-point constraints** (MPC: rigid-body links, periodic-cyclic boundary conditions across non-conforming faces, contact gap closure) where the constraint cannot be reduced to a single-DOF prescription. femorph-solver's cyclic-symmetry boundary condition is a periodic MPC and uses the **augmented complex- Hermitian harmonic decomposition** (Grimes, Lewis, Simon 1994), which is mathematically equivalent to a Lagrange multiplier on a per-harmonic basis. How femorph-solver runs the elimination --------------------------------------- The flow inside :meth:`Model.solve() ` is: 1. ``model.fix(...)`` calls accumulate the :math:`(\mathrm{node}, \mathrm{dof}, \mathrm{value})` constraints into an internal list. 2. Just before the static / modal solve, ``_bc_reduce_and_release`` (see :mod:`femorph_solver._bc_reduction`) partitions the global :math:`\mathbf{K}` / :math:`\mathbf{M}` / :math:`\mathbf{f}` into free / constrained blocks. 3. The reduced :math:`\mathbf{K}_{ff}` is handed to the active linear backend; the result is scattered back into the global vector. 4. The reaction vector :math:`\mathbf{r}_c` from :eq:`reactions` populates :attr:`StaticResult.reaction`. Modal analysis follows the same pattern, with the constrained block dropped from both :math:`\mathbf{K}` and :math:`\mathbf{M}` before the eigenvalue solve. The free-DOF eigenvectors are scattered back to the global node ordering on the way out. References ---------- * Saad, Y. (2003) *Iterative Methods for Sparse Linear Systems*, 2nd ed., SIAM, §3.3 (sparse partitioning) + §13 (saddle-point linear systems / Lagrange multipliers). * Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) *Concepts and Applications of Finite Element Analysis*, 4th ed., Wiley, §2.10 (Dirichlet partitioning), §3.6 (penalty vs Lagrange-multiplier MPC). * Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed., §3.4 (boundary-condition imposition), §10.3 (cyclic- symmetry MPC and the harmonic decomposition). * Strang, G. and Fix, G. (2008) *An Analysis of the Finite Element Method*, 2nd ed., Wellesley-Cambridge Press, §1.5 (essential vs natural BCs in the variational form). * Bunch, J. R. and Kaufman, L. (1977) "Some stable methods for calculating inertia and solving symmetric linear systems," *Mathematics of Computation* 31, 162–179 (the LDL factorisation used when an indefinite augmented system is unavoidable).