mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-19 09:54:55 +00:00
chore(test): delete modal analysis new-API experiment
Remove unmaintained modal solver API harness. - Drop `test/solvers/test_modal_analysis_new_api.jl`.
This commit is contained in:
Jukka Aho
@@ -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)
|
||||
|
||||
"""
|
||||
Reference in New Issue
Block a user