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

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).

• Do you have evidence that this is substantially slowing you down? – Bill Barth Aug 21 '15 at 10:25
• @BillBarth Yes, although there is always the chance that I am doing something inefficiently. Gennadily himself writes that while the explicit method is easy for full 2d arrays, "..it is not always easy to access matrix rows and columns when a matrix is in compact format." suggesting that the explicit method may be more challenging to implement efficiently. As my code is currently written, the explicit method can take more time then the actual solve. – James Aug 21 '15 at 17:34
• do it like Wolfgang says and apply boundary conditions to element matrices before you assemble. – Bill Barth Aug 21 '15 at 21:47
• @BillBarth Yes I think I will do that. His videos are amazing! I just left a comment/question for him about whether you need to re-assemble the global matrices at each timestep, after which I think I will accept his answer. – James Aug 21 '15 at 21:49

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.

• In the case of a time-dependent problem (say the heat equation) do you re-assemble the global matrix at every timestep? The reason I ask is that in the case of non-zero Dirichlet conditions you need information from the original global matrix when modifying the righthand side, but if you zeroed out those columns during the previous timestep then this information is lost (unless you store it in additional arrays). This wouldnt be a problem if the global matrix is re-assembled every timestep though which is what I am considering doing and what would have to be done anyway if using adaptive mesh. – James Aug 21 '15 at 21:47
• It depends on the application. All of the "big" codes solve nonlinear time dependent problems, and for these it's clear that you need to re-assemble one way or the other. For linear codes, you can just store the original matrix, and in every time step copy it somewhere else, apply boundary conditions, and then use that in the solver. This just requires more memory, but is otherwise cheap. – Wolfgang Bangerth Aug 24 '15 at 0:06
• Ah I see that is what I suspected. I will implement as you suggested. Ok thats for your help. Those deallii tutorial videos are really good btw! – James Aug 24 '15 at 0:08

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.