Matrices of the third-order shear deformation theory

The TSDT uses the in-plane matrices A, B, D, E, F, H and the uncorrected transverse shear matrices Abar_ts, Dts, Fts, see Laminate constitutive matrices per plate theory. The test below assembles the 13 x 13 constitutive matrix of the generalized strains \(\{\varepsilon^{(0)}, \varepsilon^{(1)}, \varepsilon^{(3)}, \gamma^{(0)}, \gamma^{(2)}\}\) and checks that it gives the strain energy per unit area obtained by integrating the ply strain energy density through the thickness:

def test_example_tsdt_matrices():
    # unsymmetric laminate with an offset
    laminaprop = (142e9, 8.7e9, 0.28, 5.1e9, 5.1e9, 3.2e9)
    lam = laminated_plate([0, 45, 90, -30], plyts=[0.2e-3, 0.5e-3, 0.3e-3,
        0.4e-3], laminaprop=laminaprop, offset=0.37e-3)
    # 13 x 13 constitutive matrix of the TSDT, the generalized strains being
    # {eps0, eps1, eps3, gamma0, gamma2}, with eps = eps0 + z eps1 + z^3 eps3
    # and gamma = {gamma_yz, gamma_xz} = gamma0 + z^2 gamma2
    Z3 = np.zeros((3, 3))
    Z32 = np.zeros((3, 2))
    K = np.block([[lam.A, lam.B, lam.E, Z32, Z32],
                  [lam.B, lam.D, lam.F, Z32, Z32],
                  [lam.E, lam.F, lam.H, Z32, Z32],
                  [Z32.T, Z32.T, Z32.T, lam.Abar_ts, lam.Dts],
                  [Z32.T, Z32.T, Z32.T, lam.Dts, lam.Fts]])
    assert np.allclose(K, K.T, rtol=0, atol=0)
    # strain energy per unit area integrated ply by ply, against e^T K e / 2
    rng = np.random.default_rng(1)
    h = lam.h
    scale = np.array(3*[1.] + 3*[1/h] + 3*[1/h**3] + 2*[1.] + 2*[1/h**2])
    xi, w = np.polynomial.legendre.leggauss(6)
    zbot = -h/2 + lam.offset
    for _ in range(5):
        e = 1e-3*rng.uniform(-1, 1, 13)*scale
        U = 0.
        zk = zbot
        for ply in lam.plies:
            C = np.asarray(ply.get_constitutive_matrix())
            for x, wi in zip(xi, w):
                z = zk + ply.h*(x + 1)/2
                eps = e[0:3] + z*e[3:6] + z**3*e[6:9]
                gamma = e[9:11] + z**2*e[11:13]
                U += wi*ply.h/2*0.5*(eps @ C[:3, :3] @ eps
                                     + gamma @ C[3:, 3:] @ gamma)
            zk += ply.h
        assert np.isclose(0.5*e @ K @ e, U, rtol=1e-12)