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