Just because OpenSees is open source does not mean it is a fully transparent box. This is mostly because documentation has lagged behind development. So, pessimists would say the box is semi-opaque while optimists would characterize it as semi-transparent. But a few parts of OpenSees are definitely housed in an opaque box.
Take, for instance, the stiffness matrix. We received many requests for functionality to extract the stiffness matrix for post-processing and offline calculations, e.g., in MATLAB. As a result, the printA function was added a few years ago. This function returns the left-hand side matrix in whatever linear system of equations, 
Note that the printA function only works with the FullGeneral system of equations, which stores the entire 
With no input arguments, printA prints the 
printA -file filename. You can then load this file into MATLAB. In OpenSeesPy, there is an option to return the matrix as a list to the interpreter, printA('-ret'). More on this later.
For a static analysis, printA gives the tangent stiffness matrix. For a dynamic analysis, printA gives the effective tangent stiffness matrix, which is not physically meaningful because it is a linear combination of the mass, damping, and static stiffness matrices with scalar coefficients dictated by the integrator and the time step.
Let’s suppose you’re doing a dynamic analysis and want to get 
printA. Easy. And, I recently learned that you can get printA after switching to the NewmarkExplicit integrator with 
But what if your model has viscous dampers and you want to get the non-classical damping matrix? What if you want to do state space calculations offline because OpenSees does not have state space functionality? You’ll be SOL on the damping matrix. As far as I can tell, there’s no integrator available that you can trick into forming only the damping matrix like you can with NewmarkExplicit and the mass matrix.
So, after trying some arduous and non-scalable approaches to forming printA. I named it GimmeMCK (see PR #343) because better names just didn’t come to mind.
Python
ops.integrator('GimmeMCK', m, c, kt, ki=0.0)
Tcl
integrator GimmeMCK $m $c $kt <$ki>
Assembles A = m*M + c*C + kt*KT + ki*KI for the printA command
m = factor applied to mass matrix
c = factor applied to damping matrix
kt = factor applied to tangent stiffness matrix
ki = factor applied to initial stiffness matrix (optional, default=0.0)
Because you don’t want to perform an actual analysis with this integrator, I removed the functionality to advance time in the domain and to update displacement, velocity, and acceleration like you see in the proper dynamic integrators. Although, maybe one would want this functionality to play around with new integrators without having to touch the source code. Let me know if that functionality appeals to you.
With your model defined, at any time you can switch to the GimmeMCK integrator, dump out the matrices you want, then switch back to a real dynamic integrator and proceed as if nothing happened. Don’t worry if you see analysis failed warnings because the assembled matrix is singular, e.g., massless DOFs when assembling
Below is a synopsis of using the GimmeMCK integrator to get the mass, stiffness, and damping matrices from an OpenSees model defined in Python.
import openseespy.opensees as ops
import numpy as np
import scipy.linalg as slin
#
# Define your model
#
ops.wipeAnalysis()
ops.system('FullGeneral')
ops.analysis('Transient')
# Mass
ops.integrator('GimmeMCK',1.0,0.0,0.0)
ops.analyze(1,0.0)
# Number of equations in the model
N = ops.systemSize() # Has to be done after analyze
# Convert to np array and reshape to NxN matrix
M = np.array(ops.printA('-ret')).reshape((N,N))
# Stiffness
ops.integrator('GimmeMCK',0.0,0.0,1.0)
ops.analyze(1,0.0)
K = np.array(ops.printA('-ret')).reshape((N,N))
# Damping
ops.integrator('GimmeMCK',0.0,1.0,0.0)
ops.analyze(1,0.0)
C = np.array(ops.printA('-ret')).reshape((N,N))
Continue with the following code if you want to perform state space analysis on the model. You’ll notice code to separate the DOFs with mass from those without mass using the nodeDOFs and nodeMass commands. This separation is necessary for static condensation prior to the state space analysis. You may find the node utility commands useful for other stuff too.
# Determine number of DOFs with mass
massDOFs = []
for nd in ops.getNodeTags():
for j in range(NDF): # NDF is number of DOFs/node
if ops.nodeMass(nd,j+1) > 0.0:
massDOFs.append(ops.nodeDOFs(nd)[j])
# Number of DOFs with mass
Nmass = len(massDOFs)
# DOFs without mass
masslessDOFs = np.setdiff1d(range(N),massDOFs)
Nmassless = len(masslessDOFs)
# Form matrices for D*x = -lam*B*x
B = np.zeros((2*Nmass,2*Nmass)) # = [ 0 M; M C]
D = np.zeros((2*Nmass,2*Nmass)) # = [-M 0; 0 K]
# Mass
B[:Nmass,:][:,Nmass:2*Nmass] = M[massDOFs,:][:,massDOFs]
B[Nmass:2*Nmass,:][:,:Nmass] = M[massDOFs,:][:,massDOFs]
D[:Nmass,:][:,:Nmass] = -M[massDOFs,:][:,massDOFs]
# Damping
B[Nmass:2*Nmass,:][:,Nmass:2*Nmass] = C[massDOFs,:][:,massDOFs]
# Static condensation
Kmm = K[massDOFs,:][:,massDOFs]; Kmn = K[massDOFs,:][:,masslessDOFs]
Knm = K[masslessDOFs,:][:,massDOFs]; Knn = K[masslessDOFs,:][:,masslessDOFs]
# Kc = Kmm - Kmn*inv(Knn)*Knm
if Nmassless > 0:
Kc = Kmm - np.dot(Kmn,np.linalg.solve(Knn,Knm))
else:
Kc = K
# Stiffness at DOFs with mass
D[Nmass:2*Nmass,:][:,Nmass:2*Nmass] = Kc
# State space eigenvalue analysis
[lam,x] = slin.eig(D,-B)
Let me know if you run into any problems getting the matrices out of OpenSees.
The title for this post was inspired by ZZ Top’s “Gimme All Your Lovin'”, which saw heavy play on MTV, back when MTV played music videos.
Portwood Digital 
PD,
This is awesome! Just what I needed to perform a response spectrum analysis!
I am definitely going to try it and let you know!
As an extra comment, I just would say that saving a matrix to a text file is probably not the most efficient way. I am working in a model whose system size is 3000 and I found that saving it to a binary file using numpy.save() was the best approach.
LikeLiked by 1 person
Hola Gustavo,
Thanks for the tip! Let me know how it goes.
PD
LikeLike
Hello PD,
I have a question regarding the element’s stiffness matrix. Is there a way to obtain the stiffness matrix of each element of a structure modeled in OpenSees?
Thank you.
LikeLike
Hello Mary,
Some of the frame elements have a recorder option for their basic stiffness matrix. Most elements do not have the option though.
PD
LikeLike
This is a great development! Thank you
LikeLiked by 1 person
Good job! Has this integrator “‘GimmeMCK'” been included in the release of OpenSeesPY? I can’t find this interpreter on the home page of Openseespy. Thank you!
LikeLike
Thanks, Yan!
It should be in OpenSeesPy, but not yet documented on the page.
PD
LikeLike
I just tried this command and it works well. Thanks for this great development. I found one thing confusing though. In the synopsis you provided above, the numberer isn’t specifed. In this case, from the error message I found OpenSees assumes the default RCM numberer. However, even if I specifeied the numberer as ‘Plain’, error message still prints ‘RCM default will be used’, which is confusing. In reality, Plain numberer has been used, because [K] given by Plain numberer is diffrent from the previous [K] given by RCM.
Please correct me if I am wrong. Thanks again!
LikeLike
Thanks! Whether or not you see a warning message depends on the order in which you specify the analysis commands.
LikeLike
I wanted to use this function to calculate the energy of the system.. say for viscous damping, something along the lines of the damping matrix *velocity vector * nodal displacement for a given time instant in the analysis. One problem that I’m trying to figure out is which DOF maps to which velocity vectors /displacement vectors since it seems that equalDOF (and maybe other constraints), changes the size of the MCK matrices. It’s easy enough to get the velocity/disp/acceleration vectors of the system (using recorders on every node), but I wish I could do something to map the MCK to those nodal dofs so I could perform energy calcs.
LikeLike
You can use the nodeDOFs command to get the equation numbers.
LikeLike
Ah- That’s perfect!
Thank you Prof. McKenna.
I briefly saw energy recorders mentioned in a powerpoint somewhere on the web- will those be implemented in the future?
LikeLike
Ha! I’m flattered that you think Frank writes this blog 🙂
Regarding your recorder question, you should create a GitHub issue to see if this feature can be added.
LikeLike
Just wanted to drop a note of appreciation and thanks for your blog. This and several other posts have clarified some very niche topics in Opensees for me. Also, the title on this post is brilliant. Thanks Dr. Scott.
LikeLiked by 1 person
Thank you! It’s nice to know the posts are helpful.
LikeLike
I have defined modal damping with damping ratio of 0.02 in first mode. My model has only one storey and 1 bay. But the extracted damping matrix has 0 as its elements. What could be the problem? Does the integrator only works with Rayleigh damping?
LikeLike
The damping matrix is never formed explicitly with modal damping. The damping forces are added directly to the right-hand side vector.
LikeLike
Thank you very much for the clarification.
LikeLiked by 1 person
I’m sorry, it turns out a necessary function wasn’t defined in GimmeMCK. It’s there now and should return the modal damping matrix.
https://github.com/OpenSees/OpenSees/commit/d2250047f74ab47fe7136e22769fcfb2e8b7c5a3
LikeLike
Prof. Scott Thank you very much for updating the integrator. Extracted modal damping matrix matches with manual calculation.
LikeLiked by 1 person
Dear Prof. Scott, thank you for your post about GImmeMCK. It’s a wonderful function in OpenSees and I have been extensively using it.
I was quite curious about your approach to the state space analysis because the method used in this article is alien to me. AFAIK, State Space requires the K and M only. And usually the A matrix is just = [zeros ones; -k/m -c/m]. Does your D and -B represent -k/m and -c/m, respectively? My intuition about the matrix size said it not (because the A matrix in state space would only be 30×30).
Thank you very much for your kind attention.
LikeLike
Or, to rephrase my question in a simpler term, what kind of state space formulation is [ 0 M; M C] [-M 0; 0 K]?
Thanks again for your attention.
LikeLike
It’s from the Appendix of Chapter 14 in Chopra, 5th edition
LikeLike
Dear Prof. Scott, I learned a lot from GImmeMCK. Thank you so much! I find a question when I am using the GImmeMCK: Why C≠(alphaM *M+betaKinit*K)?
LikeLike
Without any context, that question cannot be answered. I recommend you post more details on the OpenSees message board or Facebook group.
LikeLike
Dear Prof. Scott, I uploaded more details and source files on Facebook. In a two DOF rigid shear frame model, the Rayleigh damping matrix C calculated by ‘GimmeMCK’ and the result calculated by ‘M*alphaM+K*betaKinit’ are the same. However, in a two DOF rigid shear frame model with diagonal braces. I’m very confused about this. Thank you very much!
LikeLike
Thank you very much for your answer in the Facebook group! Now the result is perfect! Can I ask you two more questions: (1) Is the order of the degrees of freedom of the mass matrix obtained by GimmeMCK consistent with the order of getNodeTags0? (2) Does the node mass obtained by the nodeMass0 command include all masses (including the node mass converted from the consistent mass), or only those created with the Mass command?
LikeLike
Hi Prof. Scott,
I am doing an analysis where I applied rigid diaphragm. Because of the rigid diaphragm, the system size is no longer number of node * DOF, thus I am facing some problem with the static condensation code to remove the massless node using the above code. I don’t know how it really work, but I guess OpenSees internal probably has done some static condensation to create the rigid diaphragm.
(Note: it’s a 2D problem so rigid diaphragm may not be exactly the correct term, pardon my ignorance).
I used this code as an alternative to the first 12 lines of your state-space matrix construction code.
# Number of DOFs with mass
massDOFs = np.where(np.diag(M) != 0)[0]
Nmass = len(massDOFs)
# DOFs without mass
masslessDOFs = np.where(np.diag(M) == 0)[0]
Nmassless = len(masslessDOFs)
Do you have any comment on this? For your kind attention, thank you.
LikeLike
I used systemSize in the code, not Nnode*Ndof.
LikeLike
It is for this part
# Determine number of DOFs with mass
massDOFs = []
for nd in ops.getNodeTags():
for j in range(NDF): # NDF is number of DOFs/node
if ops.nodeMass(nd,j+1) > 0.0:
massDOFs.append(ops.nodeDOFs(nd)[j])
For the M,C,K, yes it is determined by the system size. But for the massDOFs, isn’t the loop above equals to Nnode*Ndof? Pardon me if I am reading it wrong.
Thank you for the attention.
LikeLike
Hi Prof. Scott,
In openseespy i try to get mass, stiffness and damping matrices for mdof lumped model which have isolator in midstory. GimmeMCK does not give me right damping matrix but it gives right mass and stiffness matrices. How can i get right damping matrix?
Thanks in advance
LikeLike
Please post your script and question on the OpenSees message board or Facebook group. Also, generate a MWE where you know what the “right” damping matrix is.
(I am removing the long script your posted with your comment)
LikeLike
i did
LikeLike
Dear Prof. Scott, I use OS software to analyze the dynamic response of a single-layer framework model (3-D, 6Dofs, 8 nodes, 4 nodes fix, 4 nodes with lumped mass) under earthquake loads(5000.02s, 10s), first I get the displacement(d, 50024), velocity(v, 50024) and acceleration(a, 50024) of 6 Dofs for the 4 nodes with lumped mass, and I get MCK matrix by the GimmeMCK(M,2424;K,2424;C,2424). Then, according to solving the dynamic equation, Ma+Cv+Kd=-M*g, inverted seismic acceleration. If the mass of the four nodes is the same, the seismic acceleration obtained from each node’s response inversion are the same and consistent with the original input seismic acceleration. However, if the mass of the four nodes is different (with constant stiffness), the seismic waves obtained from each node’s response inversion are different and inconsistent with the original input seismic waves. I want to ask if dynamic analysis in OS is not solved through simple dynamic equations.
LikeLike
I don’t understand what you’re attempting to do, but of course OpenSees solves dynamic equations properly. Double check what you are trying to accomplish.
LikeLike
Dear Prof. Scott, thanks a lot for your reply. In fact, I want to use the extracted MCK matrix to invert the input seismic accelerations based on the complete output dynamic response of each node. However, the seismic acceleration obtained from solving the dynamic equation is not the only value. The input is unique, and the output is also unique.
LikeLike
The matrix used in solving for displacements depends on the time step and the time integration method, e.g., the Newmark gamma and beta factors. And you can’t solve for the response of individual nodes because there is coupling between the nodes due to stiffness.
LikeLike
Dear Prof. Scott
Thank you for providing this blog with this information.
However, I want to build my model using Opensees.
I did not start as I want to confirm that I can postprocess my K,M,C matrices as follows:
I want these matrices to be in terms of notations like [k1+k2 -k2;-k2 k2]
something like that. Because I need later to somehow parametrize each one.
is that possible using the command you have provided?
Thank you
LikeLike
Hello Hassan,
I don’t understand what you mean, but you can try what you’re planning and see if it works.
Michael
LikeLike
I think it is not possible to get whatever I wanted,
I meant to say that If we run PrintA it will give a matrix which contains the assembled stiffness matrix of the system.
for example if we have 2dof system with equal stifnesses say k1=k2=1000.
then it will give us [2000 -1000;-1000 1000]?
not [k1+k2 -k2;-k2 k2].
in my case this want be helpful as i want to scale k1 and k2 individually.
however, i will look to how to do this later.
i have other issue, the PrintA command is not giving the K matrix. nor M after switching to Newmark integrator.
is there any issue I might check?
this is the code i used after specifying the model essential parts.
Run a one step gravity load with no loading (to record eigenvectors)
———————————————————————–
integrator LoadControl 0 1 0 0
Convergence test
tolerance maxIter displayCode
test EnergyIncr 1.0e-10 100 0
Solution algorithm
algorithm Newton
DOF numberer
numberer RCM
Constraint handler
constraints Transformation
System of equations solver
system ProfileSPD
analysis Static
printA -file y
LikeLike
printA only works with the FullGeneral system
LikeLike
Thank you,
Opensees closes when i run the code using FullGeneral system?!
What might be the mistake?
LikeLike
You have to analyze before printA
LikeLike
Thank you very much, it finally worked and I can view the matrix on MATLAB.
However, I could not switch to Newmark integrator with 0 /gamma.
It is not working with static analysis and when I did transient analysis opensees dosn’t run
LikeLike
Dear Dr. Scott,
Thank you for your excellent explanation of this matter! I have been grappling with OpenSees while attempting to get the mass, damping, and stiffness matrix (MCK). I’ve experimented with OpenseePy and this method to extract MCK.
Once I have the MCK matrix, I perform an eigen analysis and compare the resulting frequencies to those obtained from OpenSees/OpenseesPy. However, I’ve encountered an issue with a specific model. The problem is that I’m getting a series of zero frequencies as the first few results. I suspect this might be related to a rigid-link beam element that I have in one of the Degrees of Freedom (DOFs).
To investigate further, I replaced the rigid-link element with a very stiff one, but I still couldn’t match the frequencies with the output from the OpenseePy model. Do you have any insights or suggestions on how I might resolve this issue?
Your guidance would be greatly appreciated.
LikeLike
Sorry, no. Sounds like an issue with the rigid link, but you’ll have to make a minimal example and figure things out.
LikeLike
Dear Prof. Scott,Thank you for develping the GimmeMCK. I see that the example you provided is for OpenSEESpy, may I ask whether this command can be used in Tcl language editor? I tried to use GimmeMCK in tcl, however, the output seems to be finite tangential stiffness, and the mass matrix and stiffness matrix of the structure cannot be obtained.
The follows are my input command in Tcleditor:
integrator GimmeMCK 1.0 0.0 0.0
analyze 1 0.02
puts “The mass matrix is:”
printA
#
integrator GimmeMCK 0.0 0.0 1.0
analyze 1 0.01
puts “The tangential stiffness matrix is:”
printA
integrator GimmeMCK 0.0 0.0 0.0 1.0
analyze 1 0.01
puts “The initial stiffness matrix is:”
printA
#
The outputs are as follows:
The mass matrix is:
2.275e+08 0 2.7487e-10
7.18039e-26 1.33491e+06 -2.40283e+06
2.7487e-10 -2.40283e+06 5.7668e+06
The tangential stiffness matrix is:
2.275e+08 0 2.7487e-10
7.18039e-26 1.33491e+06 -2.40283e+06
2.7487e-10 -2.40283e+06 5.7668e+06
The initial stiffness matrix is:
2.275e+08 0 2.7487e-10
7.18039e-26 1.33491e+06 -2.40283e+06
2.7487e-10 -2.40283e+06 5.7668e+06
LikeLike
I don’t understand your question. Please post modeling questions on the OpenSees Facebook group.
LikeLike
Thank you Prof. Scott, my question is that whether the command ‘GimmeMCK’ can be used in Tcl language editor? I will try to find the OpenSEES Facebook group.
LikeLike
You can try this example in Tcl.
LikeLike
Dear Prof. Scott
Thank you for providing this blog with this information.
I have built a simple 3D concrete frame structure using OpenSeesPy and am preparing to apply a seismic wave input to obtain the structural response of the model under seismic wave excitation. During the nonlinear time history analysis, the structure is theoretically expected to sustain damage, leading to changes in structural stiffness. I would like to ask if GimmeMCK can capture the stiffness matrix of the structure at each time step during the time history analysis, or the stiffness matrix of the structure in a damaged state towards the end of the analysis? I look forward to your reply!
LikeLike
Yes, you can call GimmeMCK inside an analysis loop over number of steps with analyze(1,dt) in the loop.
LikeLike
Thank you, Prof. Scott.
LikeLike