Validity of algorithm for assembling the finite element global stiffness matrix

This blog post describes an algorithm for constructing the finite element global stiffness matrix as two basic steps:

[Step 1] Initialize a table with the elements shared by a pair of connected node vertices, as shown below:

[Step 2] For each vertex in each element copy the submatrices from the element stiffness matrix to their appropriate locations in the global stiffness matrix, as shown below:

My question is, why is step 1 necessary? I can map local element vertices (indexed 1,2,3 -- if counting from 1) to (global) nodal points without the table created in step 1. Or is there something I am missing?

ultimately the matrix you create in step 1 is where you're going to construct the global stiffness matrix, which is a direct summation of all the element stiffness matrices.

This global matrix is going to be used to solve a linear system $$A x = b$$, where $$A$$ is the global stiffness matrix you just constructed by "adding" all the element stiffness matrices in step 2 together.

You could initialize a dense matrix (i.e. have memory allocated for every possible $$i,j$$ combination of global indices) and use an implicit node numbering mapping function to satisfy the requirements of step 1, however many FEM problems are quite sparse (i.e. the number of entries in the matrix which are non-zero is small compared to the total number of entries), thus a sparse matrix which tries to allocate memory only for entries which are non-zero (along with some additional information which can be used to determine which entry has a given value).

To efficiently construct sparse matrices you typically want to pre-allocate all the memory it will need before-hand (similar to how you would pre-allocate a dense matrix and then just set the values to what you want). Doing so requires you to have some estimate of which entries in your matrix could ever possibly be non-zero, since it's better to slightly over-allocate than to under-allocate. You can determine this by figuring out which $$i,j$$ global index combinations are a part of any local element stiffness matrix, and simply mark those entries as potentially non-zero, and everything else is guaranteed to be zero. Then you can construct a sparse matrix using this information.

You don't have to take this exact approach necessarily, but this is one common way which finite element methods are implemented in practice because it's relatively easy to understand and efficient, and there are of course additional extensions/modifications you can make to this method.

• So, step 1 is only useful in order to determine the number and positions of non-zero stiffness values, in order to create a sparse matrix? (I have worked with sparse matrices and I understand their benefits.) May 25, 2021 at 1:39
• It can also be used for understanding what data dependencies need to be satisfied when you want to write parallel codes, especially those with non-shared memory like MPI, though I think most parallel sparse matrix libraries can simplify this process significantly. May 25, 2021 at 17:06
• I see how step 1 would be useful in calculating dependencies. I'm actually implementing a small FEA example (just 4 elements) and I'm not looking to do anything fancy. Therefore I can ignore step 1. May 26, 2021 at 3:16