OpenSees Spy

Original Post - 13 Mar 2022 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

A previous post on the fullGenLapack eigenvalue solver led me down a rabbit hole of sorting out all the equation solvers (linear and eigen) and equation numberers available in OpenSees.

I have yet to emerge from the rabbit hole, but I wanted to show how to create “spy” matrices from OpenSees models.

A spy matrix shows the locations of the non-zero matrix entries. From this topology, you can determine the bandwidth, profile/skyline, and sparsity pattern of the matrix. How the equations are numbered can have a huge impact on the matrix bandwidth, and thus the computational time for some equation solvers, as well as the amount of data required to store the matrix.

For instance, the solid bar model shown in this post, using a quadrilateral mesh with mesh size c=t/2, gives the 1200x1200 spy matrices shown below for the Plain, AMD, and RCM numberers available in OpenSees.

Spy matrices using different equation numberers

The figure also shows the matrix bandwidth produced by each numbering scheme. Plain numbering, i.e., numbering the DOFs in the order of user-input node tags, give a bandwidth of 1101. As expected, RCM produces a very narrow band, and is the best option for most solvers available in OpenSees.

The AMD numberer leads to a bandwidth (1200) that is equal to the number of equations–and larger than the Plain bandwidth. This is not good for the BandSPD solver, but not as bad as you might think for BandGeneral and other solvers. So don’t totally write off the AMD numberer.

The intricate relationship between numberers and solvers in OpenSees is what I need to figure out in order to get out of the rabbit hole. If you can lend a hand, that would be great.

Until then, here is the Python code to produce a spy matrix and compute bandwidth from an OpenSees model.

import openseespy.opensees as ops
import numpy as np
import matplotlib.pyplot as plt

# 
# Build your model
#

ops.analyze(1) ;# Or whatever

Neqn = ops.systemSize()

# Build spy matrix
SpyMatrix = np.identity(Neqn)
for e in ops.getEleTags():
   dofs = []
   # Build list of DOFs for this element
   for nd in ops.eleNodes(e):
      dofs += ops.nodeDOFs(nd)

      for idof in dofs:
         if idof < 0: # Constrained DOF
            continue
         for jdof in dofs:
            if jdof < 0: # Constrained DOF
               continue
            SpyMatrix[idof,jdof] = 1.0

# Determine bandwidth
bw = 0
for i in range(Neqn):
   bwi = 0
   # Find non-zero farthest from diagonal on this row
   for j in range(i,Neqn):
      kij = SpyMatrix[i,j]
      if kij != 0.0:
         bwi = j-i+1
   if bwi > bw:
      bw = bwi

# Plot       
plt.spy(SpyMatrix,markersize=.05)
plt.xlabel(f'Bandwidth={bw}')

OpenSees Cloud

OpenSees AMI

OpenSees Spy