Linear static analysis

A rectangular isotropic plate, with \(w\) constrained at the nodes of the four edges, under a uniform pressure. The stiffness matrix is assembled with update_KC0() and the consistent load vector of the pressure is integrated element by element with the rows of BFSCPlate2D.update_Sw(); \(u\) and \(v\) are constrained at the central node only, which removes the in-plane rigid-body motions. The code is extracted from one of the bfscplate2d unit tests:

def test_static_pressure(plot=False):
    # number of nodes
    nx = 5
    ny = 7
    points, weights = get_points_weights(nint=4)

    # geometry
    a = 3
    b = 7

    # material properties
    E = 200e9
    nu = 0.3
    h = 0.005
    lam = isotropic_plate(thickness=h, E=E, nu=nu)

    xtmp = np.linspace(0, a, nx)
    ytmp = np.linspace(0, b, ny)
    xmesh, ymesh = np.meshgrid(xtmp, ytmp)

    # getting nodes
    ncoords = np.vstack((xmesh.T.flatten(), ymesh.T.flatten())).T
    x = ncoords[:, 0]
    y = ncoords[:, 1]
    nids = 1 + np.arange(ncoords.shape[0])
    nid_pos = dict(zip(nids, np.arange(len(nids))))

    nids_mesh = nids.reshape(nx, ny)

    n1s = nids_mesh[:-1, :-1].flatten()
    n2s = nids_mesh[1:, :-1].flatten()
    n3s = nids_mesh[1:, 1:].flatten()
    n4s = nids_mesh[:-1, 1:].flatten()

    num_elements = len(n1s)
    print('num_elements', num_elements)

    N = DOF*nx*ny
    Kr = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=INT)
    Kc = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=INT)
    Kv = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=DOUBLE)
    init_k_KC0 = 0

    plates = []
    for n1, n2, n3, n4 in zip(n1s, n2s, n3s, n4s):
        plate = BFSCPlate2D()
        plate.c1 = DOF*nid_pos[n1]
        plate.c2 = DOF*nid_pos[n2]
        plate.c3 = DOF*nid_pos[n3]
        plate.c4 = DOF*nid_pos[n4]
        plate.ABD = lam.ABD
        plate.lex = a/(nx - 1)
        plate.ley = b/(ny - 1)
        plate.init_k_KC0 = init_k_KC0
        update_KC0(plate, points, weights, Kr, Kc, Kv)
        init_k_KC0 += KC0_SPARSE_SIZE
        plates.append(plate)

    KC0 = coo_matrix((Kv, (Kr, Kc)), shape=(N, N)).tocsc()

    # applying boundary conditions
    bk = np.zeros(KC0.shape[0], dtype=bool)

    # simply supported
    check = isclose(x, 0) | isclose(x, a) | isclose(y, 0) | isclose(y, b)
    bk[6::DOF] = check
    # point supports (middle)
    check = isclose(x, a/2) & isclose(y, b/2)
    bk[0::DOF] = check
    bk[3::DOF] = check

    bu = ~bk

    # external force vector applying consistent pressure
    f = np.zeros(KC0.shape[0], dtype=float)
    P = 10 # Pa
    nint = 4
    points, weights = get_points_weights(nint)
    for plate in plates:
        lex = plate.lex
        ley = plate.ley
        indices = []
        c1 = plate.c1
        c2 = plate.c2
        c3 = plate.c3
        c4 = plate.c4
        cs = [c1, c2, c3, c4]
        for ci in cs:
            for i in range(DOF):
                indices.append(ci + i)
        fe = np.zeros(4*DOF, dtype=float)
        for i in range(nint):
            xi = points[i]
            weight_xi = weights[i]
            for j in range(nint):
                eta = points[j]
                weight_eta = weights[j]
                weight = weight_xi*weight_eta
                plate.update_Sw(xi, eta)
                Sw = np.asarray(plate.Sw)
                fe += (lex*ley)/4*weight*Sw*P
        f[indices] += fe

    Kuu = KC0[bu, :][:, bu]
    fu = f[bu]

    # solving
    uu = spsolve(Kuu, fu)
    u = np.zeros(KC0.shape[0], dtype=float)
    u[bu] = uu

    w = u[6::DOF].reshape(nx, ny).T
    print('wmax', w.max())
    print('wmin', w.min())
    if plot:
        import matplotlib
        matplotlib.use('TkAgg')
        import matplotlib.pyplot as plt
        plt.gca().set_aspect('equal')
        levels = np.linspace(w.min(), w.max(), 300)
        plt.contourf(xmesh, ymesh, w, levels=levels)
        plt.colorbar()
        plt.show()
    assert np.isclose(w.max(), 3.9407e-3, rtol=1e-3)