Prescribed displacements imposed exactly

A prescribed displacement can be imposed exactly as a linear constraint on the degrees of freedom, \(\{c\} = [T]\{c_r\} + \{c_{p0}\} + \lambda \{c_{p1}\}\), where \([T]\) is a basis of the null space of the constraints and \(\{c_{p0}\}\), \(\{c_{p1}\}\) are particular solutions for the constant and the incremented prescribed values, as done by the null-space method of panels. The reduced residual is

\[\{R_r(c_r, \lambda)\} = [T]^T \left(\lambda \{F_{ext}\} - \{F_{int}([T]\{c_r\} + \{c_{p0}\} + \lambda \{c_{p1}\})\}\right)\]

and its derivative with respect to the load factor,

\[\frac{\partial R_r}{\partial \lambda} = [T]^T \left(\{F_{ext}\} - [K_T(c)]\{c_{p1}\}\right)\]

depends on the configuration and on the load factor. calc_fext receives both when it declares the keyword arguments c and lbd, see callbacks. Under a prescribed displacement alone, \(F_{ext} = 0\), the reference vector calc_fext(inc=1.) and the reduced internal force vector both vanish at equilibrium, and the solvers use \(\partial R/\partial \lambda\) for the displacement scale of the arc-length methods, the predictor of the first Newton-Raphson increment and the relative convergence criterion.

Three hardening springs in series, with the displacement of the last node prescribed:

class SpringsProblem(object):
    def __init__(self, k=(100., 150., 200.), a=(20., 40., 60.), d0=0.,
                 d1=1., F=0., with_lbd=True):
        self.k = np.asarray(k)
        self.a = np.asarray(a)
        self.cp0 = np.array([0., 0., d0])
        self.cp1 = np.array([0., 0., d1])
        self.Fext = np.array([F, 0., 0.])
        self.with_lbd = with_lbd
        self.an = None
        # (lbd, c_r) of each call to calc_fext with lbd, and (inc, ||R_r||,
        # number of converged increments) of each call to calc_fint
        self.fext_calls = []
        self.fint_calls = []

    def full(self, c_r, lbd):
        return T @ c_r + self.cp0 + lbd*self.cp1

    def elongations(self, c):
        return np.array([c[0], c[1] - c[0], c[2] - c[1]])

    def fint_full(self, c):
        k, a = self.k, self.a
        f = k*self.elongations(c) + a*self.elongations(c)**3
        return np.array([f[0] - f[1], f[1] - f[2], f[2]])

    def kT_full(self, c):
        ke = self.k + 3*self.a*self.elongations(c)**2
        return np.array([[ke[0] + ke[1], -ke[1], 0.],
                         [-ke[1], ke[1] + ke[2], -ke[2]],
                         [0., -ke[2], ke[2]]])

    def residual(self, c_r, lbd):
        return T.T @ (lbd*self.Fext - self.fint_full(self.full(c_r, lbd)))

    def dR_dlbd(self, c_r, lbd):
        c = self.full(c_r, lbd)
        return T.T @ (self.Fext - self.kT_full(c) @ self.cp1)

    def reaction(self, c_r, lbd):
        r"""Force of the support on node 3, `c_{p1}^T r / d_1`, with `r =
        F_{int}(c) - \lambda F_{ext}`"""
        c = self.full(c_r, lbd)
        r = self.fint_full(c) - lbd*self.Fext
        return self.cp1 @ r/self.cp1[2]

    # callables of structsolve.Analysis
    def calc_fext(self, inc=1., silent=True, c=None, lbd=None):
        if c is None:
            return inc*(T.T @ self.Fext)
        self.fext_calls.append((lbd, c.copy()))
        return inc*self.dR_dlbd(c, lbd)

    def calc_fint(self, c, inc=1., silent=True):
        step = len(self.an.increments) if self.an is not None else 0
        self.fint_calls.append((inc, np.linalg.norm(self.residual(c, inc)),
                                step))
        return T.T @ self.fint_full(self.full(c, inc))

    def calc_kC(self, c=None, NLgeom=False, inc=1., silent=True):
        if c is None or not NLgeom:
            return csr_matrix(T.T @ self.kT_full(np.zeros(3)) @ T)
        return csr_matrix(T.T @ self.kT_full(self.full(c, inc)) @ T)

    def calc_kG(self, c=None, NLgeom=False, silent=True):
        return csr_matrix((2, 2))

    def analysis(self, **kwargs):
        if self.with_lbd:
            calc_fext = self.calc_fext
        else:
            calc_fext = lambda inc=1., silent=True: self.calc_fext(inc)
        an = Analysis(calc_fext, self.calc_fint, self.calc_kC, self.calc_kG)
        for key, value in kwargs.items():
            setattr(an, key, value)
        self.an = an
        return an

    def exact(self, lbd):
        r"""Solution without the load on node 1: the three springs carry the
        same force `P`, found from `\sum_i e_i(P) = x_3`"""
        assert self.Fext[0] == 0.
        x3 = self.cp0[2] + lbd*self.cp1[2]

        def e(P, i):
            return brentq(lambda ei: self.k[i]*ei + self.a[i]*ei**3 - P,
                          -1e4, 1e4, xtol=1e-15, rtol=1e-15)

        def gap(P):
            return sum(e(P, i) for i in range(3)) - x3

        P = 0. if x3 == 0 else brentq(gap, -1e4, 1e4, xtol=1e-14,
                                      rtol=1e-15)
        e1, e2 = e(P, 0), e(P, 1)
        return np.array([e1, e1 + e2]), P

    def solve_reference(self, lbd, tol=1e-14):
        """Dense Newton-Raphson with continuation on the load factor"""
        c_r = np.zeros(2)
        for l in np.linspace(0, lbd, 21)[1:]:
            for _ in range(100):
                R = self.residual(c_r, l)
                if np.linalg.norm(R) <= tol*max(1., np.linalg.norm(
                        self.dR_dlbd(c_r, l))):
                    break
                kT = T.T @ self.kT_full(self.full(c_r, l)) @ T
                c_r = c_r + np.linalg.solve(kT, R)
        return c_r

