Practice, 3D semi-analytical plate buckling analysis via the Ritz Method¶

This notebook implements a semi-analytical framework to determine the critical buckling load of a 3D plate structure. It utilizes the Ritz method with Legendre polynomials as basis functions.

1Approximation apace and quadrature¶

We define the polynomial expansion orders ( $m_1, m_2, m_3$ ) for the displacement field $(u, v, w)$ along the $x, y, z$ axes. The total degrees of freedom is $N = 3 \times m_1 \times m_2 \times m_3$ . Integration points and weights are computed using Gauss-Legendre quadrature. The number of integration points ( $2m - 1$ ) is chosen to exactly integrate the stiffness terms formulated by the polynomial bases.

import numpy as np
from scipy.special import roots_legendre

from buckling.legendre import vecf, vecfxi

# Approximation order
m1 = m2 = 6
m3 = 5
N = 3 * m1 * m2 * m3

print(f'm1 = {m1}')
print(f'm2 = {m2}')
print(f'm3 = {m3}')
print(f'Total DOF (N) = {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)

m1 = 6
m2 = 6
m3 = 5
Total DOF (N) = 540

2Material, Geometry, and Boundary Conditions¶

Isotropic material properties and plate dimensions are defined. A fundamental compressive pre-stress state (sigxx = -1.0) is applied to evaluate the geometric stiffness.

Boundary constraints (xit1, xir1, etc.) are defined to enforce essential boundary conditions directly on the trial functions. These parameters act as boolean toggles or penalty flags within the vecf and vecfxi functions.

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

# Geometric properties
a = 0.3
b = 0.1
h = 0.003

# Fundamental stress for linear buckling
sigxx = -1.0

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

3Pre-allocation and helper functions¶

Pre-allocating arrays prevents memory fragmentation inside the deep integration loop. The addouter function performs an in-place outer product, avoiding temporary array creation during the assembly of $\mathbf{K}$ and $\mathbf{K}_G$ .

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

buff = np.zeros((N, N))
K = np.zeros((N, N))
KG = np.zeros((N, N))

def addouter(matrix, vec1, vec2):
    """In-place addition of outer product to a target matrix."""
    np.outer(vec1, vec2, out=buff)
    matrix += buff

43D Numerical integration (matrix assembly)¶

The linear stiffness matrix $\mathbf{K}$ and geometric stiffness matrix $\mathbf{K}_G$ are built via numerical integration over the 3D domain.

Instead of constructing explicit $\mathbf{B}$ and $\mathbf{D}$ matrices, the script leverages the tensor product nature of the Ritz bases. It computes the stress-like quantities (Sxx, Syy, etc.) evaluated at the integration points and outer-products them with the corresponding strain-displacement derivatives (exx, eyy, etc.).

The Jacobian determinant detJ = a*b*h/8 scales the integration from the $[-1, 1]^3$ parametric domain to the physical $[0, a] \times [0, b] \times [-h/2, h/2]$ domain.

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

            # Base shape functions
            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()

            # X-derivatives
            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()

            # Y-derivatives
            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()

            # Z-derivatives
            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()

            # Kinematics (Set NL = 0 for linear eigenvalue buckling)
            NL = 0 
            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)

            # 3D Elasticity Constitutive Relations
            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

            # Jacobian determinant for the physical domain transformation
            detJ = a * b * h / 8

            # Assemble Elastic Stiffness Matrix (K)
            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)

            # Assemble Geometric Stiffness Matrix (KG)
            addouter(KG, detJ*weight*sigxx*Swx, Swx)

5Eigenvalue problem solver¶

We solve the generalized eigenvalue problem $(\mathbf{K} - \lambda \mathbf{K}_G) \boldsymbol{\phi} = \mathbf{0}$ . The output $\lambda_{cr}$ represents the buckling load multiplier. The resulting eigenvalue is normalized against the classical plate buckling parameter $k_x = \frac{N_{xx} b^2}{\pi^2 D}$ for verification.

from scipy.sparse import csc_matrix
from structsolve import lb

# Using custom generalized eigenvalue solver
eigvals, eigvecs = lb(csc_matrix(K), csc_matrix(KG))

# Plate flexural rigidity
D = E * h**3 / (12 * (1 - nu**2))

# Non-dimensionalize the critical buckling load
lbd_cr = eigvals[0] * sigxx * b**2 * h / np.pi**2 / D

print(f"Normalized Buckling Load Factor: {lbd_cr:.6f}")
print(f"Critical Stress resultant (Nx): {sigxx * h}")

Running linear buckling analysis...
		Eigenvalue solver... 
			eigsh() solver...
				WARNING: Factor is exactly singular
			aborted!
			Removing null columns...
				300 columns removed
			finished!
			eigsh() solver...
			finished!
		finished!
	first 5 eigenvalues:
		2548232155.649012
		6452448898.655781
		14966202178.2332
		17032363524.681751
		31855249632.041416
