How to efficiently implement Dirichlet boundary conditions in global sparse finite element stiffnes matrices

Question

I am wondering how Dirichlet boundary conditions in global sparse finite element matrices are actually implemented efficiently. For example lets say that our global finite element matrix was: $$K = \begin{bmatrix} 5 & 2 & 0 & -1 & 0 \\ 2 & 4 & 1 & 0 & 0 \\ 0 & 1 & 6 & 3 & 2 \\ -1 & 0 & 3 & 7 & 0 \\ 0 & 0 & 2 & 0 & 3 \end{bmatrix} \hspace{5mm}\text{and right-hand side vector}\hspace{5mm} b = \begin{bmatrix} b1 \\ b2 \\ b3 \\ b4 \\ b5 \\ \end{bmatrix} $$

Then to apply a Dirichlet condition on the first node ($x_{1}=c$) we would zero out the first row, put a 1 at $K_{11}$, and subtract the first column from the right-hand side. For example our system would become: $$K = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 4 & 1 & 0 & 0 \\ 0 & 1 & 6 & 3 & 2 \\ 0 & 0 & 3 & 7 & 0 \\ 0 & 0 & 2 & 0 & 3 \end{bmatrix} \hspace{5mm}\text{and right-hand side vector}\hspace{5mm} b = \begin{bmatrix} c \\ b2-2\times{c} \\ b3-0\times{c} \\ b4+1\times{c} \\ b5-0\times{c} \\ \end{bmatrix} $$

This is all well and good in theory, but if our K matrix is stored in compressed row format (CRS) then moving the columns to the right-hand side becomes expensive for large systems (with many nodes being dirichlet). An alternative would be to not move the columns corresponding to a Dirichlet condition to the right-hand side, i.e. our system would become: $$K = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 2 & 4 & 1 & 0 & 0 \\ 0 & 1 & 6 & 3 & 2 \\ -1 & 0 & 3 & 7 & 0 \\ 0 & 0 & 2 & 0 & 3 \end{bmatrix} \hspace{5mm}\text{and right-hand side vector}\hspace{5mm} b = \begin{bmatrix} c \\ b2 \\ b3 \\ b4 \\ b5 \\ \end{bmatrix} $$

This however has a major draw back in that the system is no longer symmetric and so we could no longer use preconditioned conjugate gradient (or other symmetric solvers). One interesting solution that I came across is the "Method of Large Numbers" which I found in the book "Programming Finite Elements in Java" by Gennadiy Nikishkov. This method uses the fact that double precision only contains around 16 digits of accuracy. Instead of putting a 1 in the $K_{11}$ position we place a large number. For example our system becomes: $$K = \begin{bmatrix} 1.0e64 & 2 & 0 & -1 & 0 \\ 2 & 4 & 1 & 0 & 0 \\ 0 & 1 & 6 & 3 & 2 \\ -1 & 0 & 3 & 7 & 0 \\ 0 & 0 & 2 & 0 & 3 \end{bmatrix} \hspace{5mm}\text{and right-hand side vector}\hspace{5mm} b = \begin{bmatrix} c\times{1.0e64} \\ b2 \\ b3 \\ b4 \\ b5 \\ \end{bmatrix} $$

The advantages of this method are that it maintains the symmetry of the matrix while also being very efficient for sparse storage formats. My questions then are as follows:

How are Dirichlet boundary conditions typically implemented in finite element codes for heat/fluids? Do people use the method of large numbers usually or do they do something else? Is there any disadvantage to the method of large numbers that someone can see? I am assuming that there is probably some standard efficient method used in most commercial and non-commercial codes that solves this problem (obviously I not expecting people to know all the inner workings of every commercial finite element solver, but this problem seems basic/fundamental enough that someone likely has worked on such projects and could provide guidance).

Answer 1

In deal.II (http://www.dealii.org -- disclaimer: I'm one of the principal authors of that library), we do eliminate whole rows and columns, and it is not too expensive overall. The trick is to use the fact that the sparsity pattern is typically symmetric, so you know which rows you need to look into when eliminating a whole column.

The better approach, in my view, is to eliminate these rows and columns in the cell matrices, before they are added to the global matrix. There you work with full matrices, so everything is efficient.

I have never heard of the large-numbers approach and would not use it because surely it will lead to terribly ill-conditioned problems.

For reference, the algorithms we use in deal.II are described conceptually in lectures 21.6 and 21.65 at http://www.math.colostate.edu/~bangerth/videos.html . They closely match your description.

Answer 2

Zero BCs are easy. For Non zero BC's you can also use Lagrange multipliers. E.g., see here. One advantage of LMs is that you can use any constraint equation, although the system becomes indefinite so you need an appropriate solver.