Choose Your Own Topology

Original Post - 22 Jun 2025 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

I’ve been working on a sparse linear equation solver. Not anything new, just implementing the methods presented by Timothy Davis in Direct Methods for Sparse Linear Systems.

Why? Because I want to learn how sparse matrix solvers work.

So for the sake of discussion, let’s call my solver BennySparse.

Aside from implementing a standalone sparse solver for personal use, I’d like to incorporate BennySparse in OpenSees. I know BennySparse won’t be as fast as industrial strength sparse solvers like Mumps, UmfPack, SuperLU, and SparseSPD, but it should be better than the banded and profile solvers (maybe) and much better than the full general solver (hopefully).

An important part of any sparse solver is symbolic factorization–basically some pre-analysis graph theory to determine matrix fill-in from matrix topology. For finite element models, matrix topology is based on the equation numbers and element connectivity. Although the numeric factorization of the stiffness matrix changes during a nonlinear analysis, the symbolic factorization remains unchanged unless you remove an element or node from the model.

To introduce symbolic factorization and demonstrate fill-in during Cholesky decomposition, Davis shows the topology of some matrix (Figure 4.2 if you’re reading along). The columns are numbered on the diagonal while the dots indicate non-zero off-diagonals.

11x11 sparse matrix topology

In addition to the non-zeros shown above, six more non-zeros will appear in, or fill-in to, the lower triangular matrix resulting from Cholesky factorization. For larger systems, the amount of fill-in can be significant, leading to more time spent on numeric factorizations. Renumbering equations can mitigate fill-in.

But this post is not about equation numbering or symbolic or numeric factorization. No, this post is about how to make OpenSees assemble a stiffness matrix with the topology shown above–or any other specified matrix topology for that matter–so that I can test BennySparse.

Sure, I could write a main() function and manually populate a matrix for BennySparse to gnaw on. But if I shoot in the dark and populate the matrix in such a way that it is not positive definite, Cholesky factorization will fail.

Instead of guessing, I know a 1D model of zero length elements will ensure a positive definite stiffness matrix. Elements attached to a fixed base node establish the diagonal entries while elements connecting unrestrained nodes give the off-diagonal entries from a list of (i,j) tuples.

import openseespy.opensees as ops 
import numpy as np

ops.wipe()
ops.model('basic','-ndm',1,'-ndf',1)

k = 10
ops.uniaxialMaterial('Elastic',1,k)

Neqn = 11
eleTag = 0
ops.node(0,0); ops.fix(0,1)
# Diagonals
for i in range(Neqn):
    ops.node(i+1,0)
    eleTag += 1
    ops.element('zeroLength',eleTag,0,i+1,'-mat',1,'-dir',1)

# Off-diagonals
for i,j in [(1,6),(1,7),(2,3),(2,8),(3,10),(3,11),(4,6),(4,10),(5,8),(5,11),(6,9),(6,10),(7,11),(8,10),(8,11),(10,11)]:
    eleTag += 1
    ops.element('zeroLength',eleTag,i,j,'-mat',1,'-dir',1)

ops.numberer('Plain')
ops.system('FullGeneral')
# ops.system('BennySparse')

ops.analysis('Static','-noWarnings')
ops.analyze(1)

K = np.array(ops.printA('-ret')).reshape((Neqn,Neqn))
print(K)

Note that I’m using FullGeneral so printA will return the stiffness matrix. Also, I’m using the Plain numberer so the stiffness matrix is not permuted. The result matches the topology of the sparse matrix shown above, from Figure 4.2 of Davis’s book.

[[ 30.   0.   0.   0.   0. -10. -10.   0.   0.   0.   0.]
 [  0.  30. -10.   0.   0.   0.   0. -10.   0.   0.   0.]
 [  0. -10.  40.   0.   0.   0.   0.   0.   0. -10. -10.]
 [  0.   0.   0.  30.   0. -10.   0.   0.   0. -10.   0.]
 [  0.   0.   0.   0.  30.   0.   0. -10.   0.   0. -10.]
 [-10.   0.   0. -10.   0.  50.   0.   0. -10. -10.   0.]
 [-10.   0.   0.   0.   0.   0.  30.   0.   0.   0. -10.]
 [  0. -10.   0.   0. -10.   0.   0.  50.   0. -10. -10.]
 [  0.   0.   0.   0.   0. -10.   0.   0.  20.   0.   0.]
 [  0.   0. -10. -10.   0. -10.   0. -10.   0.  60. -10.]
 [  0.   0. -10.   0. -10.   0. -10. -10.   0. -10.  60.]]

The assembled stiffness matrix is symmetric and diagonally dominant. And, for what it’s worth, the matrix is also an M-matrix.

Now I can move forward and verify that BennySparse produces the correct fill-in pattern for Cholesky decomposition, then test BennySparse against other solvers in OpenSees. Stay tuned!

OpenSees Cloud

OpenSees AMI

Choose Your Own Topology