Non-linear static analysis with pyfe3d.solver#
The geometrically non-linear analyses of Analysis receive the
callables of a model, see Non-linear static analyses. A simply supported
plate with immovable edges under a uniform pressure, whose linear
deflection is about three times the thickness, such that the membrane
action stiffens the plate:
def plate(name, nx=11, ny=11):
"""Simply supported square plate with immovable edges, uniform pressure
The pressure gives a linear deflection of about three times the thickness.
"""
cls, probecls, datacls, num_nodes = SHELLS[name]
a = b = 1.
h = 0.01
E = 70.e9
nu = 0.3
D = E*h**3/(12*(1 - nu**2))
q = 3*h*D/(0.00406*a**4)
prop = isotropic_plate(thickness=h, E=E, nu=nu)
xtmp = np.linspace(0, a, nx)
ytmp = np.linspace(0, b, ny)
xmesh, ymesh = np.meshgrid(xtmp, ytmp)
ncoords = np.vstack((xmesh.T.flatten(), ymesh.T.flatten(),
np.zeros(nx*ny))).T
x = ncoords.flatten()
nids = np.arange(nx*ny).reshape(nx, ny)
n1s = nids[:-1, :-1].flatten()
n2s = nids[1:, :-1].flatten()
n3s = nids[1:, 1:].flatten()
n4s = nids[:-1, 1:].flatten()
if num_nodes == 3:
connectivity = ([(n1, n2, n3) for n1, n2, n3 in zip(n1s, n2s, n3s)]
+ [(n1, n3, n4) for n1, n3, n4 in zip(n1s, n3s, n4s)])
else:
connectivity = list(zip(n1s, n2s, n3s, n4s))
data = datacls()
probe = probecls()
elements = []
for k, nodes in enumerate(connectivity):
elem = cls(probe)
for i, n in enumerate(nodes):
setattr(elem, 'n%d' % (i + 1), n)
setattr(elem, 'c%d' % (i + 1), DOF*n)
elem.init_k_KC0 = k*data.KC0_SPARSE_SIZE
elem.init_k_KCNL = k*data.KCNL_SPARSE_SIZE
elem.init_k_KG = k*data.KG_SPARSE_SIZE
elem.update_rotation_matrix(x)
elem.update_probe_xe(x)
elements.append(elem)
N = DOF*nx*ny
X = ncoords[:, 0]
Y = ncoords[:, 1]
edge = (np.isclose(X, 0) | np.isclose(X, a) | np.isclose(Y, 0)
| np.isclose(Y, b))
bk = np.zeros(N, dtype=bool)
bk[0::DOF] = edge
bk[1::DOF] = edge
bk[2::DOF] = edge
# NOTE drilling rotations do not take part in the response of a flat plate
bk[5::DOF] = True
bu = ~bk
# uniform pressure lumped with the tributary areas of the nodes
wx = np.full(nx, a/(nx - 1))
wx[[0, -1]] /= 2
wy = np.full(ny, b/(ny - 1))
wy[[0, -1]] /= 2
fext = np.zeros(N)
fext[2::DOF] = q*np.outer(wx, wy).flatten()
center = nids[nx//2, ny//2]
return elements, data, prop, x, N, bu, fext, lambda u: u[DOF*center + 2]/h
The class below wraps a list of elements into the four callables, with the
vectors and matrices of all the dofs: calc_fint with update_fint(...,
nonlinear=1), calc_kC with KC0 and KCNL, calc_kG with
KG. The scaling S changes the units of the rotations, which leaves
Crisfield’s convergence criterion unchanged:
class Model(object):
r"""Callables of :class:`pyfe3d.solver.Analysis` for a list of pyfe3d
elements, with the vectors and matrices of all the degrees-of-freedom
``calc_kC`` gives ``KC0`` when ``NLgeom=False`` and ``KC0 + KCNL(u)``
otherwise, ``calc_kG`` gives ``KG(u)`` and ``calc_fint`` the internal
force vector of ``update_fint(..., nonlinear=1)``, whose exact Jacobian
is ``KC0 + KCNL + KG``. ``S`` scales the degrees-of-freedom, the
callables then working with ``c = S u``, the forces ``S^{-1} f`` and the
matrices ``S^{-1} K S^{-1}``. ``elem_kwargs``, e.g. the
``hgfactor_*`` of Quad4R, are passed to ``update_KC0`` and
``update_fint`` of every element.
"""
def __init__(self, elements, data, prop, x, N, fext, S=None,
elem_kwargs=None):
self.elements = elements
self.elem_kwargs = {} if elem_kwargs is None else dict(elem_kwargs)
self.prop = prop
self.x = x
self.N = N
self.fext = fext
self.S = np.ones(N) if S is None else S
num_elements = len(elements)
self.data = data
self.KC0r = np.zeros(data.KC0_SPARSE_SIZE*num_elements, dtype=INT)
self.KC0c = np.zeros(data.KC0_SPARSE_SIZE*num_elements, dtype=INT)
KC0v = np.zeros(data.KC0_SPARSE_SIZE*num_elements, dtype=DOUBLE)
self.KCNLr = np.zeros(data.KCNL_SPARSE_SIZE*num_elements, dtype=INT)
self.KCNLc = np.zeros(data.KCNL_SPARSE_SIZE*num_elements, dtype=INT)
self.KGr = np.zeros(data.KG_SPARSE_SIZE*num_elements, dtype=INT)
self.KGc = np.zeros(data.KG_SPARSE_SIZE*num_elements, dtype=INT)
for elem in elements:
elem.update_probe_xe(x)
elem.update_KC0(self.KC0r, self.KC0c, KC0v, prop,
**self.elem_kwargs)
self.KC0 = self.sparse(self.KC0r, self.KC0c, KC0v)
self.fint_calls = []
def sparse(self, r, c, v):
return coo_matrix((v, (r, c)), shape=(self.N, self.N)).tocsr()
def scaled(self, A):
Sinv = 1./self.S
return csr_matrix(A.multiply(Sinv[:, None]).multiply(Sinv[None, :]))
def internal_forces(self, u):
fint = np.zeros(self.N)
for elem in self.elements:
# NOTE the probe is shared amongst the elements
elem.update_probe_xe(self.x)
elem.update_probe_ue(u)
elem.update_fint(fint, self.prop, nonlinear=1,
**self.elem_kwargs)
return fint
def KCNL(self, u):
v = np.zeros(self.data.KCNL_SPARSE_SIZE*len(self.elements))
for elem in self.elements:
elem.update_probe_xe(self.x)
elem.update_probe_ue(u)
elem.update_KCNL(self.KCNLr, self.KCNLc, v, self.prop)
return self.sparse(self.KCNLr, self.KCNLc, v)
def KG(self, u):
v = np.zeros(self.data.KG_SPARSE_SIZE*len(self.elements))
for elem in self.elements:
elem.update_probe_xe(self.x)
elem.update_probe_ue(u)
elem.update_KG(self.KGr, self.KGc, v, self.prop)
return self.sparse(self.KGr, self.KGc, v)
# NOTE callables of pyfe3d.solver.Analysis
def calc_fext(self, inc=1., silent=True):
return inc*self.fext/self.S
def calc_fint(self, c, silent=True):
self.fint_calls.append(c.copy())
return self.internal_forces(c/self.S)/self.S
def calc_kC(self, c=None, NLgeom=False, silent=True):
if c is None or not NLgeom:
return self.scaled(self.KC0)
return self.scaled(self.KC0 + self.KCNL(c/self.S))
def calc_kG(self, c=None, NLgeom=False, silent=True):
return self.scaled(self.KG(c/self.S))
def analysis(self, bu, uk=None, **kwargs):
an = Analysis(self.calc_fext, self.calc_fint, self.calc_kC,
self.calc_kG, bu=bu, uk=uk)
for key, value in kwargs.items():
setattr(an, key, value)
return an
def zetas(self, bu):
r"""Crisfield's `\zeta` of every call of ``calc_fint``, at the load
factor 1, over the free degrees-of-freedom"""
w = crisfield_weights(self.scaled(self.KC0)[bu, :][:, bu])
out = []
for c in self.fint_calls:
fint = (self.internal_forces(c/self.S)/self.S)[bu]
fext = (self.fext/self.S)[bu]
out.append(crisfield_zeta(fint - fext, fext, fint, w))
return out
The full load in a single increment, with the free dofs given by bu:
the solution is the one of a Newton-Raphson loop written by hand with the
element matrices, and the convergence criterion falls quadratically:
@pytest.mark.parametrize('kind, name', [('plate', 'Quad4R'),
('plate', 'Tria3DSG'),
('beam', 'BeamC')])
def test_newton_raphson_reproduces_the_hand_written_loop(kind, name):
r"""A single increment at the full load, as the hand-written loop, gives
the same solution, with the quadratic convergence of `\zeta`"""
elements, data, prop, x, N, bu, fext, deflection = (
plate(name) if kind == 'plate' else beam(name))
u_linear, u_ref, residuals = solve_by_hand(elements, data, prop, x, N,
bu, fext)
model = Model(elements, data, prop, x, N, fext)
an = model.analysis(bu, initialInc=1., relTOL=1.e-9)
increments, cs = an.static(NLgeom=True, silent=True)
assert increments == [1.]
u = cs[-1]
assert u.shape == (N,)
assert np.all(u[~bu] == 0)
print(name, 'w/h', deflection(u), 'hand-written', deflection(u_ref))
assert np.isclose(deflection(u), deflection(u_ref), rtol=1.e-7)
np.testing.assert_allclose(u, u_ref, rtol=1.e-6,
atol=1.e-8*np.abs(u_ref).max())
zetas = model.zetas(bu)
print('zeta', ' '.join('%.1e' % z for z in zetas))
assert zetas[-1] <= 1.e-9
assert len(zetas) <= 8
# NOTE quadratic decay of zeta in the asymptotic range
for z0, z1 in zip(zetas[:-1], zetas[1:]):
if z0 < 1.e-2:
assert z1 < max(10*z0**1.8, 1.e-11), zetas
Prescribed displacements#
The central deflection of the load-controlled analysis, prescribed with
uk at a constrained dof, with the Newton-Raphson and the Riks methods:
the displacement field is the same and the reaction equals the load:
def point_loaded_plate():
r"""The plate of ``plate('Quad4R')`` with a central point load whose
linear deflection is three times the thickness, instead of the
pressure"""
elements, data, prop, x, N, bu, fext, deflection = plate('Quad4R')
h = 0.01
center = int(np.argmin((x[0::3] - 0.5)**2 + (x[1::3] - 0.5)**2))
dof = DOF*center + 2
f = np.zeros(N)
f[dof] = 1.
model = Model(elements, data, prop, x, N, f)
w1 = solve(model.KC0[bu, :][:, bu], f[bu])[np.flatnonzero(bu) == dof][0]
P = 3*h/w1
model.fext = P*f
return model, bu, dof, P
@pytest.mark.parametrize('NL_method', ['NR', 'arc_length_riks'])
def test_prescribed_displacement_and_reaction(NL_method):
r"""Displacement control: the central deflection of the load-controlled
analysis, prescribed with ``uk``, gives the same displacement field and
a reaction equal to the load"""
model, bu, dof, P = point_loaded_plate()
an = model.analysis(bu, initialInc=0.2, relTOL=1.e-10)
increments, cs = an.static(NLgeom=True, silent=True)
u_load = cs[-1]
w_c = u_load[dof]
print('load control: w/h = %.6f' % (w_c/0.01))
bu_d = bu.copy()
bu_d[dof] = False
uk = np.zeros(model.N)
uk[dof] = w_c
model.fext = np.zeros(model.N)
an_d = model.analysis(bu_d, uk=uk, initialInc=0.2, relTOL=1.e-10,
NL_method=NL_method)
increments_d, cs_d = an_d.static(NLgeom=True, silent=True)
assert np.isclose(increments_d[-1], 1.)
for lbd, u in zip(increments_d, cs_d):
assert np.isclose(u[dof], lbd*w_c, rtol=1.e-14)
u_disp = cs_d[-1]
np.testing.assert_allclose(u_disp, u_load, rtol=1.e-6,
atol=1.e-8*np.abs(u_load).max())
reaction = model.internal_forces(u_disp)[dof]
print('reaction %.8e, load %.8e' % (reaction, P))
assert np.isclose(reaction, P, rtol=1.e-6)
# NOTE the reaction grows faster than the deflection: membrane
# stiffening
reactions = [model.internal_forces(u)[dof]/lbd
for lbd, u in zip(increments_d, cs_d)]
assert np.all(np.diff(reactions) > 0)
Configuration-dependent load#
The aerodynamic load of the piston theory, \(\lambda \beta
\boldsymbol{K}_{A_\beta} \boldsymbol{u}\), from update_KA_beta of Quad4R, enters
calc_fint with the load factor inc, and its unsymmetric load
stiffness enters calc_kC; calc_fext(c=...) gives the load vector of
the current configuration:
def piston_theory_plate():
r"""Quad4R plate under pressure with the piston-theory aerodynamic load
`\boldsymbol{F}_a = \beta \boldsymbol{K}_{A_\beta} \boldsymbol{u}`, scaled to 20 % of the pressure
load in the linear solution"""
elements, data, prop, x, N, bu, fext, deflection = plate('Quad4R')
r = np.zeros(data.KA_BETA_SPARSE_SIZE*len(elements), dtype=INT)
c = np.zeros(data.KA_BETA_SPARSE_SIZE*len(elements), dtype=INT)
v = np.zeros(data.KA_BETA_SPARSE_SIZE*len(elements), dtype=DOUBLE)
for k, elem in enumerate(elements):
elem.init_k_KA_beta = k*data.KA_BETA_SPARSE_SIZE
elem.update_probe_xe(x)
elem.update_KA_beta(r, c, v)
KA = coo_matrix((v, (r, c)), shape=(N, N)).tocsr()
model = Model(elements, data, prop, x, N, fext)
u_lin = np.zeros(N)
u_lin[bu] = solve(model.KC0[bu, :][:, bu], fext[bu])
beta = 0.2*np.linalg.norm(fext[bu])/np.linalg.norm((KA @ u_lin)[bu])
return model, bu, beta*KA, deflection
def test_piston_theory_load_unsymmetric_tangent():
r"""The aerodynamic load is configuration-dependent, `\boldsymbol{F}_a(\boldsymbol{u}) =
\lambda \beta \boldsymbol{K}_{A_\beta} \boldsymbol{u}` with the load factor `\lambda`: it
enters ``calc_fint`` with ``inc`` and its load stiffness `-\lambda \beta
\boldsymbol{K}_{A_\beta}`, unsymmetric, enters ``calc_kC``. Newton-Raphson
converges quadratically, the residual of the full model vanishes, and
the linear static analysis with ``calc_fext(c=...)`` gives the solution
of `(\boldsymbol{K}_{C0} - \beta \boldsymbol{K}_{A_\beta}) \boldsymbol{u} = \boldsymbol{F}`"""
model, bu, BKA, deflection = piston_theory_plate()
from pyfe3d.solver import is_symmetric
assert not is_symmetric((model.KC0 - BKA)[bu, :][:, bu])
zetas = []
def calc_fint(c, inc=1., silent=True):
f = model.internal_forces(c) - inc*(BKA @ c)
R = (f - inc*model.fext)[bu]
w = crisfield_weights(model.KC0[bu, :][:, bu])
zetas.append((inc, crisfield_zeta(R, inc*model.fext[bu], f[bu], w)))
return f
def calc_kC(c=None, NLgeom=False, inc=1., silent=True):
if c is None or not NLgeom:
return model.KC0
return model.KC0 + model.KCNL(c) - inc*BKA
def calc_fext(inc=1., silent=True, c=None):
if c is None:
return inc*model.fext
return inc*(model.fext + BKA @ c)
an = Analysis(calc_fext, calc_fint, calc_kC, model.calc_kG, bu=bu)
an.initialInc = 1.
an.relTOL = 1.e-10
increments, cs = an.static(NLgeom=True, silent=True)
assert increments == [1.]
u = cs[-1]
R = (model.internal_forces(u) - (model.fext + BKA @ u))[bu]
assert np.linalg.norm(R) <= 1.e-8*np.linalg.norm(model.fext[bu])
z = [zeta for inc, zeta in zetas]
print('zeta', ' '.join('%.1e' % zi for zi in z))
assert len(z) <= 8
for z0, z1 in zip(z[:-1], z[1:]):
if z0 < 1.e-2:
assert z1 < max(10*z0**1.8, 1.e-11), z
# NOTE the aerodynamic load changes the solution
model_dead = Model(model.elements, model.data, model.prop, model.x,
model.N, model.fext)
an_dead = model_dead.analysis(bu, initialInc=1., relTOL=1.e-10)
u_dead = an_dead.static(NLgeom=True, silent=True)[1][-1]
assert abs(deflection(u) - deflection(u_dead)) > 1.e-3*deflection(u_dead)
# NOTE geometrically linear: the load correction iterations of the linear
# static analysis
an_lin = Analysis(calc_fext, None, calc_kC, None, bu=bu)
an_lin.relTOL = 1.e-12
an_lin.maxNumIter = 100
u_lin = an_lin.static(silent=True)[1][-1]
ref = np.zeros(model.N)
ref[bu] = solve((model.KC0 - BKA)[bu, :][:, bu], model.fext[bu])
np.testing.assert_allclose(u_lin, ref, rtol=1.e-8,
atol=1.e-10*np.abs(ref).max())
Limit points with the arc-length methods#
The shallow von Mises truss of tests/analytic_problems.py, loaded
through a linear spring beyond its limit load: the arc-length methods trace the
snap-through, along which the load factor becomes negative, until the load
factor 1.0 on the inverted configuration:
def von_mises_truss(frac=0.99, ks=10.):
"""Shallow von Mises truss (height h=1) loaded through a linear spring
DOF 0 is the apex deflection, DOF 1 the loaded point, connected to the
apex through a spring of stiffness ``ks``. The truss softens and reaches
a limit point at ``w = 1 - 1/sqrt(3)``; the load is ``frac`` of the
limit load.
"""
Flim = 2/(3*np.sqrt(3))
fext = np.array([0., frac*Flim])
kC = np.array([[2. + ks, -ks], [-ks, ks]])
def fint(c):
w = c[0]
return np.array([w**3 - 3*w**2 + 2*w + ks*(c[0] - c[1]),
ks*(c[1] - c[0])])
def kG(c):
w = c[0]
return np.array([[3*w**2 - 6*w, 0.], [0., 0.]])
return Problem(fext, fint, kC, kG)
@pytest.mark.parametrize('method', METHODS)
def test_arc_length_snap_through(method):
"""Load of 1.5 times the limit load of the von Mises truss
The path goes through the limit point, the load factor becomes negative
during the snap-through, and it increases again until the load factor
reaches 1.0 on the inverted configuration.
"""
p = von_mises_truss(frac=1.5)
p.analysis(NL_method=method, initialInc=0.05)
p.an.static(NLgeom=True, silent=True)
check_equilibrium(p)
inc = np.asarray(p.an.increments)
w = np.array([c[0] for c in p.an.cs])
assert inc[-1] == 1.0
# limit load factor is 1/1.5, the minimum along the path is -1/1.5
assert inc[inc < 1.].max() < 1/1.5 + 1e-6
assert inc[inc < 1.].max() > 0.6
assert inc.min() < -0.5
assert np.all(inc[:-1] < 1.)
# continuous path along the apex deflection, without jumps
assert np.all(np.diff(w) > 0)
assert np.diff(w).max() < 0.5
assert w[-1] > 2.
np.testing.assert_allclose(p.an.cs[-1],
p.solve_reference(c0=p.an.cs[-1], tol=1e-12),
rtol=1e-6)
The previous converged point is at the distance of the arc length from the last one, and the iterations of a step may converge to it, especially with long steps or near a limit point: such a step, whose increment points against the tangent predictor, is repeated with a shorter arc length, such that the path is not traced backwards:
@pytest.mark.parametrize('convergence, relTOL', [('relative', 1.e-6),
('crisfield', 1.e-9)])
@pytest.mark.parametrize('method', METHODS)
def test_arc_length_no_doubling_back(method, convergence, relTOL):
"""Long steps through the limit point at 10 times the limit load
The previous converged point is at the distance of the arc length from
the last one, and the iterations of the second step may converge to it,
which would trace the path backwards, the load factor decreasing below
zero without ever reaching 1.0. Such a step points against the tangent
predictor, and it is repeated with a shorter arc length.
"""
p = von_mises_truss(frac=10., ks=1.)
p.analysis(NL_method=method, initialInc=0.2, convergence=convergence,
relTOL=relTOL)
p.an.static(NLgeom=True, silent=True)
check_equilibrium(p)
inc = np.asarray(p.an.increments)
cs = np.asarray(p.an.cs)
assert inc[-1] == 1.0
assert np.all(np.diff(np.r_[0., cs[:, 0]]) > 0)
np.testing.assert_allclose(cs[-1],
p.solve_reference(c0=cs[-1], tol=1e-12),
rtol=1e-6)
The snap-through of a low sinusoidal arch of BeamC elements, hinged
at both ends under a sinusoidal load, against the limit loads of Fung and
Kaplan (1952), at the default tolerances: the null load between the limit
points is passed thanks to the round-off floor of the convergence criteria,
see Non-linear static analyses:
def arch(lam1, num_elements=40, L=100., t=1., b=1., E=1.e7, nu=0.3,
R_ref=1.):
r"""Pinned sinusoidal arch of BeamC elements under a sinusoidal load
Returns the :class:`Model`, the mask of the free degrees-of-freedom, the
function that gives `B_1` from the displacements and the load `R` of the
reference load vector, ``R_ref``.
"""
A = b*t
I = b*t**3/12
rA = np.sqrt(A/I)
a1L = 2*lam1/rA
q0 = R_ref*2*np.pi**4*E*I/(L**4*rA)
prop = BeamProp()
prop.A = A
prop.E = E
prop.G = E/2/(1 + nu)
prop.Iyy = I # bending in the xz plane
prop.Izz = t*b**3/12
prop.J = b*t**3/3
num_nodes = num_elements + 1
xs = np.linspace(0, L, num_nodes)
x = np.vstack((xs, 0*xs, a1L*np.sin(np.pi*xs/L))).T.flatten()
data = BeamCData()
probe = BeamCProbe()
elements = []
for k in range(num_elements):
elem = BeamC(probe)
elem.n1, elem.n2 = k, k + 1
elem.c1, elem.c2 = DOF*k, DOF*(k + 1)
elem.init_k_KC0 = k*data.KC0_SPARSE_SIZE
elem.init_k_KCNL = k*data.KCNL_SPARSE_SIZE
elem.init_k_KG = k*data.KG_SPARSE_SIZE
elem.update_rotation_matrix(0., 1., 0., x)
elem.update_probe_xe(x)
elements.append(elem)
N = DOF*num_nodes
bk = np.zeros(N, dtype=bool)
# NOTE planar problem: v, rx and rz restrained at every node
bk[1::DOF] = True
bk[3::DOF] = True
bk[5::DOF] = True
# NOTE hinged ends, u and w restrained
for n in (0, num_nodes - 1):
bk[DOF*n + 0] = True
bk[DOF*n + 2] = True
fext = np.zeros(N)
weights = np.full(num_nodes, L/num_elements)
weights[[0, -1]] /= 2
fext[2::DOF] = -q0*np.sin(np.pi*xs/L)*weights
model = Model(elements, data, prop, x, N, fext)
mid = DOF*(num_elements//2) + 2
def B1(u):
return (a1L + u[mid])/2*rA
return model, ~bk, B1
@pytest.mark.parametrize('convergence_criterion', CONVERGENCE)
@pytest.mark.parametrize('method', METHODS)
def test_fung_kaplan_beam_arch(method, convergence_criterion, monkeypatch):
r"""Symmetric path through both limit points at the default tolerances
The limit loads of 40 elements differ from those of Fung and Kaplan by
-0.0094 (upper) and +0.040 (lower), measured with arc-length increments
of 0.025; with 80 and 160 elements these differences decrease to
-0.0053 and +0.0087, and to -0.0042 and +0.0007. The converged point of
extreme load is within about `3 B_{1,lim} h^2 + h^3 = 0.006` of the
limit load of the model, with `B_{1,lim} = 1.13` and `h = 0.042` the
largest increment of `B_1` of a step, as in Eq. (23). Hence the
tolerance ``0.01*R_cr = 0.051`` of both limit loads.
The round-off floor of the convergence criteria must accept at least one
iteration, see the module documentation.
"""
lam1 = 2.2
R_cr, R_low = limit_loads(lam1)
B1_lim = np.sqrt((lam1**2 - 1)/3)
R_ref = 1.25*R_cr
model, bu, B1 = arch(lam1, R_ref=R_ref)
roundoff = []
def recording_check_convergence(*args, **kwargs):
out = convergence.check_convergence(*args, **kwargs)
if kwargs.get('full_output'):
roundoff.append(out[3])
return out
monkeypatch.setattr(newton_raphson, 'check_convergence',
recording_check_convergence)
monkeypatch.setattr(arc_length, 'check_convergence',
recording_check_convergence)
an = model.analysis(bu, NL_method=method, initialInc=0.05, maxInc=0.05,
convergence=convergence_criterion, maxArcLength=1e6)
assert an.relTOL == 1e-9 and an.absTOL is None
increments, cs = an.static(NLgeom=True, silent=True)
assert increments[-1] == 1.0
assert any(roundoff)
R = R_ref*np.asarray(increments)
B1s = np.array([B1(c) for c in cs])
# NOTE continuous path, without doubling back, to the inverted arch
assert np.all(np.diff(np.r_[lam1, B1s]) < 0)
assert B1s[-1] < -2.
kmin = np.argmin(R)
kmax = np.argmax(R[:kmin])
assert B1s[kmax] > 0 > B1s[kmin]
assert abs(B1s[kmax] - B1_lim) < 0.05
assert abs(B1s[kmin] + B1_lim) < 0.05
assert abs(R[kmax] - R_cr) <= 0.01*R_cr
assert abs(R[kmin] - R_low) <= 0.01*R_cr
# NOTE two null loads, between the limit points and after the lower one
assert np.count_nonzero(R[:-1]*R[1:] < 0) == 2