Boundary-condition elimination#

Once \(\mathbf{K}\) and \(\mathbf{f}\) are assembled (Global assembly), the prescribed-displacement DOFs supplied through 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 \(\mathbf{u}_f \in \mathbb{R}^{N_f}\) precede the constrained ones \(\mathbf{u}_c \in \mathbb{R}^{N_c}\), \(N = N_f + N_c\). The partitioned linear system reads

(1)#\[\begin{split}\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},\end{split}\]

where \(\mathbf{u}_c\) is fully prescribed (the user- supplied values from each Model.fix(node, dof, value=) call) and \(\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 (1) rearranges to

(2)#\[\mathbf{K}_{ff}\, \mathbf{u}_f \;=\; \mathbf{f}_f \,-\, \mathbf{K}_{fc}\, \mathbf{u}_c,\]

which is a square symmetric-positive-definite system in \(\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 (1):

(3)#\[\mathbf{r}_c \;=\; \mathbf{K}_{cf}\, \mathbf{u}_f \,+\, \mathbf{K}_{cc}\, \mathbf{u}_c \,-\, \mathbf{f}_c.\]

Two implementation notes follow from this:

  • \(\mathbf{K}_{cf}\) is the transpose of \(\mathbf{K}_{fc}\) (symmetry), so we don’t store it separately; the reaction product re-uses the assembled upper-triangle structure.

  • \(\mathbf{K}_{cc}\) blocks are typically tiny (a handful of constrained DOFs per node × the number of constrained nodes), so the dense \(\mathbf{K}_{cc}\, \mathbf{u}_c\) product is a free post-step.

Why partition rather than Lagrange multipliers?#

Two standard approaches handle \(\mathbf{u}_c = \bar{\mathbf{u}}_c\):

  1. Partition + reduce. Drop the constrained rows / columns and solve (2) directly. Pros: the reduced \(\mathbf{K}_{ff}\) is symmetric positive- definite — every direct backend’s optimal path applies. The full \((N + N_c) \times (N + N_c)\) linear system never materialises. Cons: rebuilding \(\mathbf{K}_{ff}\) on every step requires bookkeeping on the free / constrained partition.

  2. Lagrange multipliers. Augment with the constraint \(\mathbf{C}\, \mathbf{u} = \bar{\mathbf{u}}_c\) as an auxiliary block:

    \[\begin{split}\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}.\end{split}\]

    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 \(\mathbf{K}_{ff}\) re-uses the CSR pattern from Global 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 Model.solve() is:

  1. model.fix(...) calls accumulate the \((\mathrm{node}, \mathrm{dof}, \mathrm{value})\) constraints into an internal list.

  2. Just before the static / modal solve, _bc_reduce_and_release (see femorph_solver._bc_reduction) partitions the global \(\mathbf{K}\) / \(\mathbf{M}\) / \(\mathbf{f}\) into free / constrained blocks.

  3. The reduced \(\mathbf{K}_{ff}\) is handed to the active linear backend; the result is scattered back into the global vector.

  4. The reaction vector \(\mathbf{r}_c\) from (3) populates StaticResult.reaction.

Modal analysis follows the same pattern, with the constrained block dropped from both \(\mathbf{K}\) and \(\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).