Closing the Loop on Direct Assembly

Original Post - 30 Sep 2020 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

All structural engineering students learn the direct assembly method, where you fix all degrees of freedom (DOFs) in a structural model, then impose a unit value of displacement at and in the direction of the j^th DOF in order to get the j^th column of the stiffness matrix from the fixed-end forces of each member. Repeat this process for every DOF and you get the entire stiffness matrix. Do it one more time with the member loads and you get the equivalent nodal load vector.

The next logical step in structural analysis education is to program the stiffness method. However, direct assembly is not exactly how you implement the stiffness method, so it’s understandable that these topics are thought of as independent and unrelated, except for their tedium–one involves a lot of deflected shapes and handwriting while the other requires a lot of coding and debugging.

Sometimes direct assembly and programming the stiffness method are covered in the same course. Other times, for a variety of reasons, these topics are covered in separate classes with several academic terms intervening. This temporal distance only reinforces the notion that the topics are unrelated.

You can use direct assembly to “spot check” the reactions from a computerized structural analysis. But going the other way, have you ever contemplated how to make structural analysis software mimic hand-written direct assembly? You should be able to do this meta-analysis with any structural analysis software, but this post will show how to do it with OpenSees.

First, define your model. Make sure it is first order, linear-elastic. We will impose displacements and rotations equal to 1.0, which is certain to cause convergence problems in a nonlinear model. For demonstration, here’s a one-story, one-bay frame with a diagonal brace (not designed, just picked some shapes). All members have E=29000 ksi. There are 8 unconstrained DOFs in the model: one rotation at each pin support (nodes 1 and 4) and two displacements and one rotation at each of nodes 2 and 3.

Example frame structure with diagonal brace

Next, define static analysis options. The analysis will impose a lot of single point constraints, so use the Transformation constraint handler. Use UmfPack as the solver because there will be no unconstrained DOFs after we impose all those constraints. If you don’t care about the DOF numbering, use the RCM numberer. If you want OpenSees to number the DOFs in the order of your node tags, use the Plain numberer. Use load control with \(\Delta t\)=0 because pseudo-time is not important here. This might be the only acceptable use of the Linear algorithm.

#
# Define your model
#

ops.constraints('Transformation')
ops.system('UmfPack')
ops.numberer('Plain') # or ops.numberer('RCM')
ops.integrator('LoadControl',0.0)
ops.algorithm('Linear')
ops.analysis('Static')

To form the stiffness matrix, first initialize K based on the system size. Then, loop through the model DOFs, and, inside a load pattern, impose a unit value of displacement at and in the direction of the current DOF and zero at the other DOFs. After imposing all the single point constraints, do an analysis and find the reactions. Finally, make another loop through the nodes and assemble the stiffness matrix from the reactions at the all of the nodes. Note that the same constant time series is used over and over.

ops.analyze(1) # To make systemSize return a non-zero
N = ops.systemSize()
K = np.zeros((N,N))

# Store the tags so we don't keep reconstructing this list in the loop
nodeTags = ops.getNodeTags()

ops.timeSeries('Constant',1)
for j in range(N):
   ops.pattern('Plain',1,1)
   for nd in nodeTags:
      dof = 0 # Nodal DOF counter
      for i in ops.nodeDOFs(nd):
         dof += 1
         if i < 0: # A boundary condition
            continue
         if j == i:
            ops.sp(nd,dof,1.0)
         else:
            ops.sp(nd,dof,0.0)

   ops.analyze(1)
   ops.reactions()

   ops.remove('loadPattern',1)
   ops.analyze(1) # To reset the equation numbers

   for nd in nodeTags:
      dof = 0 # Nodal DOF counter
      for i in ops.nodeDOFs(nd):
         dof += 1
         if i < 0: # A boundary condition
            continue
         K[i,j] += ops.nodeReaction(nd,dof)

That second analyze(1) inside the loop resets the equation numbers returned by nodeDOFs to their original values. Setting the single point constraints in the first inner loop will make all the equation numbers -1. So, if we don’t undo that damage, nothing will be assembled into K in the second inner loop. However, you can avoid this extra analysis by mapping equation numbers to nodal DOFs in a Python dictionary instead of relying on calls to nodeDOFs. And with a dictionary, you can also impose your own equation numbers on the model, but that goes beyond the scope of this post.

At any rate, using the above code, here is the stiffness matrix assembled for the frame.

Output of direct assembly procedure

And here is the stiffness matrix obtained from the printA command after switching the system to FullGeneral.

Output of printA command

The meta-analysis and printA give the same result. You can apply the same meta-approach to forming the equivalent nodal load vector: fix all the DOFs, apply the member loads, then assemble the vector from the nodal reactions.

As you might expect, this meta-analysis does not scale well for large models because of all the analyses, reactions (most of which are not needed), and other function calls inside the loop. You’re better off to get the stiffness matrix out of OpenSees using printA. But where’s the sport in that?

OpenSees Cloud

OpenSees AMI

Closing the Loop on Direct Assembly