femlabpy — Dynamics, Periodic BCs & Advanced Solvers Implementation Plan#

Author: Auto-generated implementation plan
Date: 2026-04-04
Target Version: v0.6.0
Total Tasks: 82 across 12 phases

Current State Summary#

femlabpy v0.5.0 is a Python FEM teaching library with:

Category

What Exists

Elements

T3 (CST), Q4 (bilinear), T4 (tet), H8 (hex), Bar/Truss — structural + scalar/potential

Solvers

Linear elastic, nonlinear bar (arc-length/orthogonal residual), elastoplastic (Von Mises + Drucker-Prager)

BCs

Direct elimination (setbc), Lagrange multiplier (solve_lag, solve_lag_general)

Materials

Linear elastic (plane stress/strain/3D), Von Mises plasticity, Drucker-Prager plasticity

Mesh I/O

Gmsh .msh parser (v2.x ASCII + v4.x binary), 19 element types

Post-process

Reaction recovery, stress/strain at Gauss points, matplotlib contour plots

Assembly

Dense + scipy sparse (COO scatter), vectorized (T3/T4/H8/Bar) and loop-based (Q4)

What’s Missing#

  • No mass matrices — none of the element routines compute element mass matrices

  • No damping — no Rayleigh damping or any damping formulation

  • No time integration — no Newmark-beta, central difference, HHT-alpha, or any transient solver

  • No modal analysis — no eigenvalue solver for natural frequencies / mode shapes

  • No periodic BCs — no multi-point constraint tying for periodic unit cells (RVE homogenization)

  • No consistent mass vs. lumped mass options

  • No time-dependent loading — no infrastructure for load functions P(t)

Implementation Order#

  1. Phase 1 + 10 → Mass matrices + new module infrastructure (foundation)

  2. Phase 6 → Modal analysis (validates mass matrices, needed for damping)

  3. Phase 2 → Damping models (required before dynamics)

  4. Phase 3 → Newmark-beta integrator (core dynamic solver)

  5. Phase 4 → Central difference / explicit solver

  6. Phase 5 → HHT-alpha method (refinement of Newmark)

  7. Phase 7 → Periodic boundary conditions (independent track)

  8. Phase 8 → Nonlinear dynamics

  9. Phase 9 → Dynamic post-processing & visualization

  10. Phase 11 + 12 → Examples & tests (validation throughout)

Phase 1: Mass Matrices#

Add element-level mass matrix routines (consistent + lumped) for all element types. Priority: HIGH | New files: element modules get mass functions

1.1 Bar Mass Matrix#

  • [ ] Function: mebar(Xe, Ge, *, lumped=False)

  • Location: src/femlabpy/elements/bars.py

  • Consistent (4×4 for 2D, 6×6 for 3D):

    M = (rho * A * L / 6) * [2, 1; 1, 2] ⊗ I_dof

    where rho = density, A = cross-section area, L = bar length, I_dof = identity of size dof

  • Lumped (diagonal):

    M = (rho * A * L / 2) * I_{2*dof}

  • Material table: G = [A, E, rho] (extend from [A, E])

  • Inputs: Xe = nodal coordinates (2, ndim), Ge = material row [A, E, rho]

  • Outputs: Me = element mass matrix (2*dof, 2*dof)

1.2 T3 (CST Triangle) Mass Matrix#

  • [ ] Function: met3e(Xe, Ge, *, lumped=False)

  • Location: src/femlabpy/elements/triangles.py

  • Consistent (6×6 for 2-DOF):

    M = (rho * t * A / 12) * [2, 1, 1;
                           1, 2, 1;
                           1, 1, 2] ⊗ I_2

    where A = triangle area, t = thickness

  • Lumped (diagonal):

    M = (rho * t * A / 3) * I_6

  • Material table: G = [E, nu, type, t, rho] (extend from [E, nu, type] or [E, nu])

  • Inputs: Xe = nodal coordinates (3, 2), Ge = material row

  • Outputs: Me = (6, 6) mass matrix

1.3 Q4 (Bilinear Quad) Mass Matrix#

  • [ ] Function: meq4e(Xe, Ge, *, lumped=False)

  • Location: src/femlabpy/elements/quads.py

  • Consistent (8×8) via 2×2 Gauss quadrature:

    M = Σ_{gp} rho * t * N^T * N * det(J) * w

    where N = shape function row vector [N1*I2, N2*I2, N3*I2, N4*I2]

  • Lumped via row-sum:

    M_L[i,i] = Σ_j M[i,j]    (then scale to preserve total mass)

  • Material table: G = [E, nu, type, t, rho]

  • Inputs: Xe = (4, 2), Ge = material row

  • Outputs: Me = (8, 8) mass matrix

1.4 T4 (Tetrahedron) Mass Matrix#

  • [ ] Function: meT4e(Xe, Ge, *, lumped=False)

  • Location: src/femlabpy/elements/solids.py

  • Consistent (12×12) — analytical formula:

    M = (rho * V / 20) * [2I, I, I, I;
                       I, 2I, I, I;
                       I, I, 2I, I;
                       I, I, I, 2I]

    where V = tet volume, I = 3×3 identity

  • Lumped (diagonal):

    M = (rho * V / 4) * I_12

  • Material table: G = [E, nu, rho] (extend from [E, nu])

  • Inputs: Xe = (4, 3), Ge = material row

  • Outputs: Me = (12, 12) mass matrix

1.5 H8 (Hexahedron) Mass Matrix#

  • [ ] Function: meh8e(Xe, Ge, *, lumped=False)

  • Location: src/femlabpy/elements/solids.py

  • Consistent (24×24) via 2×2×2 Gauss quadrature:

    M = Σ_{gp} rho * N^T * N * det(J) * w

    where N = (3, 24) shape function matrix with N_i * I_3 blocks

  • Lumped via HRZ (Hinton-Rock-Zienkiewicz) method:

    1. Compute consistent diagonal d_i = M_c[i,i]

    2. Scale: M_L[i,i] = d_i * (total_mass / Σ d_i)

  • Material table: G = [E, nu, rho]

  • Inputs: Xe = (8, 3), Ge = material row

  • Outputs: Me = (24, 24) mass matrix

