Update affine constraints docs

docs/src/topics/constraints.md

@@ -2,51 +2,104 @@
 DocTestSetup = :(using Ferrite)
 ```
 
-# Constraints
+# Affine Constraints  
 
-PDEs can in general be subjected to a number of constraints,
+In FEM, affine (or linear) constraints are commonly used to couple different degrees of freedom (DOFs). These constraints arise in cases such as [periodic boundary conditions](@ref "Periodic boundary conditions"), where the DOFs (or nodes) on one face should deform periodically with a corresponding node on the opposing face:  
 
 ```math
-g_I(\underline{a}) = 0, \quad I = 1 \text{ to } n_c
+a_1 = a_2,
 ```
 
-where $g$ are (non-linear) constraint equations, $\underline{a}$ is a vector of the
-degrees of freedom, and $n_c$ is the number of constraints. There are many ways to
-enforce these constraints, e.g. penalty methods and Lagrange multiplier methods.
+or in hanging nodes, where the value at the hanging node should be a weighted combination of the adjacent nodes:  
 
-## Affine constraints
+```math
+a_1 = 0.5a_3 + 0.5a_4.
+```
+
+In the above equations ``a_i`` represents the degree of freedom at node ``i``. Furthermore, affine constraints can also been viewed as a generalization of Dirichlets BCs, e.g. ``a_1 = 2``, however,Dirichlet BCs are mathematically and implementation-wise much easier to handle.
+
+To enforce affine constraints, the original linear system is modified. To explain this process, consider a problem with six DOFs where we have the following linear system of equations (obtained from a standard finite element procedure): 
+
+```math
+\boldsymbol{K} \boldsymbol{a} = \boldsymbol{f},
+```
+
+which is subjected to the following constraints:
+
+```math
+a_1 = 5a_2 + 3a_3 + 1 \\
+a_4 = 2a_3 + 6a_5
+```
+
+To incorporate these constraints into the linear system, we first collect them into a constraint matrix, ``\boldsymbol{C}`` and a vector containing the inhomogeneities, ``\boldsymbol{g}``:  
+
+```math
+\begin{align}
+\boldsymbol{a} = \boldsymbol{C} \boldsymbol{a}_f + \boldsymbol{g}
+\end{align}
+```
 
-Affine or linear constraints can be handled directly in Ferrite. Such constraints can typically
-be expressed as:
+where
 
 ```math
-a_1 =  5a_2 + 3a_3 + 1 \\
-a_4 =  2a_3 + 6a_5 \\
-\dots
+C =
+\begin{bmatrix}
+5 & 3 & 0 & 0 \\
+1 & 0 & 0 & 0 \\
+0 & 1 & 0 & 0 \\
+0 & 2 & 6 & 0 \\
+0 & 0 & 1 & 0 \\
+0 & 0 & 0 & 1
+\end{bmatrix}, \quad
+g =
+\begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}
 ```
 
-where $a_1$, $a_2$ etc. are system degrees of freedom. In Ferrite, we can account for such constraint using the `ConstraintHandler`:
+and where ``\boldsymbol{a}_f = [a_2, a_3, a_5, a_6]`` contains the "free" DOFs, and ``a_1`` and ``a_4`` are dependent on the others. Next, the above equation is inserted into the original system, and if we pre-multiply with ``\boldsymbol{C}^T``, we get a reduced system of equations:  
+
+```math
+\hat{\boldsymbol{K}} \boldsymbol{a}_f = \hat{\boldsymbol{f}}, \quad \hat{\boldsymbol{K}} = \boldsymbol{C}^T \boldsymbol{K} \boldsymbol{C}, \quad \hat{\boldsymbol{f}} = \boldsymbol{C}^T(\boldsymbol{f} - \boldsymbol{K}\boldsymbol{g})
+```
+
+The reduced system of equations can used to solve for `` a_f ``, which can then be used to calculate the dependent DOFs. Ferrite has functionaliy for setting up the ``\hat{\boldsymbol{K}}`` and ``\hat{\boldsymbol{f}}`` in an efficient way.
+
+!!! note "Limitation" 
+    Ferrite currently cannot handle (untangle) constraints where a DOF is both *master* and *slave* DOF. For example, if we have two affine constraints such as: 
+    ```math
+    a_1 = 2a_2 + 4 \\
+    a_2 = 3a_3 + 1
+    ``` 
+    Ferrite will not be able to resolve this situation because `` a_2 `` is both a master and a slave DOF in different constraints. 
+
+### Affine Constraints in Ferrite  
+
+To explain how affine constraints are handled in Ferrite, we will use the same example as above. The constraint equations can be constructed with `Ferrite.AffineConstraints` and added to the `ConstraintHandler`:  
 
 ```julia
 ch = ConstraintHandler(dh)
-lc1 = AffineConstraint(1, [2 => 5.0, 3 => 3.0], 1)
-lc2 = AffineConstraint(4, [3 => 2.0, 5 => 6.0], 0)
+
+lc1 = Ferrite.AffineConstraint(1, [2 => 5.0, 3 => 3.0], 1)
+lc2 = Ferrite.AffineConstraint(4, [3 => 2.0, 5 => 6.0], 0)
+
 add!(ch, lc1)
 add!(ch, lc2)
 ```
 
-Affine constraints will affect the sparsity pattern of the stiffness matrix, and as such, it is important to also include
-the `ConstraintHandler` as an argument when creating the sparsity pattern:
+Affine constraints will affect the sparsity pattern of the stiffness matrix, and as such, it is important to also include the `ConstraintHandler` as an argument when creating the sparsity pattern:  
 
 ```julia
 K = allocate_matrix(dh, ch)
 ```
 
-### Solving linear problems
-To solve the system ``\underline{\underline{K}}\underline{a}=\underline{f}``, account for affine constraints the same way as for
-`Dirichlet` boundary conditions; first call `apply!(K, f, ch)`. This will condense `K` and `f` inplace (i.e
-no new matrix will be created). Note however that we must also call `apply!` on the solution vector after
-solving the system to enforce the affine constraints:
+To solve the system, we could create ``\boldsymbol{C}`` and ``\boldsymbol{g}`` with the function 
+
+```julia
+C, g = Ferrite.create_constraint_matrix(ch)
+```
+
+and condense ``\boldsymbol{K}`` and ``\boldsymbol{f}`` as described above. However, this is in general inefficient, since the operation `C'*K*C` will allocate temporary matrices. Instead, Ferrite provides a much more efficient way of doing this. Infact, the affine constraints can be accounted for in the same way as `Dirichlet` boundary conditions; first call `apply!(K, f, ch)` which will condense `K` and `f` 
+inplace (i.e no new matrix will be created), and secondly call `apply!` on 
+the solution vector after solving the system to enforce the affine constraints:
 
 ```julia
 # ...