!python -m pip install buckling composites structsolve plotly > tmp.txt

import numpy as np
from scipy.special import roots_legendre
from composites import isotropic_plate

from buckling.legendre import vecf, vecfxi, vecfxixi

# approximation order
m1 = m2 = 20
N = 3*m1*m2
print('m1 = %d' % m1)
print('m2 = %d' % m2)
print('N = %d' % N)
i = np.arange(m1)
j = np.arange(m2)
pts1, weights1 = roots_legendre(2*m1 - 1)
pts2, weights2 = roots_legendre(2*m2 - 1)

# Material properties
E = 200.e9
nu = 0.3
G = E/(2*(1 + nu))

# Boundary conditions w (simply supported)
wxit1 = 0
wxir1 = 1
wxit2 = 0
wxir2 = 1
wetat1 = 0
wetar1 = 1
wetat2 = 0
wetar2 = 1

# Boundary conditions phix (simply supported: constrained at y=0, b)
phix_xit1 = 1
phix_xir1 = 1
phix_xit2 = 1
phix_xir2 = 1
phix_etat1 = 0
phix_etar1 = 1  # Derivative producing twist, must be free
phix_etat2 = 0
phix_etar2 = 1  # Derivative producing twist, must be free

# Boundary conditions phiy (simply supported: constrained at x=0, a)
phiy_xit1 = 0
phiy_xir1 = 1  # Derivative producing twist, must be free
phiy_xit2 = 0
phiy_xir2 = 1  # Derivative producing twist, must be free
phiy_etat1 = 1
phiy_etar1 = 1
phiy_etat2 = 1
phiy_etar2 = 1

# Geometric properties
a = 3
b = 7
h = 0.005

# Calculating ABDEFH matrices
prop = isotropic_plate(thickness=h, E=E, nu=nu)

Swx = np.zeros(N)
Swxx = np.zeros(N)
Swyy = np.zeros(N)
Swxy = np.zeros(N)
Swy = np.zeros(N)
Sphix = np.zeros(N)
Sphiy = np.zeros(N)
Sphixx = np.zeros(N)
Sphixy = np.zeros(N)
Sphiyx = np.zeros(N)
Sphiyy = np.zeros(N)
    
K = np.zeros((N, N))

detJ = a * b / 4