1.6 Global Mass Assembly#

  • [ ] Function: assmm(M, Me, Te, dof)

  • Location: src/femlabpy/assembly.py

  • Identical scatter logic to assmk() but for mass matrix

  • Handles both dense np.ndarray and scipy.sparse.lil_matrix

  • Uses np.ix_ indexing for scatter operations

1.7 Vectorized Mass Assembly#

  • [ ] Functions: mt3e(M, T, X, G), mT4e(M, T, X, G), mh8e(M, T, X, G), mbar(M, T, X, G, dof=2)

  • Location: Respective element modules

  • Batch assembly using np.einsum matching existing vectorized k-routines

  • Sparse COO scatter for large meshes

  • Pattern: Same as kt3e, kT4e, kh8e, kbar but for mass

1.8 Material Table Extension#

  • [ ] Add rho (density) field to material table G

  • Backward compatibility: If rho not present in G row, default to rho=1.0 with a warning

  • Updated conventions:

    • Bar: [A, E] → [A, E, rho]

    • 2D structural: [E, nu] or [E, nu, type] → [E, nu, type, t, rho]

    • 3D structural: [E, nu] → [E, nu, rho]

  • Helper function: _extract_density(Ge, element_type) with smart defaults

1.9 Lumped Mass Option#

  • [ ] All mass routines accept lumped: bool = False keyword argument

  • When lumped=True:

    • Bar, T3, T4: use analytical lumped formulas (equal distribution)

    • Q4, H8: use row-sum or HRZ method

  • Document trade-offs: lumped = cheaper, no coupling, required for explicit methods; consistent = more accurate, better convergence for implicit

Phase 2: Damping Models#

Implement Rayleigh and other damping formulations. Priority: HIGH | New file: src/femlabpy/damping.py

2.1 Rayleigh Damping#

  • [ ] Function: rayleigh_damping(M, K, alpha, beta) -> C

  • Formula: C = alpha * M + beta * K

  • Convenience: rayleigh_coefficients(omega1, omega2, zeta1, zeta2) -> (alpha, beta)

    [alpha]   1        [omega2, -omega1] [zeta1]
[beta ] = ――――――――― [-1/omega2, 1/omega1] [zeta2]
          2

    Solves the 2×2 system for α, β given two target frequencies and damping ratios.

  • Support both dense and sparse matrices

  • Location: src/femlabpy/damping.py

2.3 Caughey Damping Series#

  • [ ] Function: caughey_damping(M, K, omega_targets, zeta_targets) -> C

  • Formula: C = M * Σ aₖ * (M⁻¹K)^k for k = 0, 1, …, n-1

  • Allows prescribing damping ratios at more than 2 frequencies

  • Coefficients aₖ found by solving a Vandermonde-like system

  • Priority: LOW

Phase 3: Newmark-Beta Time Integrator#

Core transient dynamic solver — the backbone of structural dynamics. Priority: HIGH | New file: src/femlabpy/dynamics.py

3.1 Time Integration Data Structures#

  • [ ] Dataclass: TimeHistory

    @dataclass
class TimeHistory:
    t: np.ndarray          # (nsteps+1,) time values
    u: np.ndarray          # (nsteps+1, ndof) displacement history
    v: np.ndarray          # (nsteps+1, ndof) velocity history
    a: np.ndarray          # (nsteps+1, ndof) acceleration history
    dt: float              # time step size
    nsteps: int            # number of time steps
    energy: dict | None    # {kinetic, strain, dissipated, external, total}

  • Dataclass: NewmarkParams

    @dataclass
class NewmarkParams:
    beta: float = 0.25     # Newmark beta parameter
    gamma: float = 0.5     # Newmark gamma parameter
    
    @classmethod
    def average_acceleration(cls): return cls(0.25, 0.5)
    
    @classmethod
    def linear_acceleration(cls): return cls(1/6, 0.5)
    
    @classmethod
    def central_difference(cls): return cls(0.0, 0.5)

3.2 Load Function Interface#

  • [ ] Abstract interface: load_func(t) -> p_vector

  • Built-in load functions:

    def constant_load(P):
    """Returns lambda t: P (constant in time)"""

def ramp_load(P, t_ramp):
    """Linear ramp: P * min(t/t_ramp, 1.0)"""

def harmonic_load(P, omega, phase=0.0):
    """Sinusoidal: P * sin(omega*t + phase)"""

def pulse_load(P, t_start, t_duration):
    """Rectangular pulse: P for t_start <= t <= t_start+t_duration, else 0"""

def tabulated_load(P, time_table, value_table):
    """Interpolated from table: P * interp1d(time_table, value_table)(t)"""

def seismic_load(M, direction, accel_record, dt_record):
    """Ground motion: p(t) = -M @ direction * a_g(t)"""

  • All return a callable f(t) -> ndarray(ndof, 1)

3.3 Newmark-Beta Solver (Main)#

  • [ ] Function: solve_newmark(M, C, K, p_func, u0, v0, dt, nsteps, *, beta=0.25, gamma=0.5, C_bc=None, dof=2) -> TimeHistory

  • Algorithm (implicit, average acceleration default):

    1. Compute initial acceleration: a0 = M^{-1} * (p(0) - C*v0 - K*u0)
2. Form effective stiffness: K_eff = K + M/(beta*dt^2) + gamma*C/(beta*dt)
3. Factorize K_eff (once, since linear)
4. For each step n = 0, 1, ..., nsteps-1:
   a. Compute effective load:
      p_eff = p(t_{n+1}) 
            + M * [u_n/(beta*dt^2) + v_n/(beta*dt) + (1/(2*beta)-1)*a_n]
            + C * [gamma*u_n/(beta*dt) + (gamma/beta-1)*v_n + dt*(gamma/(2*beta)-1)*a_n]
   b. Solve: K_eff * u_{n+1} = p_eff
   c. Update acceleration: a_{n+1} = (u_{n+1}-u_n)/(beta*dt^2) - v_n/(beta*dt) - (1/(2*beta)-1)*a_n
   d. Update velocity: v_{n+1} = v_n + dt*[(1-gamma)*a_n + gamma*a_{n+1}]
   e. Store in TimeHistory

