If you’ve used OpenSees–even if you’re a geotech–you’ve used the force-based element.

When Remo implemented the force-based element, it was the only material nonlinear frame element available in OpenSees (G3 at the time); thus, the original name `nonlinearBeamColumn`

. Only after a standard displacement-based frame element (`dispBeamColumn`

) was added did we change the name from `nonlinearBeamColumn`

to `forceBeamColumn`

. Other material nonlinear formulations have come since.

That the force-based element in OpenSees is based on the non-iterative formulation developed by Neuenhofer and Filippou (1997) is technically true. However, the non-iterative formulation is not the default behavior. By default, the maximum number of internal iterations, `maxIter`

, is set to 10, not 1. You have to use the `-iter`

option to make `maxIter=1`

.

However, there are some strange boundary cases where using `maxIter=1`

can grind the force-based element to a halt, even for a simple model.

Run this example and watch the convergence test churn on the last dozen or so time steps.

```
wipe
model basic -ndm 2 -ndf 3
node 1 0 0; fix 1 1 1 0
node 2 100 0; fix 2 1 1 0
uniaxialMaterial Steel02 1 50 29000 0.005
section Fiber 1 {
patch rect 1 20 1 -10.0 -5.0 10.0 5.0
}
geomTransf Linear 1
element forceBeamColumn 1 1 2 4 1 1 -iter 1 1e-12
integrator LoadControl 500.0
timeSeries Linear 1
pattern Plain 1 1 {
load 2 0 0 1.0
}
test NormDispIncr 1e-6 10 1
analysis Static
analyze 100
```

I chose Tcl because it’s easier to use diagnostic tools like `gprof`

and `valgrind`

with a stand-alone executable. You’ll see the same churn on the last dozen time steps in an equivalent OpenSeesPy model.

Anyway, I set out to determine why `maxIter=1`

causes huge problems for this simple example.

But first, some background on the OpenSees force-based element state determination is necessary. Given basic deformations, coming from the nodal displacements by way of the geometric transformation, the element attempts three internal Newton algorithms:

- Regular Newton-Raphson, always using the current element flexibility.
- Newton-Raphson with the initial element flexibility on the first iteration, then current flexibility on subsequent iterations.
- Newton-Raphson with the initial element flexibility on every iteration.

A few things to note about these algorithms within the force-based element:

- Despite the code comments, the sequence of algorithms is programmed as 1, 3, then 2, i.e., the element will try iteration with the initial flexibility (3) before trying initial then current flexibility (2).
- To give algorithm 3 a fighting chance, the number of internal iterations is set to
`10*maxIter`

for this algorithm only. - If
`maxIter=1`

, algorithms 2 and 3 are the same, i.e., algorithm 2 will never make it to "subsequent iterations".

So, with `maxIter=1`

and the 1-3-2 order of algorithms, algorithm 2, the best backup plan for failure of algorithm 1, never gets enacted. If 10 iterations with constant flexibility fails, one iteration with constant flexibility will also fail.

The failure of the 1-3-2 algorithm sequence in the force-based element state determination sets off subdivisions of those basic deformations coming from the geometric transformation. Four subdivisions are allowed, each cutting the basic deformations by an order of magnitude: 1/10, 1/100, 1/1000, then 1/10000.

If the 1-3-2 algorithm continues failing, the element will chug on the 1/10000 subdivision and the run-time for simple models goes nuts when `maxIter=1`

.

Of course none of this is an issue with the default `maxIter=10`

. But we want the element to work when `maxIter=1`

, right?

The fix is simple:

- Change the order of algorithms from 1-3-2 to 1-2-3 putting the initial flexibility iteration last.
- Allow algorithm 2 a few more iterations, e.g.,
`maxIter+5`

, to do its thing.

I’m not really sure you can call it a fix though. It’s just a more sensible backup plan for failure of algorithm 1, which is inevitable when there is material yielding and `maxIter=1`

. So, the element is still not really non-iterative, but I’m not ready to make a go at Excalibur and all the responsibility that comes with it.

Anyway, before Frank dismisses PR #1168 and calls me a “tosser”, pull my branch and see how the updates affect your IDAs of 2D and 3D reinforced concrete frames using force-based elements with fiber sections of *Concrete23* and *Steel08*.

To get the branch, issue these two commands in your favorite `git`

client:

```
% git checkout -b mhscott-fbc-update master
% git pull git@github.com:mhscott/OpenSees.git fbc-update
```

Then compile OpenSees (Tcl or Py) and run your IDAs. Do the analyses run faster? Do the results change?

If the results do change, your model was fickle to begin with–consider retracting published papers and writing errata. Or, for works in progress, change the narrative. No, seriously, if your results change, I’d love to hear about it.

It’s never too late to fix a problem. And it’s never too late to revert a pull request.

If the `maxIter=10`

default was an IDA ghost containment unit, I’m Walter Peck, the EPA inspector played by William Atherton (who also played a prickish role in Die Hard) shutting down the Ghostbusters. However, unlike Peck, I promise not to *blame you* if all IDA-hell breaks loose with this change in the force-based element.

You do not need element level iterations, you can find the stiffness of the element in one step (not like the above or the mentioned papers). Please look at http://parker.name.tr/ForceBasedInternalForce.zip (run force.m), I know that you know it but you did not really examined it

LikeLike