OpenSees Direct Stiffness

Import packages

 import openseespy.opensees as ops
 import numpy as np

Define your model

 ops.wipe()
 ops.model('basic',...)

 ops.node(1,...)
 ops.fix(1,...)
   
 ops.element('elasticBeamColumn',1,...)

Define DOF and CDOF as dictionaries

Note: Number keys sequentially starting at 1. Values are lists of [node,dof] lists

 DOF = {}
 DOF[1] = [[1,3]]
 DOF[2] = [[2,1]]
 ...
 DOF[8] = [[10,3]]

 CDOF = {}
 CDOF[1] = [[1,1]]
 CDOF[2] = [[1,2]]
 ...
 CDOF[5] = [[4,2]]

Find numbers of DOFs and basic forces
```
 Nfdof = len(DOF)
 Ncdof = len(CDOF)
 Nbf = NumBasicForces(whichEleTags=None)
```
Note: Specify whichEleTags as a list of element tags. If unspecified, whichEleTags = ops.getEleTags()

Form load-independent objects

 Af = np.zeros((Nbf,Nfdof))
 FormAf(Af, DOF, whichEleTags=None)

 Ac = np.zeros((Nbf,Ncdof))
 FormAc(Ac, DOF, CDOF, whichEleTags=None)

 Kb = np.zeros((Nbf,Nbf))
 FormKb(Kb, whichEleTags=None)

Note: Use FormAf2, FormAc2, and FormKb2 for cantilever basic system

Form stiffness matrix

Either using objects

 Kff = np.zeros((Nfdof,Nfdof))
 Kff = Af.T@Kb@Af

or from the OpenSees model

 Kff = np.zeros((Nfdof,Nfdof))
 FormKff(Kff, DOF)

The results should be the same

Define a constant time series
```
 ops.timeSeries('Constant',1)
```

Define nodal loads

 ops.pattern('Plain',1,1)
 ops.load(1,...)
 ops.load(4,...)

Form nodal load vector
```
 Pf = np.zeros(Nfdof)
 FormPf(Pf, DOF)
```

Define member loads

 ops.pattern('Plain',2,1)
 ops.eleLoad(2,...)
 ops.eleLoad(3,...)

Form fixed-end basic force vector
```
 Pb0 = np.zeros(Nbf)
 FormPb0(Pb0, DOF, whichEleTags=None)
```
Note: Use FormPb02 for cantilever basic system

Define support displacements

 ops.pattern('Plain',3,1)
 ops.sp(1,2,...)
 ops.sp(4,2,...)

Form support displacement vectors

 Uc = np.zeros(Ncdof)
 FormUc(Uc, CDOF)

 Ubc = np.zeros(Nbf)
 FormUbc(Ubc, DOF, whichEleTags=None)

Note: Use FormUbc2 for cantilever basic system

Solve for nodal displacements (using OpenSees)

 Uf = np.zeros(Nfdof)
 SolveUf(Uf, DOF)

Recover reactions and basic forces

 Pc = np.zeros(Ncdof)
 FormPc(Pc, CDOF)

 Pb = np.zeros(Nbf)
 FormPb(Pb, whichEleTags=None)

Note: Use FormPb2 for cantilever basic system

Form equivalent nodal load vectors

 Pfw = np.zeros(Nfdof)
 Pfw = Pf - Af.T@Pb

 Pcw = np.zeros(Ncdof)
 Pcw = Pc - Ac.T@Pb

Remove nodal and member loads

 ops.remove('loadPattern',1)
 ops.remove('loadPattern',2)

Define lists of axially and flexurally rigid elements
```
 axialRigid = [1,2]
 flexuralRigid = [2]
```

Find number of constraint equations

 Nconstr = NumConstraintEquations(axialRigid, flexuralRigid)

Form constraint matrix and vector

 Cf = np.zeros((Nconstr,Nfdof))
 FormCf(Cf, DOF, axialRigid, flexuralRigid)
   
 R = np.zeros(Nconstr)
 FormR(R, DOF, axialRigid, flexuralRigid)   

Remove support displacements
```
 ops.remove('loadPattern',3)
```