3.4 Effective Stiffness Formation#

  • [ ] Function: _newmark_effective_stiffness(K, M, C, dt, beta, gamma) -> K_eff

  • Formula:

    K_eff = K + (1 / (beta * dt^2)) * M + (gamma / (beta * dt)) * C

  • Handles sparse and dense matrices

  • Returns factorized form (scipy.sparse.linalg.splu or scipy.linalg.lu_factor) for repeated solves

3.5 Newmark Update Formulas#

  • [ ] Function: _newmark_update(u_new, u_old, v_old, a_old, dt, beta, gamma) -> (a_new, v_new)

  • Acceleration update:

    a_{n+1} = (1/(beta*dt^2)) * (u_{n+1} - u_n) - (1/(beta*dt)) * v_n - (1/(2*beta) - 1) * a_n

  • Velocity update:

    v_{n+1} = v_n + dt * [(1 - gamma) * a_n + gamma * a_{n+1}]

3.6 Newmark Parameter Presets#

  • [ ] Named presets as class methods on NewmarkParams:

    • average_acceleration() → (beta=0.25, gamma=0.5) — unconditionally stable, no numerical dissipation, 2nd order accurate

    • linear_acceleration() → (beta=1/6, gamma=0.5) — conditionally stable (dt < 0.551 * T_min), 2nd order accurate

    • central_difference() → (beta=0, gamma=0.5) — explicit, conditionally stable (dt < T_min/π), 2nd order accurate

    • fox_goodwin() → (beta=1/12, gamma=0.5) — conditionally stable, 4th order accurate for undamped SDOF

  • Document stability regions and accuracy order for each

3.7 Adaptive Time Stepping#

  • [ ] Function: solve_newmark_adaptive(M, C, K, p_func, u0, v0, dt0, t_end, *, tol=1e-4) -> TimeHistory

  • Error estimator: Run dual solutions with beta=0.25 and beta=1/6, compare:

    err = ||u_025 - u_016|| / ||u_025||

  • Step adjustment: dt_new = dt * (tol / err)^(1/3) with safety factor 0.9

  • Priority: LOW (nice-to-have)

3.8 Boundary Conditions in Dynamic Context#

  • [ ] Function: _apply_dynamic_bc(K_eff, p_eff, C_bc, dof, t, ks)

  • Handle three BC types:

    1. Fixed zero: Same as static setbc — zero row/col, penalty diagonal

    2. Prescribed displacement history: u_bc(t) — modify effective load:

      p_eff[free] -= K_eff[free, bc] * u_bc(t)

    3. Velocity BC: Convert to displacement BC via integration

  • BC constraint matrix C_bc can include time-dependent values via callable

  • Re-use existing setbc pattern where possible

3.9 Energy Balance Check#

  • [ ] Function: _compute_energy(M, C, K, u, v, p_func, t, dt) -> dict

  • Quantities at each step:

    E_kinetic    = 0.5 * v^T * M * v
E_strain     = 0.5 * u^T * K * u
E_dissipated += 0.5 * (v_n + v_{n+1})^T * C * (v_n + v_{n+1}) * dt  (incremental)
W_external   += 0.5 * (p_n + p_{n+1})^T * (u_{n+1} - u_n)           (incremental)
E_total      = E_kinetic + E_strain + E_dissipated

  • Validation: For undamped system (C=0), E_total should be constant

  • Store in TimeHistory.energy dict

3.10 Stability Check#

  • [ ] Function: _check_stability(K, M, dt, beta, gamma) -> (is_stable, dt_critical)

  • For conditionally stable schemes (beta < 0.25 or gamma != 0.5):

    omega_max = sqrt(max(eig(M^{-1} K)))      # maximum natural frequency
dt_critical = 2 / omega_max               # for central difference

    For linear acceleration: dt_critical = 0.551 * (2*pi / omega_max)

  • Warning: Emit UserWarning if dt > dt_critical

  • For unconditionally stable schemes (beta >= gamma/2 >= 0.25): no warning needed

Phase 4: Central Difference (Explicit) Solver#

For wave propagation, impact, and blast loading. Priority: MEDIUM | Location: src/femlabpy/dynamics.py

4.1 Central Difference Solver#

  • [ ] Function: solve_central_diff(M_lumped, C, K, p_func, u0, v0, dt, nsteps, *, dof=2) -> TimeHistory

  • Requirements: M_lumped must be diagonal (1D array or diagonal matrix)

  • Algorithm:

    1. Initial: a0 = M_L^{-1} * (p(0) - C*v0 - K*u0)
2. Compute: u_{-1} = u0 - dt*v0 + 0.5*dt^2*a0  (backward extrapolation)
3. For each step n = 0, 1, ..., nsteps-1:
   a. p_eff = p(t_n) - (K - 2*M_L/dt^2)*u_n - (M_L/dt^2 - C/(2*dt))*u_{n-1}
   b. u_{n+1} = (M_L/dt^2 + C/(2*dt))^{-1} * p_eff
   c. Velocity: v_n = (u_{n+1} - u_{n-1}) / (2*dt)
   d. Acceleration: a_n = (u_{n+1} - 2*u_n + u_{n-1}) / dt^2

  • Key advantage: No matrix factorization when M is diagonal and C is also diagonal (or zero)

4.2 Central Difference Algorithm Details#

  • [ ] Implement the three-level time marching scheme

  • Handle the special initial step (u_{-1} computation)

  • For diagonal M and diagonal C: pure element-by-element, O(ndof) per step

  • For diagonal M and full C: still requires M/C combination but no full solve