@pytest.mark.parametrize('method', ['NR', 'arc_length_riks',
                                    'arc_length_crisfield'])
def test_displacement_control_relative_criterion(method):
    r"""The prescribed displacement alone, with the relative convergence
    criterion only: every converged increment is the exact solution, and
    ``calc_fext`` receives the load factor of the configuration"""
    pb = SpringsProblem(d1=0.8)
    an = pb.analysis(NL_method=method, initialInc=0.2, relTOL=1e-10,
                     absTOL=None)
    an.static(NLgeom=True, silent=True)
    assert np.isclose(an.increments[-1], 1.)
    assert len(an.increments) > 2
    for lbd, c_r in zip(an.increments, an.cs):
        ref, P = pb.exact(lbd)
        assert np.allclose(c_r, ref, rtol=1e-8, atol=1e-12), (lbd, c_r, ref)
        assert np.isclose(pb.reaction(c_r, lbd), P, rtol=1e-8)
        assert (np.linalg.norm(pb.residual(c_r, lbd))
                <= 1e-9*np.linalg.norm(lbd*pb.dR_dlbd(c_r, lbd)))
    # the arc-length solvers evaluate dR/dlbd at the configurations and load
    # factors of their iterations
    if method != 'NR':
        assert len(pb.fext_calls) > len(an.increments)
        lbds = [lbd for lbd, _ in pb.fext_calls]
        assert lbds[0] == 0.
        assert all(lbd is not None for lbd in lbds)
        assert len(set(lbds)) > 2

A constant prescribed displacement together with an incremented one and a load:

@pytest.mark.parametrize('method', ['NR', 'arc_length_riks',
                                    'arc_length_crisfield'])
def test_constant_and_incremented_prescriptions(method):
    r"""A constant prescribed displacement `d_0`, held during the analysis,
    with a load `\lambda F` on node 1 and an incremented displacement
    `\lambda d_1`"""
    pb = SpringsProblem(d0=-0.3, d1=0.6, F=40.)
    an = pb.analysis(NL_method=method, initialInc=0.2, relTOL=1e-11)
    an.static(NLgeom=True, silent=True)
    assert np.isclose(an.increments[-1], 1.)
    for lbd, c_r in zip(an.increments, an.cs):
        ref = pb.solve_reference(lbd)
        assert np.allclose(c_r, ref, rtol=1e-8, atol=1e-12), (lbd, c_r, ref)
    # at lbd = 0 the constant displacement alone is not in equilibrium with
    # c_r = 0: the first increment starts from a non-null residual
    assert np.linalg.norm(pb.residual(np.zeros(2), 0.)) > 1.