Is It Close Enough?

Original Post - 05 Nov 2024 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

The locations and weights for Gauss-Lobatto beam integration, the de facto beam integration for force-based elements, are hard-coded in the OpenSees source code. For most cases in the LobattoBeamIntegration class, the locations and weights are written to only ten significant figures instead of 16 or more.

Although I am certain that leaving six sig-figs on the table is the least of your concerns when running IDAs of reinforced concrete shear walls with regularized material properties and what not, those six sig-figs will cause problems for simple tests and verifications–things you should do before all those IDAs.

Suppose you wanted to verify that the free end deflection of an end-loaded cantilever, modeled with a force-based element with elastic sections, matches the exact solution PL³/(3EI).

A test script using five Gauss-Lobatto points and concluding with an assertion of logical comparison between the computed and exact deflections is shown below.

import openseespy.opensees as ops

E = 29000
L = 48
A = 20
I = 800

P = 20

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.geomTransf('Linear',1)

ops.section('Elastic',1,E,A,I)
Np = 5
ops.beamIntegration('Lobatto',1,1,Np)

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

ops.timeSeries('Constant',1)
ops.pattern('Plain',1,1)
ops.load(2,0,P,0)

ops.analysis('Static','-noWarnings')
ops.analyze(1)

u = ops.nodeDisp(2,2)
uexact = P*L**3/(3*E*I)

assert u == uexact

The assertion will fail due to round-off from the truncated hard-coded values in the OpenSees implementation of the Gauss-Lobatto integration rule. The computed deflection is 0.03177931033115315 while the exact deflection is 0.03177931034482759, with differences showing up in the tenth sig-fig, as expected.

A better assertion uses the isclose function from Python’s math library.

import openseespy.opensees as ops

from math import isclose

E = 29000
L = 48
A = 20
I = 800

P = 20

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.geomTransf('Linear',1)

ops.section('Elastic',1,E,A,I)
Np = 5
ops.beamIntegration('Lobatto',1,1,Np)

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

ops.timeSeries('Constant',1)
ops.pattern('Plain',1,1)
ops.load(2,0,P,0)

ops.analysis('Static','-noWarnings')
ops.analyze(1)

u = ops.nodeDisp(2,2)
uexact = P*L**3/(3*E*I)

#assert u == uexact
assert isclose(u,uexact)

As described in its documentation, isclose uses relative and absolute tolerances to determine if two values are close. Try the above script with four instead of five Lobatto points and you’ll see that you need to increase the relative tolerance for isclose from its default value to 1e-8.

assert isclose(u,uexact,rel_tol=1e-8)

Note that to compare computed values to exact solutions of zero, you have to set an absolute tolerance when using isclose. For example, the computed vertical reaction at the free end of the cantilever is 1.7763568394002505e-14, not 0.0. So, set a non-zero absolute tolerance, for which 1e-12 is fine.

ops.reactions()
R = ops.nodeReaction(2,2)

assert isclose(R,0.0,abs_tol=1e-12)

OpenSees Cloud

OpenSees AMI

Is It Close Enough?