4.3 Critical Time Step Computation#

  • [ ] Function: critical_timestep(K, M, *, method='power') -> dt_cr

  • Methods:

    • 'power': Power iteration to find omega_max = sqrt(lambda_max(M^{-1}K)) then dt_cr = 2 / omega_max

    • 'element': Estimate from smallest element size and wave speed dt_cr = h_min / c where c = sqrt(E/rho) and h_min = min element characteristic length

  • Usage: Called automatically by solve_central_diff with warning if dt > dt_cr

4.4 Lumped Mass Validation#

  • [ ] Auto-check in solve_central_diff:

    • If M is not diagonal, raise ValueError("Central difference requires lumped (diagonal) mass matrix")

    • Provide helper: is_lumped(M) -> bool — checks if off-diagonal entries are zero (or negligible)

    • Suggest: "Use lumped=True in mass matrix computation"

Phase 5: HHT-Alpha (Hilber-Hughes-Taylor) Method#

Numerical dissipation control for high-frequency noise. Priority: MEDIUM | Location: src/femlabpy/dynamics.py

5.1 HHT-Alpha Solver#

  • [ ] Function: solve_hht(M, C, K, p_func, u0, v0, dt, nsteps, *, alpha=-0.05) -> TimeHistory

  • Valid range: alpha ∈ [-1/3, 0]

    • alpha = 0 → standard Newmark (trapezoidal)

    • alpha = -0.05 → mild dissipation (recommended default)

    • alpha = -1/3 → maximum dissipation

5.2 HHT Parameters and Equilibrium#

  • [ ] Derived Newmark parameters:

    beta  = (1 - alpha)^2 / 4
gamma = 0.5 - alpha

  • Modified equilibrium equation at step n+1:

    M * a_{n+1} + (1+alpha) * C * v_{n+1} - alpha * C * v_n
             + (1+alpha) * K * u_{n+1} - alpha * K * u_n
             = (1+alpha) * p_{n+1} - alpha * p_n

  • Effective stiffness:

    K_eff = (1+alpha)*K + M/(beta*dt^2) + (1+alpha)*gamma*C/(beta*dt)

  • Effective load includes alpha-weighted previous step terms

5.3 Recovery & Validation#

  • [ ] Verify that alpha=0 produces identical results to solve_newmark(beta=0.25, gamma=0.5)

  • Demonstrate numerical dissipation: run high-frequency oscillation and show amplitude decay increases with |alpha|

  • Show that low-frequency response is preserved (2nd order accuracy maintained)

Phase 6: Modal Analysis (Eigenvalue Solver)#

Natural frequencies and mode shapes. Priority: HIGH | New file: src/femlabpy/modal.py

6.1 Eigenvalue Solver#

  • [ ] Function: solve_modal(K, M, n_modes=10, *, C_bc=None, dof=2, sigma=0.0) -> ModalResult

  • Dataclass: ModalResult

    @dataclass
class ModalResult:
    eigenvalues: np.ndarray      # (n_modes,) omega^2 values
    omega: np.ndarray            # (n_modes,) natural frequencies in rad/s
    freq_hz: np.ndarray          # (n_modes,) frequencies in Hz
    period: np.ndarray           # (n_modes,) periods in seconds
    mode_shapes: np.ndarray      # (ndof, n_modes) eigenvectors (columns)
    participation: np.ndarray    # (n_modes, ndim) modal participation factors
    effective_mass: np.ndarray   # (n_modes, ndim) effective modal mass

  • Solver: scipy.sparse.linalg.eigsh(K, k=n_modes, M=M, sigma=sigma, which='SM') (shift-invert mode for smallest eigenvalues)

6.2 Boundary Condition Handling in Eigenvalue Problem#

  • [ ] Strategy: System Reduction

    def _reduce_system(K, M, C_bc, dof):
    """Eliminate constrained DOFs before eigensolve."""
    free_dofs = _get_free_dofs(C_bc, dof, ndof)
    K_red = K[np.ix_(free_dofs, free_dofs)]
    M_red = M[np.ix_(free_dofs, free_dofs)]
    return K_red, M_red, free_dofs

  • After eigensolve, expand mode shapes back to full DOF space (zeros at constrained DOFs)

  • Alternative: Large penalty on constrained DOFs: K[j,j] += 1e20 * max(diag(K)) (simpler but less accurate)

6.3 Mode Shape Normalization#

  • [ ] Mass normalization (default):

    phi_i = phi_i / sqrt(phi_i^T * M * phi_i)

    so that Phi^T * M * Phi = I

  • Displacement normalization:

    phi_i = phi_i / max(|phi_i|)

    so that maximum component is 1.0

6.4 Frequency Output#

  • [ ] Compute and store:

    omega_i = sqrt(lambda_i)         # rad/s
f_i     = omega_i / (2 * pi)     # Hz
T_i     = 1 / f_i                # seconds (period)

  • Print formatted table:

    Mode  omega (rad/s)   f (Hz)    T (s)
----  ------------   --------  --------
  1      12.34        1.964     0.509
  2      48.67        7.747     0.129
  ...

6.6 Mode Shape Plotting#

  • [ ] Function: plot_modes(T, X, phi, dof, mode_indices=None, *, scale=1.0, rows=None, cols=None)

  • Subplot grid showing deformed mesh for each requested mode

  • Overlay undeformed (dashed) and deformed (solid) mesh

  • Title: “Mode {i}: f = {f_i:.3f} Hz”

  • Optional: colormap by displacement magnitude

Phase 7: Periodic Boundary Conditions#

For RVE/unit cell analysis and computational homogenization. Priority: HIGH | New file: src/femlabpy/periodic.py

7.1 Periodic Node Pair Identification#

  • [ ] Function: find_periodic_pairs(X, axis, tol=1e-6) -> ndarray

  • Algorithm:

    1. Find bounding box: x_min, x_max along specified axis
