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

import numpy as np
from scipy.special import roots_legendre

from buckling.legendre import vecf, vecfxi


# approximation order
m1 = m2 = 10
m3 = 5

# number of degrees of freedom
N = 3*m1*m2*m3

print('m1 = %d' % m1)
print('m2 = %d' % m2)
print('m3 = %d' % m3)
print('N = %d' % N)

i = np.arange(m1)
j = np.arange(m2)
k = np.arange(m3)
pts1, weights1 = roots_legendre(2*m1 - 1)
pts2, weights2 = roots_legendre(2*m2 - 1)
pts3, weights3 = roots_legendre(2*m3 - 1)

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

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

# Boundary conditions u v w
xit1 = 0
xir1 = 1
xit2 = 0
xir2 = 1
etat1 = 0
etar1 = 1
etat2 = 0
etar2 = 1
zetat1 = 1
zetar1 = 1
zetat2 = 1
zetar2 = 1

Su = np.zeros(N)
Sv = np.zeros(N)
Sw = np.zeros(N)
Sux = np.zeros(N)
Svx = np.zeros(N)
Swx = np.zeros(N)
Suy = np.zeros(N)
Svy = np.zeros(N)
Swy = np.zeros(N)
Suz = np.zeros(N)
Svz = np.zeros(N)
Swz = np.zeros(N)

# memory buffers for stiffness matrix assembly
buff = np.zeros((N, N))
K = np.zeros((N, N))

def addouter(matrix, vec1, vec2):
    np.outer(vec1, vec2, out=buff)
    matrix += buff

# numerical integration in 3D using Gauss quadrature
for xi, wxi in zip(pts1, weights1):
    P_xi = vecf(m1, xi, xit1, xir1, xit2, xir2)
    Px_xi = vecfxi(m1, xi, xit1, xir1, xit2, xir2)
    for eta, weta in zip(pts2, weights2):
        P_eta = vecf(m2, eta, etat1, etar1, etat2, etar2)
        Px_eta = vecfxi(m2, eta, etat1, etar1, etat2, etar2)
        for zeta, wzeta in zip(pts3, weights3):
            P_zeta = vecf(m3, zeta, zetat1, zetar1, zetat2, zetar2)
            Px_zeta = vecfxi(m3, zeta, zetat1, zetar1, zetat2, zetar2)

            weight = wxi*weta*wzeta

            Pi, Pj, Pk = np.meshgrid(P_xi, P_eta, P_zeta, indexing='ij')
            Su[:m1*m2*m3] = (Pi*Pj*Pk).flatten()
            Sv[m1*m2*m3:2*m1*m2*m3] = (Pi*Pj*Pk).flatten()
            Sw[2*m1*m2*m3:] = (Pi*Pj*Pk).flatten()

            Pxi, Pj, Pk = np.meshgrid(Px_xi, P_eta, P_zeta, indexing='ij')
            Sux[:m1*m2*m3] = (Pxi*Pj*Pk*(2/a)).flatten()
            Svx[m1*m2*m3:2*m1*m2*m3] = (Pxi*Pj*Pk*(2/a)).flatten()
            Swx[2*m1*m2*m3:] = (Pxi*Pj*Pk*(2/a)).flatten()

            Pi, Pxj, Pk = np.meshgrid(P_xi, Px_eta, P_zeta, indexing='ij')
            Suy[:m1*m2*m3] = (Pi*Pxj*Pk*(2/b)).flatten()
            Svy[m1*m2*m3:2*m1*m2*m3] = (Pi*Pxj*Pk*(2/b)).flatten()
            Swy[2*m1*m2*m3:] = (Pi*Pxj*Pk*(2/b)).flatten()

            Pi, Pj, Pxk = np.meshgrid(P_xi, P_eta, Px_zeta, indexing='ij')
            Suz[:m1*m2*m3] = (Pi*Pj*Pxk*(2/h)).flatten()
            Svz[m1*m2*m3:2*m1*m2*m3] = (Pi*Pj*Pxk*(2/h)).flatten()
            Swz[2*m1*m2*m3:] = (Pi*Pj*Pxk*(2/h)).flatten()

            NL = 0 # NOTE set to 1 for full nonlinear relations, 0 for linear relations
            exx = Sux + NL*1/2*(Sux**2 + Svx**2 + Swx**2)
            eyy = Svy + NL*1/2*(Suy**2 + Svy**2 + Swy**2)
            ezz = Swz + NL*1/2*(Suz**2 + Svz**2 + Swz**2)
            gxy = Suy + Svx + NL*(Sux*Suy + Svx*Svy + Swx*Swy)
            gxz = Suz + Swx + NL*(Sux*Suz + Svx*Svz + Swx*Swz)
            gyz = Svz + Swy + NL*(Suy*Suz + Svy*Svz + Swy*Swz)

            C = E/((1 - 2*nu)*(1 + nu))
            Sxx = C*((1 - nu)*exx + nu*eyy + nu*ezz)
            Syy = C*(nu*exx + (1 - nu)*eyy + nu*ezz)
            Szz = C*(nu*exx + nu*eyy + (1 - nu)*ezz)
            Txy = G*gxy
            Txz = G*gxz
            Tyz = G*gyz

            detJ = a*b*h/8

            # stiffness matrix
            addouter(K, detJ*weight*Sxx, exx)
            addouter(K, detJ*weight*Syy, eyy)
            addouter(K, detJ*weight*Szz, ezz)
            addouter(K, detJ*weight*Txy, gxy)
            addouter(K, detJ*weight*Txz, gxz)
            addouter(K, detJ*weight*Tyz, gyz)

