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 DOF in order to get the 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 =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.

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

And here is the stiffness matrix obtained from the `printA`

command after switching the system to `FullGeneral`

.

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?