2. Select nodes on "left" face: X[:, axis] ≈ x_min
3. Select nodes on "right" face: X[:, axis] ≈ x_max
4. For each left node, find matching right node by remaining coordinates
5. Return pairs as (n_pairs, 2) array: [[left_node, right_node], ...]

  • Matching: KD-tree or brute-force distance for remaining coordinate dimensions

  • Tolerance: tol relative to element size for floating-point robustness

  • Output: pairs = [[n_left_1, n_right_1], [n_left_2, n_right_2], ...] (1-based node numbers)

7.2 Multi-Axis Periodicity#

  • [ ] Function: find_all_periodic_pairs(X, periodic_axes, tol=1e-6) -> dict

  • Support:

    • 1D: periodic_axes=[0] → x-periodic only

    • 2D: periodic_axes=[0, 1] → x+y periodic

    • 3D: periodic_axes=[0, 1, 2] → fully periodic

  • Returns dict: {0: pairs_x, 1: pairs_y, 2: pairs_z}

  • Handle corner nodes that appear in multiple pair sets (avoid double-constraining)

  • Handle edge nodes in 3D (on intersection of two periodic faces)

7.3 Periodic Constraint Matrix#

  • [ ] Function: periodic_constraints(X, pairs, dof, *, eps_macro=None) -> (G, Q)

  • Constraint equation for each pair (L, R):

    u_R - u_L = eps_macro * (x_R - x_L)

    In matrix form: G * u = Q

    • G[row, R_dof] = +1, G[row, L_dof] = -1

    • Q[row] = eps_macro * (x_R - x_L) projected onto the DOF direction

  • If eps_macro=None: Purely periodic (zero fluctuation): Q = 0, constraints become u_R = u_L

  • Size: G is (n_constraints, n_total_dof), Q is (n_constraints, 1)

  • n_constraints = n_pairs * dof

7.4 Integration with Lagrange Multiplier Solver#

  • [ ] Function: solve_periodic(K, p, X, T, G, dof, *, eps_macro=None) -> (u, lam)

  • Build G, Q from periodic_constraints()

  • Pass to existing solve_lag_general(K, p, G, Q) to solve the augmented system:

    [K    G^T] [u]   [p]
[G    0  ] [λ] = [Q]

  • Return displacement field u and Lagrange multipliers λ (constraint forces = tractions)

7.5 MPC Alternative (Direct Elimination)#

  • [ ] Function: apply_mpc(K, p, pairs, dof) -> (K_red, p_red, slave_dofs)

  • Strategy: For each pair, eliminate “right” (slave) DOF:

    u_slave = u_master + Q/G

    Condense K and p by substituting slave DOFs

  • Advantage: Smaller system, no Lagrange multipliers

  • Disadvantage: More complex implementation, harder to extract reactions

  • Priority: MEDIUM

7.6 Macro Strain Application#

  • [ ] Function: apply_macro_strain(X, pairs, eps_macro, dof) -> Q

  • 2D macro strain tensor (Voigt): eps_macro = [exx, eyy, gamma_xy]

    eps_tensor = [[exx,       gamma_xy/2],
              [gamma_xy/2, eyy       ]]

  • 3D macro strain tensor (Voigt): eps_macro = [exx, eyy, ezz, gamma_xy, gamma_yz, gamma_xz]

  • RHS for each pair:

    Q_pair = eps_tensor @ (x_R - x_L)   → vector of size dof

  • Assembled into full Q vector matching constraint ordering

7.7 Computational Homogenization#

  • [ ] Function: homogenize(K, T, X, G, pairs, dof) -> C_eff

  • Algorithm:

    1. For each canonical unit strain e_i (i = 1..n_voigt):
   a. Build G, Q_i = apply_macro_strain(X, pairs, e_i, dof)
   b. Solve: u_i = solve_periodic(K, p=0, X, T, G, dof, eps_macro=e_i)
   c. Compute: sigma_avg_i = volume_average_stress(T, X, G, u_i, dof)
2. Assemble C_eff column-wise: C_eff[:, i] = sigma_avg_i

  • 2D: 3 load cases (exx, eyy, gamma_xy) → 3×3 C_eff

  • 3D: 6 load cases → 6×6 C_eff

  • Verify: C_eff should be symmetric and positive definite

7.8 Volume Averaging Utilities#

  • [ ] Function: volume_average_stress(T, X, G, u, dof, *, element_type='q4') -> sigma_avg

    sigma_avg = (1/V_total) * Σ_e sigma_e * V_e

    where sigma_e = stress at element centroid (or Gauss point average), V_e = element volume/area

  • [ ] Function: volume_average_strain(T, X, G, u, dof, *, element_type='q4') -> eps_avg

    eps_avg = (1/V_total) * Σ_e eps_e * V_e

  • Support T3, Q4, T4, H8 elements

  • For elements with multiple Gauss points (Q4, H8), average over Gauss points first

7.9 Corner Node Fixing#

  • [ ] Function: fix_corner(X, C, dof) -> C_extended

  • Find the node at the minimum coordinates (bottom-left corner for 2D)

  • Add zero displacement constraints for all DOFs at that node

  • Removes rigid body translation while preserving periodicity

  • For 2D: fix 2 DOFs (remove translation)

  • For 3D: fix 3 DOFs at corner + additional rotational fix if needed

7.10 Periodic Mesh Validation#

  • [ ] Function: check_periodic_mesh(X, T, axis, tol=1e-6) -> (is_valid, report)

  • Checks:

    1. Equal number of nodes on opposite faces

    2. Matching coordinate patterns (sorted by remaining coordinates)

    3. Compatible element topology on faces (matching edge patterns)

    4. No orphan nodes (not paired)

  • Report: Human-readable dict with pass/fail for each check + suggested fixes

  • Warn if mesh is non-periodic and suggest re-meshing with Gmsh structured mesh

Phase 8: Nonlinear Dynamics#

Combine existing nonlinear solvers with time integration. Priority: MEDIUM | Location: src/femlabpy/dynamics.py

