From fb589d27c40c09dbf598579152da72940bdae137 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 9 May 2026 18:30:06 +0300 Subject: [PATCH] chore(test): delete modal analysis new-API experiment Remove unmaintained modal solver API harness. - Drop `test/solvers/test_modal_analysis_new_api.jl`. --- test/solvers/test_modal_analysis_new_api.jl | 560 -------------------- 1 file changed, 560 deletions(-) delete mode 100644 test/solvers/test_modal_analysis_new_api.jl diff --git a/test/solvers/test_modal_analysis_new_api.jl b/test/solvers/test_modal_analysis_new_api.jl deleted file mode 100644 index d7ba957..0000000 --- a/test/solvers/test_modal_analysis_new_api.jl +++ /dev/null @@ -1,560 +0,0 @@ -""" -# Modal Analysis - NEW API (Test-Driven Development) - -**What:** Shows how modal (eigenvalue) analysis SHOULD work with the NEW API - -**Why:** -- **Eigenvalue problems** - Different solver type (not Newton-Krylov) -- **Natural frequencies** - Validates dynamic behavior -- **Geometric stiffness** - Prestress effects on frequencies -- **Multiple physics** - Works for elasticity AND heat transfer -- **Mass matrix required** - First time we need M, not just K - -**NEW API Concepts:** -1. **ModalSolver** - Eigenvalue solver (not LinearSolver or NewtonSolver) -2. **EigenvalueProblem** - K φ = λ M φ (generalized eigenvalue problem) -3. **Geometric stiffness** - K_g from initial stress/displacement -4. **Modal shapes** - φ_i (eigenvectors = mode shapes) -5. **Natural frequencies** - ω_i = √λ_i - -**Test Problems:** - -## Test 1: Single Tet4 Element -- 4-node tetrahedron with 3 nodes fixed -- Validates basic eigenvalue computation -- Tests with and without geometric stiffness - -## Test 2: Heat Transfer Modal Analysis -- Two Quad4 elements (plane heat) -- Eigenvalues of thermal diffusion operator -- Validates modal analysis works for temperature field - -## Test 3: Vibrating Beam (Future) -- Fixed-fixed beam -- Compare natural frequencies to analytical solution -- ω_i = λ_i² √(EI/ρA) - -**Expected Behavior (when implemented):** -✅ Eigenvalues computed correctly -✅ Mode shapes orthogonal (φ_i^T M φ_j = δ_ij) -✅ Geometric stiffness affects frequencies -✅ Works for both elasticity and heat transfer -✅ Sparse eigenvalue solver (only need first k modes) - -**Status:** 🚧 VISIONARY TEST - Implementation in progress -""" - -using Test -using JuliaFEM -using Tensors -using LinearAlgebra -using Statistics - -@testset "Modal Analysis - NEW API (TDD)" begin - - # ============================================================================= - # SINGLE ELEMENT: TET4 WITHOUT GEOMETRIC STIFFNESS - # ============================================================================= - - @testset "Single Tet4 Modal Analysis (Visionary)" begin - @test_skip begin # Skip until implemented - - # Geometry - X = Dict( - 1 => Vec(2.0, 3.0, 4.0), - 2 => Vec(6.0, 3.0, 2.0), - 3 => Vec(2.0, 5.0, 1.0), - 4 => Vec(4.0, 3.0, 6.0) - ) - - # Initial displacement (for geometric stiffness later) - u0 = Dict( - 1 => Vec(0.0, 0.0, 0.0), - 2 => Vec(0.0, 0.0, 0.0), - 3 => Vec(0.0, 0.0, 0.0), - 4 => Vec(0.25, 0.25, 0.25) - ) - - # Material - material = LinearElastic(E=96.0, ν=1 / 3, ρ=420.0) - - # Physics - elastic_physics = ContinuumPhysics{Displacement}( - material=material, - formulation=FullThreeD(), - finite_strain=false, - geometric_stiffness=false # No K_g initially - ) - - # Create element - tet_element = create_element( - Tet4, (1, 2, 3, 4), - geometry=X, - displacement=u0 - ) - - # Domain - domain = Domain( - name="TET", - elements=[tet_element], - physics=elastic_physics - ) - - # Boundary conditions (3 nodes fixed) - bc_fixed = DirichletBC( - name="FIXED_FACE", - nodes=[1, 2, 3], - dof=:displacement, - values=[0.0, 0.0, 0.0] - ) - - # NEW: EigenvalueProblem (not static or transient) - problem = EigenvalueProblem( - domains=[domain], - boundary_conditions=[bc_fixed] - ) - - # NEW: ModalSolver - solver = ModalSolver( - n_modes=2, # Number of modes to compute - which=:SM, # Smallest Magnitude (or :LM for largest) - method=:arpack # ARPACK, KrylovKit, or direct - ) - - # Solve: K φ = λ M φ - solution = solve!(problem, solver) - - # Extract eigenvalues and eigenvectors - λ = solution.eigenvalues # [λ_1, λ_2] - φ = solution.eigenvectors # [φ_1, φ_2] - - # Validate eigenvalues (without geometric stiffness) - @test isapprox(λ, [4 / 3, 1 / 3], rtol=1e-6) - - # Validate orthogonality: φ_i^T M φ_j = δ_ij - for i in 1:2 - for j in 1:2 - orthogonality = φ[i]' * solution.mass_matrix * φ[j] - expected = (i == j) ? 1.0 : 0.0 - @test isapprox(orthogonality, expected, atol=1e-6) - end - end - - end - end - - # ============================================================================= - # GEOMETRIC STIFFNESS EFFECT - # ============================================================================= - - @testset "Tet4 with Geometric Stiffness (Visionary)" begin - @test_skip begin - - # Same setup as before, but enable geometric stiffness - elastic_physics = ContinuumPhysics{Displacement}( - material=LinearElastic(E=96.0, ν=1 / 3, ρ=420.0), - formulation=FullThreeD(), - finite_strain=false, - geometric_stiffness=true # NEW: Enable K_g - ) - - problem = EigenvalueProblem( - domains=[domain], - boundary_conditions=[bc_fixed], - initial_displacement=u0 # Required for K_g - ) - - solver = ModalSolver(n_modes=2, which=:SM) - solution = solve!(problem, solver) - - λ = solution.eigenvalues - - # With geometric stiffness: K_eff = K + K_g(u0) - # Eigenvalues should be different! - @test isapprox(λ, [5 / 3, 2 / 3], rtol=1e-6) - - # Geometric stiffness changes natural frequencies - @test λ[1] > 4 / 3 # Stiffened by prestress - @test λ[2] > 1 / 3 - - end - end - - # ============================================================================= - # HEAT TRANSFER MODAL ANALYSIS - # ============================================================================= - - @testset "Heat Transfer Modal (Visionary)" begin - @test_skip begin - - # Geometry: Two quad elements stacked vertically - X = Dict( - 1 => Vec(0.0, 0.0), - 2 => Vec(1.0, 0.0), - 3 => Vec(1.0, 3.0), - 4 => Vec(0.0, 3.0), - 5 => Vec(0.0, 3.0), - 6 => Vec(1.0, 3.0), - 7 => Vec(1.0, 9.0), - 8 => Vec(0.0, 9.0) - ) - - # Material (thermal) - thermal_material = ThermalMaterial( - conductivity=36.0, - density=6.0, - specific_heat=1.0 - ) - - # Physics - thermal_physics = ContinuumPhysics{Temperature}( - material=thermal_material, - formulation=PlaneHeat() - ) - - # Elements - elem1 = create_element(Quad4, (1, 2, 3, 4), geometry=X) - elem2 = create_element(Quad4, (4, 3, 7, 8), geometry=X) - - # Domain - domain = Domain( - name="THERMAL", - elements=[elem1, elem2], - physics=thermal_physics - ) - - # Boundary conditions (fixed temperature at ends) - bc_bottom = DirichletBC( - nodes=[1, 2], - dof=:temperature, - value=0.0 - ) - - bc_top = DirichletBC( - nodes=[7, 8], - dof=:temperature, - value=0.0 - ) - - # Eigenvalue problem for heat equation - # M ∂T/∂t + K T = 0 - # Modal: K φ = λ M φ - problem = EigenvalueProblem( - domains=[domain], - boundary_conditions=[bc_bottom, bc_top] - ) - - solver = ModalSolver(n_modes=1, which=:SM) - solution = solve!(problem, solver) - - λ = solution.eigenvalues - - # First eigenvalue should be 1.0 (analytical) - @test isapprox(λ[1], 1.0, rtol=1e-6) - - end - end - - # ============================================================================= - # NATURAL FREQUENCIES FROM EIGENVALUES - # ============================================================================= - - @testset "Natural Frequency Computation (Visionary)" begin - @test_skip begin - - # Solve modal problem - solution = solve_modal_problem() - - # Extract eigenvalues - λ = solution.eigenvalues # [rad²/s²] - - # Natural frequencies: ω = √λ [rad/s] - ω = sqrt.(λ) - - # Convert to Hz: f = ω / (2π) - f = ω ./ (2π) - - # Validate units and values - @test all(λ .≥ 0) # Eigenvalues non-negative - @test all(ω .≥ 0) # Frequencies non-negative - @test all(f .≥ 0) - - # For simple beam: f_1 ≈ 1-10 Hz (typical) - @test 1.0 < f[1] < 100.0 - - end - end - - # ============================================================================= - # MODE SHAPE VISUALIZATION (FUTURE) - # ============================================================================= - - @testset "Mode Shape Properties (Visionary)" begin - @test_skip begin - - solution = solve_modal_problem() - - φ = solution.eigenvectors # Mode shapes - M = solution.mass_matrix - K = solution.stiffness_matrix - λ = solution.eigenvalues - - # Property 1: Orthogonality w.r.t. mass matrix - # φ_i^T M φ_j = δ_ij - for i in 1:length(φ) - for j in 1:length(φ) - orth_M = φ[i]' * M * φ[j] - expected = (i == j) ? 1.0 : 0.0 - @test isapprox(orth_M, expected, atol=1e-6) - end - end - - # Property 2: Eigenvalue equation - # K φ_i = λ_i M φ_i - for i in 1:length(φ) - Kφ = K * φ[i] - λMφ = λ[i] * (M * φ[i]) - @test Kφ ≈ λMφ rtol = 1e-6 - end - - # Property 3: Mode shapes normalized - # φ_i^T M φ_i = 1 - for i in 1:length(φ) - norm_M = φ[i]' * M * φ[i] - @test isapprox(norm_M, 1.0, rtol=1e-6) - end - - end - end - - # ============================================================================= - # MASS MATRIX ASSEMBLY - # ============================================================================= - - @testset "Mass Matrix Assembly (Visionary)" begin - # Pseudo-code showing mass matrix assembly - - println("\n" * "="^70) - println("MASS MATRIX ASSEMBLY (NODAL APPROACH)") - println("="^70) - - mass_assembly_pseudo = """ - # Similar to stiffness, but uses density ρ and N^T N - for node_i in nodes - for elem in node_to_elements[node_i] - for node_j in elem.nodes - # Mass matrix block (3×3 for displacement) - M_ij = ∫_Ω ρ N_i N_j dΩ - - # For displacement (3D): - # M_ij = m_ij * I_3×3 (often lumped) - - # Consistent mass (full integration) - M_block = compute_mass_block(elem, node_i, node_j, ρ) - - # Or lumped mass (diagonal only) - if lumped - M_block = (i == j) ? diag(M_block) : 0 - end - - M_nodal[node_i, node_j] += M_block - end - end - end - """ - - println(mass_assembly_pseudo) - println("="^70) - println("✓ Mass matrix M = ∫_Ω ρ N^T N dΩ") - println("✓ Consistent mass: Full integration (more accurate)") - println("✓ Lumped mass: Diagonal only (faster, explicit dynamics)") - println("✓ Nodal assembly works same as stiffness!") - println("="^70) - end - - # ============================================================================= - # KEY ARCHITECTURAL INSIGHTS - # ============================================================================= - - println("\n" * "="^70) - println("MODAL ANALYSIS ARCHITECTURE INSIGHTS (NEW API)") - println("="^70) - println("✓ Modal analysis is EigenvalueProblem (K φ = λ M φ)") - println("✓ Requires both stiffness K AND mass M matrices") - println("✓ ModalSolver uses sparse eigenvalue methods (ARPACK)") - println("✓ Geometric stiffness K_g affects natural frequencies") - println("✓ Works for ANY physics (displacement, temperature, etc.)") - println("✓ Mode shapes φ orthogonal w.r.t. mass matrix") - println("✓ Natural frequencies ω = √λ, in Hz: f = ω/(2π)") - println("✓ Same nodal assembly pattern for M as for K!") - println("="^70) - -end - -""" -# IMPLEMENTATION NOTES - -## Generalized Eigenvalue Problem - -**Mathematical formulation:** - -K φ = λ M φ - -where: -- K = stiffness matrix (from elastic/thermal energy) -- M = mass matrix (from kinetic/thermal energy) -- φ = eigenvector (mode shape) -- λ = eigenvalue (ω² for vibrations, α for heat) - -**For structural vibrations:** -- ω = √λ (natural frequency in rad/s) -- f = ω/(2π) (natural frequency in Hz) -- φ = displacement mode shape - -**For heat diffusion:** -- α = λ (thermal diffusivity eigenvalue) -- φ = temperature mode shape - -## Stiffness Matrix - -**Elasticity:** - -K = ∫_Ω B^T C B dΩ - -where: -- B = strain-displacement matrix -- C = elasticity tensor - -**Heat transfer:** - -K = ∫_Ω k ∇N^T ∇N dΩ - -where: -- k = thermal conductivity - -## Mass Matrix - -**Elasticity (consistent):** - -M = ∫_Ω ρ N^T N dΩ - -where ρ = density. - -**Elasticity (lumped):** - -M_ii = ∑_elements ∫_Ω_e ρ N_i dΩ_e - -(Diagonal only, faster for explicit dynamics) - -**Heat transfer:** - -M = ∫_Ω ρc N^T N dΩ - -where: -- ρ = density -- c = specific heat - -## Geometric Stiffness - -**Definition:** Stiffness contribution from initial stress state. - -K_g = ∫_Ω G^T σ₀ G dΩ - -where: -- G = geometric matrix (relates δε to ∇(δu)) -- σ₀ = initial stress - -**Effect:** Prestress changes natural frequencies: -- Tension → higher frequencies (stiffening) -- Compression → lower frequencies (softening) - -**Total effective stiffness:** - -K_eff = K + K_g - -Then solve: K_eff φ = λ M φ - -## Sparse Eigenvalue Solvers - -**Problem:** For large systems, computing ALL eigenvalues is expensive. - -**Solution:** Sparse eigenvalue methods (compute only k << n eigenvalues). - -**Methods:** -1. **ARPACK** - Arnoldi iteration (standard in Julia) -2. **KrylovKit.jl** - Modern Krylov methods -3. **Lanczos** - For symmetric problems (K and M symmetric) - -**JuliaFEM approach:** - -```julia -using Arpack - -# Solve K φ = λ M φ for k smallest eigenvalues -λ, φ = eigs(K, M, nev=k, which=:SM) -``` - -## Mode Shape Normalization - -**Goal:** Normalize eigenvectors for convenience. - -**Mass normalization (standard):** - -φ_i^T M φ_i = 1 - -**Advantages:** -- Orthogonality: φ_i^T M φ_j = δ_ij -- Energy interpretation clear -- Modal damping ratios well-defined - -**Implementation:** - -```julia -function normalize_modes!(φ, M) - for i in 1:length(φ) - # Compute φ_i^T M φ_i - norm_M = φ[i]' * M * φ[i] - - # Normalize - φ[i] ./= sqrt(norm_M) - end -end -``` - -## Nodal Assembly for Mass Matrix - -```julia -function assemble_mass_nodal!(M_nodal, ρ, elements, node_to_elements) - Threads.@threads for node_i in 1:n_nodes - for elem_idx in node_to_elements[node_i] - elem = elements[elem_idx] - - for node_j in elem.nodes - # Mass matrix block (3×3 for displacement) - M_ij = compute_mass_block(elem, node_i, node_j, ρ) - - # Add to global (no atomic on diagonal if i = j) - if node_i == node_j - M_nodal[node_i, node_i] += M_ij # Direct write - else - atomic_add!(M_nodal[node_i, node_j], M_ij) - end - end - end - end -end -``` - -**Key:** SAME PATTERN as stiffness assembly! - -## Next Steps - -1. Implement `EigenvalueProblem` type -2. Implement `ModalSolver` with ARPACK integration -3. Implement mass matrix assembly (nodal) -4. Implement geometric stiffness computation -5. Add mode shape normalization -6. Validate against analytical solutions -7. Performance benchmarks (large systems) - -"""