xi = 0 # x = a/2
eta = 0 # y = b/2
zeta = 1 # z = +h/2
P_xi = vecf(m1, xi, xit1, xir1, xit2, xir2)
P_eta = vecf(m2, eta, etat1, etar1, etat2, etar2)
P_zeta = vecf(m3, zeta, zetat1, zetar1, zetat2, zetar2)
Pi, Pj, Pk = np.meshgrid(P_xi, P_eta, P_zeta, indexing='ij')
Sw[2*m1*m2*m3:] = (Pi*Pj*Pk).flatten()
Pforce = 1.
Fext = Pforce*Sw

from scipy.sparse import csc_matrix
from structsolve import solve

u = solve(csc_matrix(K), Fext)

nx, ny, nz = 40, 20, 3
X, Y, Z = np.meshgrid(
    np.linspace(0, a, nx),
    np.linspace(0, b, ny),
    np.linspace(-h/2, h/2, nz),
    indexing='ij'
)

U = np.zeros_like(X)
V = np.zeros_like(Y)
W = np.zeros_like(Z)

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

    P_xi = vecf(m1, xi, xit1, xir1, xit2, xir2)
    P_eta = vecf(m2, eta, etat1, etar1, etat2, etar2)
    P_zeta = vecf(m3, zeta, zetat1, zetar1, zetat2, zetar2)

    Pi, Pj, Pk = np.meshgrid(P_xi, P_eta, P_zeta, indexing='ij')
    Su[:m1*m2*m3] = (Pi*Pj*Pk).flatten()
    Sv[m1*m2*m3:2*m1*m2*m3] = (Pi*Pj*Pk).flatten()
    Sw[2*m1*m2*m3:] = (Pi*Pj*Pk).flatten()

    U[ijk] = Su @ u
    V[ijk] = Sv @ u
    W[ijk] = Sw @ u

print('W.min(), W.max()', W.min(), W.max())
print('W.max() CLPT clamped', 2.774985496231098e-05)

import plotly.graph_objects as go

scale = 1e4  # deformation scaling factor

X_def = X + scale * U
Y_def = Y + scale * V
Z_def = Z + scale * W

# Pre-compute the 3D hover text array
# Using scientific notation for displacements since they are typically small (e.g., CLPT W.max ~ 2.7e-5)
hover_text = np.empty((nx, ny, nz), dtype=object)
for i in range(nx):
    for j in range(ny):
        for k in range(nz):
            hover_text[i, j, k] = (
                f"U: {U[i,j,k]:.4e}<br>"
                f"V: {V[i,j,k]:.4e}<br>"
                f"W: {W[i,j,k]:.4e}"
            )

# Clean hovertemplate formatting
# <extra></extra> hides the redundant secondary trace-name box
htemp = (
    "x: %{x:.4f}<br>"
    "y: %{y:.4f}<br>"
    "z (def): %{z:.4f}<br>"
    "---<br>"
    "%{text}"
    "<extra></extra>" 
)

fig = go.Figure()

# Plot the in-plane layers (Z-planes)
for k in range(nz):
    fig.add_trace(go.Surface(
        x=X_def[:, :, k],
        y=Y_def[:, :, k],
        z=Z_def[:, :, k],
        surfacecolor=W[:, :, k],
        text=hover_text[:, :, k],
        hovertemplate=htemp,
        colorscale='Viridis',
        opacity=0.9,
        showscale=(k == 0), 
        name=f'Z-Layer {k}'
    ))

# Plot the Y-boundary faces (front and back)
for j in [0, ny - 1]:
    fig.add_trace(go.Surface(
        x=X_def[:, j, :],
        y=Y_def[:, j, :],
        z=Z_def[:, j, :],
        surfacecolor=W[:, j, :],
        text=hover_text[:, j, :],
        hovertemplate=htemp,
        colorscale='Viridis',
        showscale=False
    ))

# Plot the X-boundary faces (left and right)
for i in [0, nx - 1]:
    fig.add_trace(go.Surface(
        x=X_def[i, :, :],
        y=Y_def[i, :, :],
        z=Z_def[i, :, :],
        surfacecolor=W[i, :, :],
        text=hover_text[i, :, :],
        hovertemplate=htemp,
        colorscale='Viridis',
        showscale=False
    ))

# Configure layout 
fig.update_layout(
    title='Plate deflection using 3D elasticity (hover for displacements)',
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z',
        aspectmode='data'
    ),
    margin=dict(l=0, r=0, b=0, t=40)
)

# Force contour lines for the wireframe effect
fig.update_traces(
    contours=dict(
        x=dict(show=True, color='black', width=2),
        y=dict(show=True, color='black', width=2),
        z=dict(show=True, color='black', width=2)
    )
)

fig.show()