8.1 Nonlinear Newmark Solver#

  • [ ] Function: solve_newmark_nl(M, C, tangent_func, internal_force_func, p_func, u0, v0, dt, nsteps, *, beta=0.25, gamma=0.5, tol=1e-6, max_iter=20) -> TimeHistory

  • Arguments:

    • tangent_func(u, state) -> K_t — computes tangent stiffness at current config

    • internal_force_func(u, state) -> (q, state_new) — computes internal forces + updates state

    • state — can carry plastic history (S, E) or geometric config

  • Algorithm at each time step:

    1. Predictor: u_pred = u_n + dt*v_n + 0.5*dt^2*(1-2*beta)*a_n
              v_pred = v_n + dt*(1-gamma)*a_n
2. Newton iteration k = 0, 1, ...:
   a. K_t = tangent_func(u^{k}, state)
   b. q = internal_force_func(u^{k}, state)
   c. K_eff = K_t + M/(beta*dt^2) + gamma*C/(beta*dt)
   d. R = p(t_{n+1}) - M*a^{k} - C*v^{k} - q
   e. Δu = K_eff^{-1} * R
   f. u^{k+1} = u^{k} + Δu
   g. Update a, v from Newmark formulas
   h. Check convergence: ||R|| / ||p|| < tol
3. Accept: u_{n+1} = u^{k+1}, update state

8.2 Tangent-Based Iteration Details#

  • [ ] Effective tangent: K_eff(u) = K_t(u) + a0*M + a1*C where a0 = 1/(beta*dt^2), a1 = gamma/(beta*dt)

  • Option: modified_newton=True → re-use K_t from first iteration (cheaper, slower convergence)

  • Option: initial_stiffness=True → use K_0 (elastic) throughout (simplest, worst convergence)

8.3 Convergence Criteria#

  • [ ] Multiple criteria (configurable):

    1. Residual norm: ||R||₂ / ||p(t_{n+1})||₂ < tol_r

    2. Energy norm: |ΔuᵀR| / |u_{n+1}ᵀp_{n+1}| < tol_e

    3. Displacement increment: ||Δu||₂ / ||u_{n+1}||₂ < tol_u

  • Default: use residual norm with tol_r = 1e-6

  • Log iteration history: {step, iter, ||R||, ||Δu||, energy_error}

8.4 Line Search and Arc-Length#

  • [ ] Line search: After Newton direction Δu, find optimal step η ∈ (0, 1]:

    min η: G(η) = Δu^T * R(u + η*Δu)

    Using bisection or cubic interpolation

  • Arc-length within time steps: For highly nonlinear snap-through in dynamics

    • Constraint: ||Δu||^2 + (β*ψ*dt)^2*||Δλ*p||^2 = Δs^2

    • Rarely needed but important for explosive/impact problems

  • Priority: LOW

Phase 9: Dynamic Post-Processing & Visualization#

Priority: MEDIUM | Location: src/femlabpy/plotting.py (extend)

9.1 Time History Plots#

  • [ ] Function: plot_time_history(result: TimeHistory, dof_index: int | list, *, quantity='displacement', ax=None)

  • Quantities: 'displacement', 'velocity', 'acceleration'

  • Support single DOF (line plot) or multiple DOFs (overlaid or subplot)

  • Auto-label axes: “Displacement (m)”, “Time (s)”, etc.

  • Return matplotlib axes for further customization

9.2 Frequency Response Function (FRF)#

  • [ ] Function: compute_frf(M, C, K, input_dof, output_dof, freq_range, *, n_points=500) -> (freq, H)

  • Formula: At each frequency ω:

    H(ω) = [(-ω^2*M + iω*C + K)^{-1}]_{output_dof, input_dof}

  • Plotting: plot_frf(freq, H, *, log_scale=True) — magnitude and phase vs frequency

  • Mark peaks (resonances) with vertical lines and frequency labels

  • Support dB scale: 20*log10(|H|/|H_ref|)

9.3 Dynamic Animation#

  • [ ] Function: animate_dynamic(T, X, u_history, dof, *, dt=1.0, scale=1.0, interval=50, save_path=None)

  • Uses matplotlib.animation.FuncAnimation

  • Each frame: deformed mesh at time step n

  • Optional: color by displacement magnitude or stress

  • Save to .gif or .mp4 via save_path

  • Show time counter in title

9.4 Energy Plots#

  • [ ] Function: plot_energy(result: TimeHistory, *, ax=None)

  • Plot on same axes: E_kinetic (blue), E_strain (red), E_dissipated (green), W_external (black dashed), E_total (thick black)

  • Verify: for undamped, E_total should be flat (horizontal line)

  • Legend with max values

9.5 Phase Plots#

  • [ ] Function: plot_phase(result: TimeHistory, dof_index: int, *, ax=None)

  • Plot velocity vs displacement at specified DOF

  • For SDOF harmonic: should show ellipse

  • For nonlinear: shows limit cycle or chaotic attractor

  • Priority: LOW

9.6 FFT of Response#

  • [ ] Function: fft_response(result: TimeHistory, dof_index: int, *, n_fft=None, window='hann') -> (freq, amplitude)

  • Apply window function, compute np.fft.rfft

  • Plot amplitude spectrum with peaks labeled

  • Compare peaks with modal frequencies (if available)

  • Priority: LOW

Phase 10: Infrastructure & Integration#

Priority: HIGH

10.1 Create src/femlabpy/dynamics.py#

  • [ ] New module housing:

    • TimeHistory and NewmarkParams dataclasses

    • solve_newmark() — implicit Newmark-beta solver

    • solve_central_diff() — explicit central difference

    • solve_hht() — HHT-alpha method

    • solve_newmark_nl() — nonlinear Newmark

    • Load function builders (constant, ramp, harmonic, pulse, tabulated, seismic)

    • critical_timestep() — critical dt computation

    • Internal helpers (_newmark_effective_stiffness, _newmark_update, _apply_dynamic_bc, _compute_energy, _check_stability)

