OpenSees Cloud

OpenSees AMI

The Worst that Could Happen

07 Jun 2026 - Michael H. Scott


A lot can go wrong in a nonlinear structural analysis. You could have an ill-defined model or incompatible analysis options, among many other scenarios.

OpenSees does not do a great job of telling you something’s wrong. Usually you get a cryptic output message like “Failed to solve” or “Failed to get compatible deformations”. In some cases, the software exits or the kernel dies.

These failures at least make it clear something went wrong.

But the worst thing that could happen is there’s absolutely no output message and the analysis continues as if nothing’s wrong. Then you plot the output and see strange, sometimes physically impossible, results.

I came across such a case when using OpenSees to solve the examples in Ann Jeffers’s new book Linear and Nonlinear Methods of Matrix Structural Analysis. The fixed-fixed beam shown below is analyzed in section 10.3.4 of the book.

The load-displacement results shown in the book are based on the following modeling and analysis choices:

Other than recommending a maximum of 100 Newton-Raphson iterations per load step (10 iterations is plenty), there is nothing obviously wrong with this example. Overall, the book contains many good examples that are suited for verification of OpenSees analyses.

Anyway, using OpenSees to run the analysis described above gives the following load-displacement response.

How did we get to the maximum load of P=200 kip, which is physically impossible for this model? There were no warnings or errors. Here’s the script so you can try it yourself with OpenSeesPy version 3.8.0.0 or earlier.

import openseespy.opensees as ops

kip = 1
inch = 1

ksi = kip/inch**2

L = 120*inch
b = 4*inch
h = 12*inch

E = 30000*ksi
sigY = 50*ksi
Mp = b*h*h/4*sigY

Nfibers = 12

Pmax = 200*kip
dP = 1*kip
Nsteps = int(Pmax/dP)

ops.wipe()
ops.model('basic','-ndm',2,'-ndf',3)

ops.node(1,0,0); ops.fix(1,1,1,1)
ops.node(2,L,0)
ops.node(3,2*L,0)
ops.node(4,3*L,0); ops.fix(4,1,1,1)

ops.uniaxialMaterial('ElasticPP',1,E,sigY/E)
ops.section('Fiber',1)
ops.patch('rect',1,Nfibers,1,-h/2,-b/2,h/2,b/2)

ops.beamIntegration('Lobatto',1,1,5)

ops.geomTransf('Linear',1)

ops.element('forceBeamColumn',1,1,2,1,1)
ops.element('forceBeamColumn',2,2,3,1,1)
ops.element('forceBeamColumn',3,3,4,1,1)

ops.timeSeries('Linear',1)
ops.pattern('Plain',1,1)
ops.load(3,0,-1.0,0)

ops.test('NormUnbalance',1e-5,10,0)
ops.algorithm('Newton')
ops.integrator('LoadControl',dP)
ops.analysis('Static','-noWarnings')

for i in range(Nsteps):
    ok = ops.analyze(1)
    if ok < 0:
        break

As expected, the stiffness matrix at multiple sections becomes singular when the sections reach their plastic moment capacity, i.e., all fibers have yielded and return zero for their tangent modulus.

The problem is we do not detect and act on a failure to invert the section stiffness matrix to flexibility in the SectionForceDeformation::getSectionFlexibility() method.

//
// Code is shortened for clarity of this post
//
const Matrix &
SectionForceDeformation::getSectionFlexibility()
{
   const Matrix &ks = this->getSectionTangent();

   ks.Invert(fs);

   return fs;
}

The Matrix::Invert() method will return a negative integer if the matrix inversion failed. We do not even store the return value! In this example, the force-based element goes off and does strange things. You will get the same response with the mixedBeamColumn element, so it’s not a specific element implementation issue.

Ultimately, I don’t think we should call exit(-1) if inversion of the section stiffness matrix fails, so I only added an error message in PR #1759.

//
// Code is shortened for clarity of this post
//
const Matrix &
SectionForceDeformation::getSectionFlexibility()
{
   const Matrix &ks = this->getSectionTangent();

   int ok = ks.Invert(fs);
   if (ok < 0) {
      opserr << "ERROR - failed to invert section stiffness matrix" << endln;
   }

   return fs;
}

Now OpenSees at least tells you something went wrong in the analysis.