Normalized Buckling Load Factor: -15.663520
Critical Stress resultant (Nx): -0.003

6Post-Processing: 3D visualization of the buckling mode shapes using plotly¶

This algorithm renders 3D deformed configurations of plate buckling modes using Plotly. It utilizes a tensor contraction approach for kinematic mapping, which is highly efficient for spectral or Ritz-based continuous displacement fields where spatial variables separate.

6.11D Domain Discretization and Natural Mapping¶

The script initializes by defining the structural domain. Instead of generating a full 3D mesh immediately, it isolates the physical axes into 1D vectors: $x$ , $y$ , and $z$ . These vectors are mapped linearly to natural parametric coordinates $\xi, \eta, \zeta \in [-1, 1]$ corresponding to the plate dimensions $a$ , $b$ , and $h$ .

6.2Basis Function Evaluation¶

The user-defined basis functions (vecf) are evaluated purely along the 1D parametric vectors.

Input: 1D parametric coordinates ( $\xi_{\text{1d}}$ , etc.) and modal parameters ( $m_1$ , etc.).
Output: Dense 2D arrays ( $P_\xi, P_\eta, P_\zeta$ ) containing the evaluated basis functions for every node along that specific axis.
Dimensions: $P_\xi$ has the shape $(N_x, m_1)$ , where $N_x$ is the spatial resolution and $m_1$ is the number of assumed modes/basis terms.

6.3Kinematic Field Assembly via Tensor Contraction¶

To map the generalized degrees of freedom (eigenvectors) to the physical continuous field, the algorithm avoids $O(N_x N_y N_z)$ loops.

Reshaping: The 1D eigenvector segment for a specific displacement component (e.g., $U$ ) is extracted and reshaped into a 3D coefficient tensor $C$ of shape $(m_1, m_2, m_3)$ .
Einstein Summation: np.einsum reconstructs the 3D spatial field via the following contraction:

U(x,y,z) = \sum_{i=1}^{m_1}\sum_{j=1}^{m_2}\sum_{k=1}^{m_3} C_{U,ijk} P_{\xi,i}(x) P_{\eta,j}(y) P_{\zeta,k}(z)

Syntax Mapping: The string 'ijk, xi, yj, zk -> xyz' instructs the C-backend to align the mode coefficient tensor (ijk) with the spatial evaluation matrices (xi, yj, zk) and output the full 3D spatial grid (xyz).

6.4Wireframe Topology Extraction¶

Plotly does not natively support volumetric structured grids. To render the mesh as a wireframe:

Deformation: The physical coordinates are perturbed by the scaled displacement fields (e.g., $X_{\text{def}} = X + \text{scale} \cdot U$ ).
Linearization: The 3D coordinate matrices are iterated over and flattened into 1D lists (lines_x, lines_y, lines_z).
Pen-up Delimiters: A None object is appended after every continuous nodal line segment along the $x$ , $y$ , and $z$ grid directions. This forces Plotly to break the line, preventing diagonal artifacts across the bounding box.

6.5Contour Mapping and Type Nullification¶

The out-of-plane displacement $W$ is mapped to the line color.

Type Enforcement: Plotly’s scatter3d.line.color strictly requires numerical arrays to apply the Jet colorscale. It throws an error if it encounters the None delimiters used in the spatial arrays.
Nullification: The float 0.0 is injected into the lines_w array wherever a None exists in the spatial arrays. Since the $X, Y, Z$ coordinates at that index are None, the rendering engine skips the vertex entirely, neutralizing the 0.0 so it does not distort the contour limits.

6.6Interactive State Management¶

Trace Generation: A separate go.Scatter3d trace is generated for each precomputed buckling mode. Only the first trace is set to visible=True at initialization.
Widget Logic: An updatemenus dropdown is built. Each button passes a specific boolean array (e.g., [True, False, False]) to the visible attribute of the figure, toggling the respective mode traces seamlessly without requiring a Python backend to recompute the matrices.

import numpy as np
import plotly.graph_objects as go

num_modes = 3  # Set the number of buckling modes to extract and toggle
scale = 100    # Deformation scaling factor for visibility

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'
)

# --- 1. Define 1D Coordinate Arrays ---
x_1d = np.linspace(0, a, nx)
y_1d = np.linspace(0, b, ny)
z_1d = np.linspace(-h/2, h/2, nz)

# Map to natural coordinates [-1, 1]
xi_1d = 2 * x_1d / a - 1
eta_1d = 2 * y_1d / b - 1
zeta_1d = 2 * z_1d / h - 1

# --- 2. Evaluate 1D Basis Functions ---
# Assuming vecf returns an array of shape (m,) for a scalar input.
# The resulting matrices will have shapes: (nx, m1), (ny, m2), (nz, m3)
P_xi = np.array([vecf(m1, xi, xit1, xir1, xit2, xir2) for xi in xi_1d])
P_eta = np.array([vecf(m2, eta, etat1, etar1, etat2, etar2) for eta in eta_1d])
P_zeta = np.array([vecf(m3, zeta, zetat1, zetar1, zetat2, zetar2) for zeta in zeta_1d])