10.2 Create src/femlabpy/modal.py#

  • [ ] New module housing:

    • ModalResult dataclass

    • solve_modal() — eigenvalue solver

    • _reduce_system() — BC elimination for eigenproblems

    • _modal_participation() — participation factors

    • plot_modes() — mode shape visualization

10.3 Create src/femlabpy/periodic.py#

  • [ ] New module housing:

    • find_periodic_pairs() — node pair identification

    • find_all_periodic_pairs() — multi-axis periodicity

    • periodic_constraints() — constraint matrix builder

    • solve_periodic() — Lagrange multiplier periodic solver

    • apply_mpc() — direct elimination alternative

    • apply_macro_strain() — RHS from macro strain

    • homogenize() — full homogenization pipeline

    • volume_average_stress() / volume_average_strain() — averaging

    • fix_corner() — rigid body removal

    • check_periodic_mesh() — mesh validation

10.4 Extend init() in core.py#

  • [ ] Updated signature: init(nn, dof, *, dynamic=False, use_sparse=None)

  • When dynamic=True: return (K, M, p, q) — allocate mass matrix M in addition to K

  • Damping matrix C is not pre-allocated (constructed from M, K via Rayleigh coefficients)

  • Backward compatible: dynamic=False (default) returns (K, p, q) as before

10.5 Update __init__.py Exports#

  • [ ] Add all new public symbols to __all__:

    # Mass matrices
from .elements.bars import mebar, mbar
from .elements.triangles import met3e, mt3e
from .elements.quads import meq4e, mq4e
from .elements.solids import meT4e, mT4e, meh8e, mh8e

# Damping
from .damping import rayleigh_damping, rayleigh_coefficients, modal_damping

# Dynamics
from .dynamics import (solve_newmark, solve_central_diff, solve_hht,
                       solve_newmark_nl, TimeHistory, NewmarkParams,
                       constant_load, ramp_load, harmonic_load, pulse_load,
                       tabulated_load, seismic_load, critical_timestep)

# Modal
from .modal import solve_modal, ModalResult, plot_modes

# Periodic
from .periodic import (find_periodic_pairs, find_all_periodic_pairs,
                       periodic_constraints, solve_periodic, apply_macro_strain,
                       homogenize, volume_average_stress, volume_average_strain,
                       check_periodic_mesh, fix_corner)

10.6 Update CLI --info#

  • [ ] Extend __main__.py to list:

    • Dynamic solvers: Newmark-beta, central difference, HHT-alpha, nonlinear Newmark

    • Modal analysis capabilities

    • Periodic BC support

    • Mass matrix element types

    • Damping models

  • Priority: LOW

Phase 11: Examples & Validation#

Priority: HIGH for core examples, MEDIUM for advanced

11.1 SDOF Spring-Mass-Damper (Newmark Validation)#

  • [ ] File: src/femlabpy/examples/dynamic_sdof.py

  • Setup: Single DOF: m=1, k=100, c=2, p(t) = 10*sin(5t)

  • Analytical solution: u(t) = A*sin(ωt - φ) + transient where ω_n = sqrt(k/m), ζ = c/(2*m*ω_n), ω_d = ω_n*sqrt(1-ζ²)

  • Verification: Compare Newmark solution with analytical at all time steps

  • Convergence: Show error decreases as O(dt²)

11.2 Cantilever Beam Free Vibration (Modal Validation)#

  • [ ] File: src/femlabpy/examples/dynamic_cantilever.py

  • Setup: Cantilever beam (Q4 elements), fixed left end, free right end

  • Analytical frequencies (Euler-Bernoulli):

    f_n = (beta_n * L)^2 / (2*pi*L^2) * sqrt(EI / (rho*A))

    where beta_1*L = 1.875, beta_2*L = 4.694, beta_3*L = 7.855, ...

  • Procedure: Assemble K and M, run solve_modal(), compare first 3 frequencies

  • Expected accuracy: Within 5% for coarse Q4 mesh, <1% for refined mesh

11.3 Bar Wave Propagation (Explicit Validation)#

  • [ ] File: src/femlabpy/examples/dynamic_wave.py

  • Setup: Long bar, impulse at one end, free at other

  • Analytical: Wave travels at c = sqrt(E/rho), arrives at other end at t = L/c

  • Solver: solve_central_diff with lumped mass

  • Verification: Wave front position matches analytical speed

  • Priority: MEDIUM

11.4 2D RVE with Circular Inclusion (Homogenization)#

  • [ ] File: src/femlabpy/examples/periodic_rve.py

  • Setup: Square unit cell (Q4 mesh) with circular inclusion of different material

  • Materials: Matrix (E1, nu1), Inclusion (E2, nu2), volume fraction f

  • Procedure: homogenize() → C_eff

  • Validation: Compare with Mori-Tanaka prediction:

    E_eff = E_m * (1 + f*(E_i/E_m - 1) / (1 + (1-f)*(E_i/E_m - 1)*S))

    and check that C_eff lies within Hashin-Shtrikman bounds

11.5 Periodic Unit Cell Under Shear#

  • [ ] File: src/femlabpy/examples/periodic_shear.py

  • Setup: Homogeneous unit cell, apply macro shear strain gamma_xy = 0.01

  • Expected result: Uniform shear stress throughout, sigma_xy = G * gamma_xy

  • Verification: Volume-averaged stress matches analytical, displacement field is affine

  • Priority: MEDIUM

11.6 Forced Vibration of Plate#

  • [ ] File: src/femlabpy/examples/dynamic_plate.py

  • Setup: Simply-supported plate (Q4), harmonic point load at center

  • Sweep frequency and compute response amplitude

  • Verification: FRF peaks coincide with modal analysis frequencies

  • Priority: MEDIUM

