Non-linear static analysis with the arc-length methods

The arc-length methods, described in arc_length._solver_arc_length(), trace equilibrium paths that include limit points, where the Newton-Raphson method with load control fails. The shallow von Mises truss below, loaded through a linear spring, softens until a limit point, snaps through to an inverted configuration, and stiffens again:

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)

The callables are defined as in the example Non-linear static analysis with the Newton-Raphson method. With a load of 1.5 times the limit load, the arc-length methods, selected with NL_method='arc_length_riks' or NL_method='arc_length_crisfield', trace the path through the limit point, where the load factor decreases and becomes negative, until the load factor reaches exactly 1.0 on the inverted configuration. The code is extracted from one of the structsolve unit tests:

@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)