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)

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)