# --- 3. Allocate Displacement Arrays ---
U = np.zeros((num_modes, nx, ny, nz))
V = np.zeros((num_modes, nx, ny, nz))
W = np.zeros((num_modes, nx, ny, nz))

n_dof = m1 * m2 * m3

# --- 4. Assemble 3D Fields via Tensor Contraction ---
for mode in range(num_modes):
    
    # Extract eigenvector segments and reshape them into 3D coefficient tensors
    # Order assumes eigvecs is packed as [U_dofs, V_dofs, W_dofs]
    c_U = eigvecs[0 : n_dof, mode].reshape((m1, m2, m3))
    c_V = eigvecs[n_dof : 2*n_dof, mode].reshape((m1, m2, m3))
    c_W = eigvecs[2*n_dof : 3*n_dof, mode].reshape((m1, m2, m3))
    
    # Contract the coefficients with the 1D basis matrices to yield the 3D spatial field
    # ijk: coefficient indices corresponding to m1, m2, m3
    # xi, yj, zk: point evaluation matrices (points vs basis index)
    # xyz: output target shape corresponding to nx, ny, nz
    U[mode] = np.einsum('ijk, xi, yj, zk -> xyz', c_U, P_xi, P_eta, P_zeta)
    V[mode] = np.einsum('ijk, xi, yj, zk -> xyz', c_V, P_xi, P_eta, P_zeta)
    W[mode] = np.einsum('ijk, xi, yj, zk -> xyz', c_W, P_xi, P_eta, P_zeta)

# Generate the 3D grid only for the final spatial plotting step
X, Y, Z = np.meshgrid(x_1d, y_1d, z_1d, indexing='ij')


fig = go.Figure()

# 2. Build traces for each precomputed mode
for mode in range(num_modes):
    X_def = X + scale * U[mode]
    Y_def = Y + scale * V[mode]
    Z_def = Z + scale * W[mode]
    W_amp = W[mode] # W displacement array for contour mapping
    
    lines_x, lines_y, lines_z, lines_w = [], [], [], []

    # X-direction lines
    for j in range(ny):
        for k in range(nz):
            lines_x.extend(X_def[:, j, k].tolist() + [None])
            lines_y.extend(Y_def[:, j, k].tolist() + [None])
            lines_z.extend(Z_def[:, j, k].tolist() + [None])
            lines_w.extend(W_amp[:, j, k].tolist() + [0.0]) # Valid float placeholder

    # Y-direction lines
    for i in range(nx):
        for k in range(nz):
            lines_x.extend(X_def[i, :, k].tolist() + [None])
            lines_y.extend(Y_def[i, :, k].tolist() + [None])
            lines_z.extend(Z_def[i, :, k].tolist() + [None])
            lines_w.extend(W_amp[i, :, k].tolist() + [0.0]) # Valid float placeholder

    # Z-direction lines
    for i in range(nx):
        for j in range(ny):
            lines_x.extend(X_def[i, j, :].tolist() + [None])
            lines_y.extend(Y_def[i, j, :].tolist() + [None])
            lines_z.extend(Z_def[i, j, :].tolist() + [None])
            lines_w.extend(W_amp[i, j, :].tolist() + [0.0]) # Valid float placeholder

    # Add the trace for the current mode
    fig.add_trace(go.Scatter3d(
        x=lines_x, 
        y=lines_y, 
        z=lines_z,
        mode='lines',
        line=dict(
            color=lines_w,            # Apply W amplitude as color
            colorscale='Jet',         # Standard structural mapping scale
            width=3,
            showscale=True,
            colorbar=dict(title=f"W (Mode {mode+1})", x=0.85)
        ),
        name=f'Mode {mode+1}',
        opacity=0.9,
        visible=(mode == 0)           # Only the first mode is visible initially
    ))

# 3. Create the toggle widget (Dropdown)
buttons = []
for mode in range(num_modes):
    # Boolean mask for Plotly visibility toggle
    visibility = [False] * num_modes
    visibility[mode] = True
    
    button = dict(
        label=f'Mode {mode+1}',
        method='update',
        args=[
            {'visible': visibility},
            {'title': f'Buckling Mode Shape: Mode {mode+1}'}
        ]
    )
    buttons.append(button)

# 4. Final layout configuration
fig.update_layout(
    updatemenus=[dict(
        active=0,
        buttons=buttons,
        x=0.05, xanchor='left',
        y=1.05, yanchor='top',
        direction='down',
        showactive=True
    )],
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z',
        aspectmode='data' 
    ),
    title='Buckling Mode Shape: Mode 1',
    showlegend=False,
    margin=dict(l=0, r=0, b=0, t=40)
)

fig.show()