# stiffness matrix
# numerical integration in 2D using Legendre-Gauss quadrature
for xi, wxi in zip(pts1, weights1):
    # 1D functions for w
    wP_xi   = vecf(m1, xi, wxit1, wxir1, wxit2, wxir2)
    wPx_xi  = vecfxi(m1, xi, wxit1, wxir1, wxit2, wxir2)
    wPxx_xi = vecfxixi(m1, xi, wxit1, wxir1, wxit2, wxir2)
    
    # 1D functions for phix
    P_phix_xi  = vecf(m1, xi, phix_xit1, phix_xir1, phix_xit2, phix_xir2)
    Px_phix_xi = vecfxi(m1, xi, phix_xit1, phix_xir1, phix_xit2, phix_xir2)
    
    # 1D functions for phiy
    P_phiy_xi  = vecf(m1, xi, phiy_xit1, phiy_xir1, phiy_xit2, phiy_xir2)
    Px_phiy_xi = vecfxi(m1, xi, phiy_xit1, phiy_xir1, phiy_xit2, phiy_xir2)
    
    for eta, weta in zip(pts2, weights2):
        # 1D functions for w
        wP_eta   = vecf(m2, eta, wetat1, wetar1, wetat2, wetar2)
        wPx_eta  = vecfxi(m2, eta, wetat1, wetar1, wetat2, wetar2)
        wPxx_eta = vecfxixi(m2, eta, wetat1, wetar1, wetat2, wetar2)

        # 1D functions for phix
        P_phix_eta  = vecf(m2, eta, phix_etat1, phix_etar1, phix_etat2, phix_etar2)
        Px_phix_eta = vecfxi(m2, eta, phix_etat1, phix_etar1, phix_etat2, phix_etar2)

        # 1D functions for phiy
        P_phiy_eta  = vecf(m2, eta, phiy_etat1, phiy_etar1, phiy_etat2, phiy_etar2)
        Px_phiy_eta = vecfxi(m2, eta, phiy_etat1, phiy_etar1, phiy_etat2, phiy_etar2)

        # --- Broadcasting applied to w DOFs ---
        Swx[:m1*m2]  = (wPx_xi[:, None] * wP_eta[None, :] * (2/a)).flatten()
        Swxx[:m1*m2] = (wPxx_xi[:, None] * wP_eta[None, :] * (2/a)**2).flatten()
        Swy[:m1*m2]  = (wP_xi[:, None] * wPx_eta[None, :] * (2/b)).flatten()
        Swyy[:m1*m2] = (wP_xi[:, None] * wPxx_eta[None, :] * (2/b)**2).flatten()
        Swxy[:m1*m2] = (wPx_xi[:, None] * wPx_eta[None, :] * (2/a) * (2/b)).flatten()
                
        # --- Broadcasting applied to phix and phiy DOFs strictly mapped ---
        Sphix[m1*m2:2*m1*m2]  = (P_phix_xi[:, None] * P_phix_eta[None, :]).flatten()
        Sphiy[2*m1*m2:]       = (P_phiy_xi[:, None] * P_phiy_eta[None, :]).flatten()
        
        Sphixx[m1*m2:2*m1*m2] = (Px_phix_xi[:, None] * P_phix_eta[None, :] * (2/a)).flatten()
        Sphiyx[2*m1*m2:]      = (Px_phiy_xi[:, None] * P_phiy_eta[None, :] * (2/a)).flatten()
        
        Sphixy[m1*m2:2*m1*m2] = (P_phix_xi[:, None] * Px_phix_eta[None, :] * (2/b)).flatten()
        Sphiyy[2*m1*m2:]      = (P_phiy_xi[:, None] * Px_phiy_eta[None, :] * (2/b)).flatten()
        
        e1xx = Sphixx
        e1yy = Sphiyy
        e1xy = Sphixy + Sphiyx

        e3xx = (-4/(3*h**2))*(Sphixx + Swxx)
        e3yy = (-4/(3*h**2))*(Sphiyy + Swyy)
        e3xy = (-4/(3*h**2))*(Sphixy + Sphiyx + 2*Swxy)
        
        g0yz = Sphiy + Swy
        g0xz = Sphix + Swx

        g2yz = (-4/(h**2))*(Sphiy + Swy)
        g2xz = (-4/(h**2))*(Sphix + Swx)

        # NOTE symmetric plate, ignoring the effect of the Bij and Eij constitutive terms
        Mxx = prop.D11*e1xx + prop.D12*e1yy + prop.D16*e1xy + prop.F11*e3xx + prop.F12*e3yy + prop.F16*e3xy
        Myy = prop.D12*e1xx + prop.D22*e1yy + prop.D26*e1xy + prop.F12*e3xx + prop.F22*e3yy + prop.F26*e3xy
        Mxy = prop.D16*e1xx + prop.D26*e1yy + prop.D66*e1xy + prop.F16*e3xx + prop.F16*e3yy + prop.F66*e3xy

        Pxx = prop.F11*e1xx + prop.F12*e1yy + prop.F16*e1xy + prop.H11*e3xx + prop.H12*e3yy + prop.H16*e3xy
        Pyy = prop.F12*e1xx + prop.F22*e1yy + prop.F26*e1xy + prop.H12*e3xx + prop.H22*e3yy + prop.H26*e3xy
        Pxy = prop.F16*e1xx + prop.F26*e1yy + prop.F66*e1xy + prop.H16*e3xx + prop.H16*e3yy + prop.H66*e3xy
        
        Qy = prop.A44*g0yz + prop.A45*g0xz + prop.D44*g2yz + prop.D45*g2xz
        Qx = prop.A45*g0yz + prop.A55*g0xz + prop.D45*g2yz + prop.D55*g2xz

        Ry = prop.D44*g0yz + prop.D45*g0xz + prop.F44*g2yz + prop.F45*g2xz
        Rx = prop.D45*g0yz + prop.D55*g0xz + prop.F45*g2yz + prop.F55*g2xz
        
        weight = wxi * weta

        # stiffness matrix
        K += np.outer(detJ*weight*Mxx, e1xx)
        K += np.outer(detJ*weight*Myy, e1yy)
        K += np.outer(detJ*weight*Mxy, e1xy)
        K += np.outer(detJ*weight*Pxx, e3xx)
        K += np.outer(detJ*weight*Pyy, e3yy)
        K += np.outer(detJ*weight*Pxy, e3xy)
        K += np.outer(detJ*weight*Qy, g0yz)
        K += np.outer(detJ*weight*Qx, g0xz)
        K += np.outer(detJ*weight*Ry, g2yz)
        K += np.outer(detJ*weight*Rx, g2xz)
       

xi = 0 # x = a/2
eta = 0 # y = b/2
P_xi = vecf(m1, xi, wxit1, wxir1, wxit2, wxir2)
P_eta = vecf(m2, eta, wetat1, wetar1, wetat2, wetar2)
Sw = np.zeros(N)
Sw[:m1*m2] = (P_xi[:, None] * P_eta[None, :]).flatten()

Pforce = 1.

Fext = Pforce*Sw

from scipy.sparse import csc_matrix
from structsolve import solve

u = solve(csc_matrix(K), Fext)

nx, ny = 40, 20
X, Y = np.meshgrid(
    np.linspace(0, a, nx),
    np.linspace(0, b, ny),
    indexing='ij'
)

W = np.zeros_like(X)

Sw = np.zeros(N)

for ij, x in np.ndenumerate(X):
    xi = 2*x/a - 1
    eta = 2*Y[ij]/b - 1

    wP_xi = vecf(m1, xi, wxit1, wxir1, wxit2, wxir2)
    wP_eta = vecf(m2, eta, wetat1, wetar1, wetat2, wetar2)

    Sw[:m1*m2] = (wP_xi[:, None] * wP_eta[None, :]).ravel()

    W[ij] = Sw @ u

print('W.min(), W.max()', W.min(), W.max())
print('W.max() FE reference, simply supported', 6.594931610258557e-05) 

import plotly.graph_objects as go

fig = go.Figure(data=[go.Surface(x=X, y=Y, z=W)])

fig.update_layout(
    title='Plate Deflection (TSDT)',
    scene=dict(
        xaxis_title='X (m)',
        yaxis_title='Y (m)',
        zaxis_title='W (m)',
        aspectmode='data',
        camera=dict(
            eye=dict(x=0, y=0, z=100*a/2),
            center=dict(x=a/2, y=b/2, z=0),
            up=dict(x=0, y=1, z=0)
        )
    ),
    width=800,
    height=600
)

fig.show()