Consider the linear system
\[ A x = b , \]
arising from a Poisson-like operator without any boundary treatment ( \(A\) is the "Neumann" operator), with Dirichlet conditions prescribed on a subset of degrees of freedom \(d\):
\[ x_d = g_d . \]
Reorder unknowns as
\[ x = \begin{pmatrix} x_i \\ x_d \end{pmatrix}, \quad A = \begin{pmatrix} A_{ii} & A_{id} \\ A_{di} & A_{dd} \end{pmatrix}, \quad b = \begin{pmatrix} b_i \\ b_d \end{pmatrix}. \]
Define the lifting vector
\[ g = \begin{pmatrix} 0 \\ g_d \end{pmatrix}, \]
i.e. zero on interior degrees of freedom and equal to the prescribed Dirichlet values on boundary degrees of freedom. Also define
\[ \hat x = \begin{pmatrix} x_{i} \\ 0 \end{pmatrix}, \]
i.e., \(\hat x_d = 0\), and \(x = \hat x + g\). Substitution into the linear system gives the lifted equation
\[ \begin{aligned} Ax &= A(\hat x + g) = A \hat x + A g = b \\ A \hat x &= b - A g . \end{aligned} \]
To enforce \(\hat x_d=0\) strongly while preserving symmetry, eliminate all off-diagonal couplings to Dirichlet degrees of freedom:
\[ A_{\mathrm{elim}} = \begin{pmatrix} A_{ii} & 0 \\ 0 & A_{dd} \end{pmatrix}. \]
Define the modified right-hand side
\[ b_{\mathrm{elim}} = b - A g , \]
and then overwrite its Dirichlet components by
\[ b_{\mathrm{elim},d} = A_{dd} g_d, \]
s.t., overall
\[ b_{\mathrm{elim}} = \begin{pmatrix} b_{\mathrm{elim},i} \\ b_{\mathrm{elim},d} \end{pmatrix} = \begin{pmatrix} b_i - A_{id} g_d \\ A_{dd} g_d \end{pmatrix}. \]
The eliminated system is
\[ A_{\mathrm{elim}} x = b_{\mathrm{elim}}. \]
The interior equations read
\[ A_{ii} x_i = b_i - A_{id} g_d , \]
while the boundary equations enforce
\[ A_{dd} x_d = A_{dd} g_d \;\Rightarrow\; x_d = g_d . \]
Thus, the Dirichlet conditions are imposed strongly, and the interior solution satisfies the standard reduced system.
This procedure is equivalent to a symmetric Gaussian elimination of Dirichlet degrees of freedom. It preserves:
The diagonal entries corresponding to Dirichlet nodes are retained, which is advantageous for conditioning and iterative solvers.
Let \(g\) be the lifting vector. Then:
\[ y = A g , \]
\[ b_{\mathrm{elim}} = b - y , \]
\[ b_{\mathrm{elim},d} = \mathrm{diag}(A)\, g_d . \]
The operator \(A_{\mathrm{elim}}\) is applied by symmetrically zeroing all off-diagonal entries in rows and columns corresponding to Dirichlet degrees of freedom. In a matrix-free context, this symmetric zeroing can be implemented by applying it to each local element matrix. In fact, this is exactly what is happening in the operators implemented in this framework (more details below).
Solving \(A_{\mathrm{elim}} x = b_{\mathrm{elim}}\) yields the solution of the original Dirichlet problem.
To solve the problem
\[ A_{\mathrm{elim}} x = b_{\mathrm{elim}} \]
that has been introduced above, we need several ingredients:
Having the operators instantiated, it only remains to build the right-hand side vector \(b_{\mathrm{elim}}\) and solve the system. To build the vector, helper functions such as terra::fe::strong_algebraic_dirichlet_enforcement_poisson_like should be used that perform all the necessary steps. Concretely, the steps are:
It then remains to pass the system:
\[ A_{\mathrm{elim}} x = b_{\mathrm{elim}} \]
so a linear solver.
\( x \) is the solution of the original problem. No further boundary corrections are necessary.
The elimination procedure above is implemented in the framework by zeroing out the off-diagonal entries in the local element matrices. Let, for instance, \(W\) be a wedge element with 6 nodes. Of those 6 nodes, 3 are interior ( \(j = 1, 2, 3\)), and the remaining 3 are on the Dirichlet boundary ( \(j = 4, 5, 6\)). The local element matrix \(A^W\) is of size \(6\times 6\), and has the following block-structure (similar to the global matrix after rearranging the degrees of freedom above):
\[ A^W = \begin{pmatrix} A^W_{ii} & A^W_{id} \\ A^W_{di} & A^W_{dd} \end{pmatrix}. \]
The symmetric elimination process simply zeroes out the off-diagonal entries in rows and columns corresponding to the boundary nodes:
\[ A_{\mathrm{elim}}^W = \begin{pmatrix} A^W_{ii} & 0 \\ 0 & \mathrm{diag}(A^W_{dd}) \end{pmatrix}. \]
This approach can be realized on-the-fly in a matrix-free context, without the need to store the full local element matrices.
Often we need to enforce Dirichlet boundary conditions on the velocity field in the context of a Stokes-like problem. The approach is very similar to the one described above, but with a few subtleties:
Assuming an inf-sup stable discretization, we consider the discrete Stokes problem
\[ Kx = \begin{pmatrix} A & B \\ C & 0 \end{pmatrix} \begin{pmatrix} u \\ p \end{pmatrix} = \begin{pmatrix} f_u \\ f_p \end{pmatrix} =b. \]
and a lifting vector \(g\) (with non-zero entries on the Dirichlet boundary nodes in the velocity component). In this definition, we explicitly allow for non-symmetric versions of \(K\) due to the possible compressibility of the fluid.
The treatment of the RHS vector, i.e., the computation of \(b_\mathrm{elim}\), works just as in the Poisson-like case. Note that although only velocity components are prescribed, the elimination may also alter the values in the pressure component of the RHS.
We perform the steps just as above:
The other subtlety is that of course the elimination process must also be carried out on the gradient and divergence-like operators \(B\) and \(C\). Specifically, we need to zero out rows in \(B\) and columns in \(C\) corresponding to the boundary nodes. Just as in the Poisson-like case, this can be realized on-the-fly in a matrix-free context.
Plate boundary conditions are just Dirichlet boundary conditions for the velocity on the plate surface. Nothing special here from the algebraic point of view. Refer to the Stokes-like case above.
Free-slip boundary conditions model fluid flow in which the velocity components at a boundary are not constrained in the tangential direction, while the normal velocity is fixed to zero.
Concretely, at free-slip boundaries, we want to enforce the following condition:
\[ \begin{aligned} u_n &= 0, &\qquad \text{no normal flow} \\ (\sigma n)_t &= 0 &\qquad \text{no tangential stress} \end{aligned} \]
where \(u_n\) is the velocity component in the normal direction, and
\[ \sigma = 2 \eta \dot \varepsilon(u) - p I \]
with
\[ \dot \varepsilon (u) = \frac{1}{2}(\nabla u + (\nabla u)^T). \]
To implement the free-slip boundary condition, we only enforce the first condition on the velocity component in the normal direction, since the second one is naturally satisfied by the weak formulation.
Let's look at the simple case of the isoviscous Stokes system:
\[ \begin{aligned} - \Delta u + \nabla p &= 0, \\ \nabla \cdot u &= 0 \end{aligned} \]
Multiplying the momentum equation by a test function \( v \) and integrating by parts yields
\[ \int_{\Omega} \nabla u : \nabla v - \int_{\Omega} p \nabla \cdot v - \int_{\partial \Omega} (\partial_n u - pn) \cdot v = 0. \]
Since we impose \(u \cdot n = 0\) we also have vanishing normal components of the test functions, i.e., \(v \cdot n = 0\), and therefore on the free-slip boundary
\[ (\partial_n u - pn) \cdot v_t = (\partial_n u)_t \cdot v_t \]
because \(n \cdot v_t = 0\). Since we prescribe no tangential traction, the weak formulation drops that boundary term and thus enforces that the traction is really 0.
We implement the free-slip condition by only strongly imposing zero normal flow at the boundary. To do this, the implementation follows the idea of Engelman (1982): "The implementation of normal and/or tangential boundary conditions in finite element codes for incompressible fluid flow". The imposition is accomplished by pre- and post-multiplication of the entire linear system with rotation matrices that transform between the canonical and a local normal-tangential coordinate system. In that latter coordinate system, we can modify the normal velocity component directly. Let \(R\) be the rotation matrix that transforms from the canonical to the local normal-tangential coordinate system, and \(R^{-1} = R^T\) the inverse. Then, the modified Stokes system reads
\[ \begin{pmatrix} RAR^T & RB^T \\ BR^T & 0 \end{pmatrix} \begin{pmatrix} Ru \\ p \end{pmatrix} = \begin{pmatrix} Rf_u \\ f_p \end{pmatrix} \]
Note that this system is equivalent to the original Stokes system. However, we can now easily access the normal velocity components of the vectors and therefore eliminate them strongly, just as in the sections above. Particularly, we want to set them to zero.
For efficiency, the rotations are only executed at the boundary nodes. So, technically \(R\) is a block-diagonal matrix with blocks of size \(3\times 3\) that are either the local rotation matrices or identity matrices.
To enhance usability, the rotation of \( u \) is executed on-the-fly in a matrix-free context, such that before and after the solution step, the solution vector \( u \) is stored in canonical coordinates.
1.9.8 