11.7 Seismic Ground Motion Response#

  • [ ] File: src/femlabpy/examples/dynamic_seismic.py

  • Setup: Multi-story frame (bar/beam elements) subjected to El Centro earthquake

  • Load: p(t) = -M @ r * a_g(t) where a_g(t) = ground acceleration record

  • Solver: solve_newmark with Rayleigh damping (5% at modes 1 and 3)

  • Output: Roof displacement time history, max inter-story drift

  • Priority: MEDIUM

Phase 12: Tests#

Priority: HIGH for core, MEDIUM for advanced

12.1 tests/test_mass_matrices.py#

  • [ ] Tests:

    • test_bar_mass_consistent — verify shape, symmetry, positive semi-definite

    • test_bar_mass_lumped — verify diagonal, total mass = rhoAL

    • test_t3_mass_consistent — verify shape (6,6), symmetry, positive definite

    • test_t3_mass_total — sum of diagonal = total mass = rhotA

    • test_q4_mass_consistent — verify (8,8), symmetry, SPD

    • test_q4_mass_lumped_preserves_mass — row-sum lumped preserves total mass

    • test_t4_mass_consistent — verify (12,12), total mass = rho*V

    • test_h8_mass_consistent — verify (24,24), total mass = rho*V

    • test_global_mass_assembly — verify assembled M has correct total mass

    • test_mass_sparse_dense_agree — sparse and dense assembly give same result

12.2 tests/test_newmark.py#

  • [ ] Tests:

    • test_sdof_undamped_free — compare with u(t) = u0*cos(omega_n*t), error < 1e-3

    • test_sdof_damped_forced — compare with analytical forced response

    • test_energy_conservation_undamped — E_total variation < 1e-6 over 1000 steps

    • test_unconditional_stability — large dt still converges (doesn’t blow up)

    • test_2nd_order_accuracy — error ratio at dt vs dt/2 ≈ 4.0

    • test_linear_acceleration_preset — beta=1/6 produces correct results

    • test_zero_initial_conditions — static load converges to static solution

    • test_bc_handling — constrained DOFs remain zero throughout

12.3 tests/test_central_diff.py#

  • [ ] Tests:

    • test_sdof_explicit_matches_implicit — for small dt, same as Newmark

    • test_conditional_stability_blowup — dt > dt_cr causes amplitude growth

    • test_critical_timestep_computation — dt_cr matches analytical for SDOF

    • test_requires_lumped_mass — error raised for consistent mass input

    • test_wave_speed — wave front position in bar matches c = sqrt(E/rho)

    • Priority: MEDIUM

12.4 tests/test_modal.py#

  • [ ] Tests:

    • test_sdof_frequency — omega = sqrt(k/m) for single spring-mass

    • test_bar_frequencies — compare with analytical f_n = n/(2L) * sqrt(E/rho)

    • test_mass_orthogonality — Phi^T M Phi = I within tolerance

    • test_stiffness_orthogonality — Phi^T K Phi = diag(omega_i^2)

    • test_bc_elimination — constrained DOFs have zero in mode shapes

    • test_participation_sum — sum of effective masses = total mass (per direction)

    • test_n_modes_parameter — requesting n modes returns exactly n

12.5 tests/test_periodic.py#

  • [ ] Tests:

    • test_pair_detection_simple — 4x4 mesh, known pairs on boundaries

    • test_pair_detection_symmetric_count — left face nodes = right face nodes

    • test_constraint_matrix_rank — G has correct rank (n_pairs * dof)

    • test_homogeneous_unit_strain — uniform material under unit strain gives expected C

    • test_C_eff_symmetry — C_eff = C_eff^T within tolerance

    • test_C_eff_positive_definite — all eigenvalues of C_eff > 0

    • test_periodic_vs_large_domain — C_eff matches large domain (homogeneous case)

    • test_corner_fixing — solution doesn’t have rigid body modes

12.6 tests/test_hht.py#

  • [ ] Tests:

    • test_alpha_zero_is_newmark — solve_hht(alpha=0) matches solve_newmark(beta=0.25, gamma=0.5)

    • test_high_freq_dissipation — high-frequency component amplitude decays with alpha < 0

    • test_low_freq_preservation — low-frequency response unaffected by alpha

    • test_alpha_range_validation — error for alpha outside [-1/3, 0]

    • Priority: MEDIUM

12.7 tests/test_rayleigh_damping.py#

  • [ ] Tests:

    • test_rayleigh_formula — C = alphaM + betaK exactly

    • test_coefficients_at_target_frequencies — damping ratios match at omega1, omega2

    • test_sparse_dense_agree — sparse and dense versions give same C

    • test_zero_coefficients — alpha=0, beta=0 gives C=0

    • Priority: MEDIUM

12.8 tests/test_energy_balance.py#

  • [ ] Tests:

    • test_undamped_energy_constant — E_k + E_s = const (within 1e-8)

    • test_damped_energy_decreasing — E_total strictly decreasing for free vibration with damping

    • test_external_work_balance — W_ext = ΔE_k + ΔE_s + E_dissipated

    • test_energy_at_rest — E_k = 0, E_s = 0.5u^TK*u when velocity = 0

    • Priority: MEDIUM

Summary Statistics#

Phase

Tasks

Priority

New Files

1. Mass Matrices

9

HIGH

— (extend existing)

2. Damping Models

3

HIGH/MED/LOW

damping.py

3. Newmark-Beta

10

HIGH

dynamics.py

4. Central Difference

4

MEDIUM

(in dynamics.py)

5. HHT-Alpha

3

MEDIUM

(in dynamics.py)

6. Modal Analysis

6

HIGH

modal.py

7. Periodic BCs

10

HIGH

periodic.py

8. Nonlinear Dynamics

4

MEDIUM

(in dynamics.py)

9. Post-Processing

6

MEDIUM/LOW

(in plotting.py)

10. Infrastructure

6

HIGH

— (integration)

11. Examples

7

HIGH/MED

examples/dynamic_*.py, examples/periodic_*.py

12. Tests

8

HIGH/MED

tests/test_*.py

TOTAL

76 subtasks

4 new modules
On this page
Show Source