Nonlinear static analysis with structsolve

The callables of structsolve.Analysis, of structsolve, for the cylinder of Nonlinear static analysis, see Solving with structsolve. They loop over the elements with update_fint(), update_KCNL() and update_KG(), assemble, and work on the unknown degrees of freedom bu, the partition of the free degrees of freedom:

class Callables(object):
    """The callables of structsolve.Analysis for the ``elements`` of a model
    under the external forces ``fext``.

    They loop over the elements with the update functions and assemble, and
    work on the unknown degrees of freedom ``bu`` only, the partition of the
    free dofs of structsolve: the solution vector of the solvers is ``c =
    u[bu]``, and the force vectors and matrices are restricted to ``bu``.
    The prescribed degrees of freedom are zero.
    """
    def __init__(self, elements, points, weights, KC0, bu, fext):
        self.elements = elements
        self.points = points
        self.weights = weights
        self.N = KC0.shape[0]
        self.KC0 = KC0
        self.KC0uu = KC0[bu, :][:, bu]
        self.bu = bu
        self.fext = fext
        num_elements = len(elements)
        self.KCNLr = np.zeros(KCNL_SPARSE_SIZE*num_elements, dtype=INT)
        self.KCNLc = np.zeros(KCNL_SPARSE_SIZE*num_elements, dtype=INT)
        self.KCNLv = np.zeros(KCNL_SPARSE_SIZE*num_elements, dtype=DOUBLE)
        self.KGr = np.zeros(KG_SPARSE_SIZE*num_elements, dtype=INT)
        self.KGc = np.zeros(KG_SPARSE_SIZE*num_elements, dtype=INT)
        self.KGv = np.zeros(KG_SPARSE_SIZE*num_elements, dtype=DOUBLE)

    def expand(self, c):
        """The displacements of the model from those of the unknown dofs."""
        u = np.zeros(self.N, dtype=DOUBLE)
        u[self.bu] = c
        return u

    def calc_fext(self, inc=1., silent=False):
        return inc*self.fext[self.bu]

    def calc_fint(self, c, silent=False):
        fint = np.zeros(self.N, dtype=DOUBLE)
        u = self.expand(c)
        for shell in self.elements:
            update_fint(u, shell, self.points, self.weights, fint)
        return fint[self.bu]

    def calc_kC(self, c=None, NLgeom=False, silent=False):
        if not NLgeom:
            return self.KC0uu
        # the generated code accumulates with "+=" over the integration
        # points, so the value arrays must be zeroed on every call. The row
        # and column arrays are rewritten identically each time, and could
        # be filled once for a fixed mesh
        self.KCNLv *= 0
        u = self.expand(c)
        for shell in self.elements:
            update_KCNL(u, shell, self.points, self.weights, self.KCNLr,
                        self.KCNLc, self.KCNLv)
        KCNL = coo_matrix((self.KCNLv, (self.KCNLr, self.KCNLc)),
                          shape=(self.N, self.N)).tocsc()
        # KCNL also carries the geometric term of the nonlinear membrane
        # stress, so that KG stays homogeneous of degree one in u and keeps
        # its linear buckling meaning. KC0 + KCNL + KG is d(fint)/du
        return (self.KC0 + KCNL)[self.bu, :][:, self.bu]

    def calc_kG(self, c=None, NLgeom=False, silent=False):
        self.KGv *= 0
        u = self.expand(c)
        for shell in self.elements:
            update_KG(u, shell, self.points, self.weights, self.KGr,
                      self.KGc, self.KGv)
        KG = coo_matrix((self.KGv, (self.KGr, self.KGc)),
                        shape=(self.N, self.N)).tocsc()
        return KG[self.bu, :][:, self.bu]

The analysis, at the end of the test of Nonlinear static analysis, applies the load in one increment, Analysis.initialInc = 1., with at most seven corrections, Analysis.maxNumIter = 7, and checks that the increment was not cut:

callables = Callables(elements, points, weights, KC0, bu, fext)
an = Analysis(calc_fext=callables.calc_fext,
              calc_fint=callables.calc_fint,
              calc_kC=callables.calc_kC,
              calc_kG=callables.calc_kG)
an.initialInc = 1.
# at most 7 corrections, otherwise structsolve cuts the load step
an.maxNumIter = 7
increments, cs = an.static(NLgeom=True)
if increments != [1.]:
    raise RuntimeError('Not converged!')

