Calculation of the sparsity structure for finite element matrices

Question

Question: What methods are available to accurately and efficiently calculate the sparsity structure of a finite element matrix?

Info: I'm working on a Poisson Pressure Equation solver, using Galerkin's method with a quadratic Lagrange basis, written in C, and using PETSc for sparse matrix storage and KSP routines. To use PETSc efficiently, I need to pre-allocate memory for the global stiffness matrix.

Currently, I am doing a mock assembly to estimate the number of nonzeros per row as follows(pseudocode)

int nnz[global_dim]
for E=1 to NUM_ELTS
  for i=1 to 6
    gi = global index of i 
    if node gi is free
      for j=1 to 6
        gj = global index of j
        if node gj is free 
          nnz[i]++

This, however overestimates nnz because some node-node interactions can occur in multiple elements.

I've considered trying to keep track of which i,j interactions I have found, but I am unsure how to do this without using a lot of memory. I could also loop over the nodes, and find the support of the basis function centered at that node, but then I'd have to search through all elements for each node, which seems inefficient.

I found this recent question, which contained some useful information, especially from Stefano M, who wrote

my advice is to implement it in python or C, applying some graph theoretical concepts, i.e. consider elements in the matrix as edges in a graph and compute the sparsity structure of the adjacency matrix. List of lists or dictionary of keys are common choices.

I'm looking for more details and resources on this. I admittedly don't know much graph theory, and I'm not familiar with all the CS tricks that might be useful(I'm approaching this from the mathematical side).

Thanks!

Answer 1

Your idea of keeping track of which i,j interactions you have found can work, I think that's the "CS trick" that you and Stefano M are referring to. This amounts to constructing your sparse matrix in list of lists format.

Not sure how much CS you have so I apologize if this is already known to you: in a linked list data structure, every entry stores a pointer to the entry after it and the entry before. It's cheap to add and delete entries from, but not as simple to find items in it -- you might have to look through all of them.

So, for each node i, you store a linked list. Then you iterate through all the elements; if you find two nodes i and j are connected, you go look in i's linked list. If j isn't already there, you add it to the list, and likewise add i to j's list. It's easiest if you add them in order.

Once you've populated your list of lists, you now know the number of non-zero entries in each row of the matrix: it's the length of that node's list. This information is exactly what you need to preallocate a sparse matrix in PETSc's matrix data structure. Then you can free your list of lists because you don't need it anymore.

However, this approach assumes that all you have is the list of which nodes each element contains.

Some mesh generation packages -- Triangle for example -- can output not just a list of elements and which nodes they contain, but also a list of every edge in your triangulation. In that case, you run no risk of overestimating the number of non-zero entries: for piecewise linear elements, each edge gives you exactly 2 stiffness matrix entries. You're using piecewise quadratic, so each edge counts for 4 entries, but you get the idea. In that case, you can find the number of non-zero entries per row with one pass through the edge list using an ordinary array.

With that approach, you have to read through an extra big file from the hard disk, which could actually be slower than using the element list if your actual computation isn't that big. Nonetheless, I think it's simpler.

Answer 2

If you specify your mesh as a DMPlex and your data layout as a PetscSection, then DMCreateMatrix() will give you the correctly preallocated matrix automatically. Here are PETSc examples for the Poisson Problem and Stokes Problem.

Answer 3

Upvoted

I personally dont know of any cheap way of doing this so I simply overestimate the number i.e., use a reasonably large value for all rows.

E.g., for a perfectly structured mesh made of linear 8 node hex elements the nnzs per row in both the diagonal and off diagonal blocks are dof*27. For most fully unstructured automatically generated hex meshes the number rarely exceeds dof*54. For linear tets I have never had the need to go beyond dof*30. For some meshes with very badly shaped/low aspect ratio elements you may have to use slightly larger values.

The penalty is that local (on rank) memory consumption is between 2x-5x so you may have to use more compute nodes on your cluster than usual.

Btw I did try using searchable lists but the time taken to determine the sparsity structure was more than the assembly/solve. But my implementation was very simple and did not use information about edges.

The other option is to use routines like DMMeshCreateExodus as shown in this example.

Answer 4

You're looking to enumerate all unique (gi,gj) connections, which suggests placing them all into an (non-duplicating) associative container and then counting its cardinality - in C++ this would be a std::set < std::pair < int,int > >. In your pseudocode, you'd replace "nnz[i]++" with "s.insert[pair(gi,gj)]", and then the final number of nonzeros is s.size(). It should run in O(n-log-n) time, where n is the number of nonzeros.

Since you probably already know the range of possible gi's, you can "splay" the table by the gi index to improve performance. This replaces your set with a std::vector < std::set < int > >. You populate that with "v[gi].insert(gj)", then the total number of nonzeros comes from summing v[gi].size() for all gi's. This should run in O(n-log-k) time, where k is the number of unknowns per element (six for you - essentially a constant for most pde codes, unless you're talking about hp-methods).

(Note - wanted this to be a comment on the selected answer, but was too long - sorry!)

Answer 5

Start from sparse matrix $E^T$ of size elements$\times$dofs. $$ E_{ij}^T = \left\{ \begin{array}{cc} 1 & \mathrm{if\ dof}\ j \in \mathrm{element}\ i \\ 0 & \mathrm{elsewhere} \end{array} \right. $$ Matrix $A = E\, E^T$ has the sparsity pattern you are looking for. Keep in mind that implementing $E^T$ is easier, this is why i defined $E^T$ instead of $E$.