The solution is u = callables.expand(cs[-1]). With the default convergence criterion of Crisfield, \(\zeta \le 10^{-9}\), the iterations give \(\zeta\) = 2.5e-2, 2.2e-4, 2.0e-9 and then the round-off, about 1e-15, the quadratic convergence of a consistent tangent stiffness matrix.

Follower pressure

The callables of a follower pressure, FollowerCallables, which receive the load factor inc, are given in Solving with structsolve. A cylinder clamped at one end and free at the other, with an unsymmetric load stiffness matrix, solved with these callables for both kinematics, checking the equilibrium with the whole follower pressure and the quadratic convergence of the Newton-Raphson iterations:

@pytest.mark.parametrize('kinematics', sorted(KINEMATICS))
def test_clamped_free_cylinder(kinematics):
    r"""Cylinder of radius 1 m, length 1 m and thickness 10 mm, clamped at
    `x = 0` and free at `x = L`, under the external follower pressure `p =
    -2 (1 + \cos(3 y/R)/2)` bar, in two equal load increments"""
    R, L, h = 1., 1., 0.01
    prop = isotropic_plate(thickness=h, E=E, nu=NU)
    mod, elements, points, weights, x, y, N = cylinder_model(kinematics, L,
            R, prop, nx=7, ny=24)
    DOF = mod.DOF
    bk = np.zeros(N, dtype=bool)
    for i in range(DOF):
        bk[i::DOF] |= isclose(x, 0)
    bu = ~bk
    nint = len(points)
    # p[i, j] at xi = points[i] and eta = points[j], from the circumferential
    # coordinate of node 1 of each element
    pressures = []
    for shell in elements:
        yj = y[shell.n1] + shell.ley*(points + 1)/2
        p = -2.e5*(1 + 0.5*np.cos(3*yj/R))
        pressures.append(np.ascontiguousarray(np.ones((nint, 1))*p[None, :]))

    callables = FollowerCallables(mod, elements, points, weights, N, bu,
                                  pressures)
    an = Analysis(calc_fext=callables.calc_fext,
                  calc_fint=callables.calc_fint,
                  calc_kC=callables.calc_kC,
                  calc_kG=callables.calc_kG)
    an.initialInc = 0.5
    an.maxInc = 0.5
    increments, cs = an.static(NLgeom=True)
    # every increment converged, none was cut
    assert np.allclose(increments, [0.5, 1.]), increments
    u = callables.expand(cs[-1])
    print('max |w|/h', np.abs(u[6::DOF]).max()/h)

    # equilibrium with the whole follower pressure, F_p(u) of the pressure p
    fint = np.zeros(N)
    fp = np.zeros(N)
    for shell, p in zip(elements, pressures):
        mod.update_fint(u, shell, points, weights, fint)
        mod.update_fext_follower(u, shell, points, weights, p, fp)
    residual = np.linalg.norm((fint - fp)[bu])/np.linalg.norm(fp[bu])
    print('relative residual', residual)
    assert residual < 1.e-8
    # the follower part of the load is not negligible
    assert np.linalg.norm((fp - callables.fp0)[bu]) > 1.e-3*np.linalg.norm(
            callables.fp0[bu])

    # the free end makes the load stiffness matrix unsymmetric
    KC = callables.calc_kC(cs[-1], NLgeom=True, inc=1.)
    KCnl = callables.calc_kC(cs[-1], NLgeom=True, inc=0.)
    KCf = (KC - KCnl).toarray()
    asym = np.abs(KCf - KCf.T).max()/np.abs(KCf).max()
    print('asymmetry of the load stiffness matrix', asym)
    assert asym > 1.e-3

    # quadratic convergence in each increment: Crisfield's zeta, the
    # criterion of structsolve, is squared by the last iteration, down to
    # the round-off
    w = crisfield_weights(callables.KC0[bu, :][:, bu])
    f = callables.fp0[bu]
    for lbd in increments:
        history = [crisfield_zeta(inc*f - fi, inc*f, fi, w)
                   for inc, fi in callables.iterations if inc == lbd]
        print('load factor', lbd, 'zeta', history)
        assert history[-1] < 1.e-9
        assert len(history) <= 6
        prev, curr = history[-2], history[-1]
        assert curr < max(10.*prev**2, 1.e-13), history