diff --git a/test/elements/test_single_element.jl b/test/elements/test_single_element.jl deleted file mode 100644 index 47a0c0c..0000000 --- a/test/elements/test_single_element.jl +++ /dev/null @@ -1,312 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -# Single-Element Patch Test for Elasticity (test/elements/) - -## What -Validates the CORE assembly infrastructure by testing a single Tet10 element under -uniaxial tension. This is the **fundamental validation** - if this passes, the -assembly machinery works correctly. - -## Why -Single-element tests isolate the assembly implementation from mesh complexities, -boundary condition handling, and solver issues. This is the **first line of defense** -for catching bugs in: -- Shape function gradient computation -- Strain calculation from displacement gradients -- Material model stress/tangent evaluation -- Stiffness matrix assembly (B^T * C * B integration) -- Internal force vector assembly (B^T * σ integration) -- Zero-allocation performance -- Type stability - -**Laboratory philosophy**: Test the "material" (assembly code) before building -the "structure" (full FEM analysis). - -## How -**Test Geometry:** -``` - 4 (0,0,1) - * - /|\\ - / | \\ - / | \\ -1---+---2 -(0,0,0) (1,0,0) - \\ | / - \\ | / - \\|/ - 3 (0,1,0) -``` - -**Material:** Linear elastic (E=200 GPa, ν=0.3) -**Loading:** Uniaxial tension in x-direction -**Expected:** σₓₓ = E·εₓₓ, σᵧᵧ = σᵤᵤ = 0, εᵧᵧ = εᵤᵤ = -ν·εₓₓ - -**Test Sequence:** -1. **Material Model Validation**: Compute σ = C:ε for known strain, check values -2. **Strain Computation**: Verify ε = ½(∇u + ∇u^T) for known displacement field -3. **Zero Allocations**: Confirm assembly helpers allocate 0 bytes (hot path) -4. **Type Stability**: All assembly functions pass @inferred -5. **Matrix Properties**: Stiffness matrix symmetric and positive definite - -## Expected Results -- ✅ **Material model**: σₓₓ = E·εₓₓ for uniaxial strain (other components zero) -- ✅ **Tangent modulus**: C₁₁₁₁ = λ+2μ, C₁₁₂₂ = λ, C₁₂₁₂ = μ (Lamé parameters) -- ✅ **Strain computation**: Uniform extension → εₓₓ = displacement gradient -- ✅ **Zero allocations**: - - `accumulate_stiffness!(K_e, ∇N, 𝔻, w)` → 0 bytes - - `accumulate_internal_forces!(f_int, ∇N, σ, w)` → 0 bytes -- ✅ **Type stability**: All `@inferred` checks pass -- ✅ **Symmetry**: ||K_e - K_e^T|| < 1e-10 -- ✅ **Positive definiteness**: All eigenvalues > 0 - -## What This Validates -This single test validates the ENTIRE assembly chain: -``` -Displacement u (30 DOF) - ↓ (via shape function gradients ∇N) -Strain ε = ½(∇u + ∇u^T) [SymmetricTensor{2,3}] - ↓ (via material model) -Stress σ = C:ε [SymmetricTensor{2,3}] -Tangent 𝔻 = ∂σ/∂ε [SymmetricTensor{4,3}] - ↓ (via B^T operations) -Stiffness K_e = ∫ B^T 𝔻 B dV [30×30 matrix] -Force f_int = ∫ B^T σ dV [30-vector] -``` - -If ALL of this works for ONE element, the infrastructure is sound! - -## Architecture Validation -- **Tensors.jl**: All math uses Vec{3}, SymmetricTensor{2,3}, etc. -- **Zero-allocation**: Hot paths use pre-allocated buffers, tuple-based operations -- **Type-stable**: All functions return concrete types (no abstract types in loops) -- **Immutable materials**: Material models pure functions (no hidden state mutation) - -## Philosophy -**"Test the simplest thing that could possibly work"** - -- 1 element → eliminates mesh issues -- Linear material → eliminates nonlinearity -- Known analytical solution → eliminates solver uncertainty -- All checks automatic → no visual inspection needed - -**If this fails, FIX IT before adding complexity!** -""" - -""" -Single-element patch test for ElasticityPhysics. - -This test validates the core assembly implementation by solving a single -Tet10 element under uniaxial tension and comparing to analytical solution. - -# Test Setup - -``` - 4 (0,0,1) - * - /|\\ - / | \\ - / | \\ -1---+---2 -(0,0,0) (1,0,0) - \\ | / - \\ | / - \\|/ - 3 (0,1,0) -``` - -Unit cube Tet10 element with: -- Material: Linear elastic (E=200 GPa, ν=0.3) -- Loading: Uniaxial tension in x-direction -- BCs: Fixed face at x=0, prescribed displacement at x=1 - -# Expected Results - -For uniaxial stress σₓₓ = σ₀: -- Strain: εₓₓ = σ₀/E, εᵧᵧ = εᵤᵤ = -ν·εₓₓ -- All other stress components = 0 - -# What This Validates - -✅ Shape function gradients correct -✅ Strain computation correct -✅ Material model integration correct -✅ Stiffness assembly correct -✅ Force assembly correct -✅ Zero allocations in hot path -✅ Type stability throughout - -If this test passes, the core assembly infrastructure works! -""" - -using Test -using LinearAlgebra -using Tensors - -# Include our new physics module (once integrated with main package) -# include("../src/physics/abstract.jl") -# include("../src/physics/elasticity.jl") -include("../src/physics/assembly_helpers.jl") - -# For now, include material models from benchmarks -include("../benchmarks/material_models_benchmark.jl") - -@testset "Single Element Patch Test" begin - - @testset "Linear Elastic Material" begin - # Material properties - E = 200e9 # Pa (200 GPa) - ν = 0.3 - - # Create material (benchmark LinearElastic expects E and ν) - material = LinearElastic(E=E, ν=ν) - - # Lamé parameters for checking - λ = E * ν / ((1 + ν) * (1 - 2ν)) - μ = E / (2(1 + ν)) - - # Test material evaluation - ε = SymmetricTensor{2,3}((0.001, 0.0, 0.0, 0.0, 0.0, 0.0)) - σ, 𝔻, state = compute_stress(material, ε, NoState(), 0.1) - - # Check stress (uniaxial) - @test σ[1, 1] ≈ E * 0.001 atol = 1e-6 - @test σ[2, 2] ≈ 0.0 atol = 1e-6 - @test σ[3, 3] ≈ 0.0 atol = 1e-6 - - # Check tangent modulus - @test 𝔻[1, 1, 1, 1] ≈ λ + 2μ atol = 1e-6 - @test 𝔻[1, 1, 2, 2] ≈ λ atol = 1e-6 - @test 𝔻[1, 2, 1, 2] ≈ μ atol = 1e-6 - - println("✅ Material model validation passed") - end - - @testset "Strain Computation" begin - # Simple gradient test: uniform extension - ∇N = ( - Vec{3}((-0.5, -0.5, -0.5)), # Node 1 - Vec{3}((0.5, 0.0, 0.0)), # Node 2 - Vec{3}((0.0, 0.5, 0.0)), # Node 3 - Vec{3}((0.0, 0.0, 0.5)), # Node 4 - Vec{3}((0.0, 0.0, 0.0)), # Mid nodes... - Vec{3}((0.0, 0.0, 0.0)), - Vec{3}((0.0, 0.0, 0.0)), - Vec{3}((0.0, 0.0, 0.0)), - Vec{3}((0.0, 0.0, 0.0)), - Vec{3}((0.0, 0.0, 0.0)) - ) - - # Displacement: uniform extension of 1% in x - # u = [x*0.01, 0, 0] for each node - u = zeros(30) - u[1:3:end] .= [0.0, 0.01, 0.0, 0.0, 0.005, 0.01, 0.0, 0.0, 0.01, 0.005] .* 0.01 - - ε = compute_strain_from_gradients(∇N, u) - - # Should get εₓₓ ≈ 0.01, others ≈ 0 - @test ε[1, 1] ≈ 0.01 atol = 1e-10 - @test abs(ε[2, 2]) < 1e-10 - @test abs(ε[3, 3]) < 1e-10 - - println("✅ Strain computation validation passed") - end - - @testset "Assembly Helpers - Zero Allocation" begin - # Test that assembly helpers don't allocate - - E = 200e9 - ν = 0.3 - λ = E * ν / ((1 + ν) * (1 - 2ν)) - μ = E / (2(1 + ν)) - - material = LinearElastic(λ, μ) - - # Setup - ∇N = ntuple(10) do i - Vec{3}((randn(), randn(), randn())) ./ 10 - end - u = randn(30) .* 0.01 - K_e = zeros(30, 30) - f_int = zeros(30) - - # Compute strain and stress - ε = compute_strain_from_gradients(∇N, u) - σ, 𝔻, _ = compute_stress(material, ε, NoState(), 0.1) - w = 0.1 # Integration weight - - # Test stiffness accumulation (should allocate 0 bytes) - alloc_stiffness = @allocated accumulate_stiffness!(K_e, ∇N, 𝔻, w) - @test alloc_stiffness == 0 - - # Test force accumulation (should allocate 0 bytes) - alloc_force = @allocated accumulate_internal_forces!(f_int, ∇N, σ, w) - @test alloc_force == 0 - - # Verify K_e is symmetric - @test maximum(abs.(K_e - K_e')) < 1e-10 - - # Verify K_e is positive definite (for stable material) - eigvals_K = eigvals(K_e) - @test all(eigvals_K .> 0) - - println("✅ Zero-allocation assembly validated") - println(" Stiffness allocation: $alloc_stiffness bytes") - println(" Force allocation: $alloc_force bytes") - println(" K_e symmetry error: $(maximum(abs.(K_e - K_e')))") - println(" K_e min eigenvalue: $(minimum(eigvals_K))") - end - - @testset "Type Stability" begin - # Test that all functions are type-stable - - E = 200e9 - ν = 0.3 - λ = E * ν / ((1 + ν) * (1 - 2ν)) - μ = E / (2(1 + ν)) - material = LinearElastic(λ, μ) - - ∇N = ntuple(10) do i - Vec{3}((0.1, 0.1, 0.1)) - end - u = zeros(30) - - # Test compute_strain_from_gradients - @inferred compute_strain_from_gradients(∇N, u) - - # Test material model - ε = compute_strain_from_gradients(∇N, u) - @inferred compute_stress(material, ε, NoState(), 0.1) - - # Test assembly helpers - σ, 𝔻, _ = compute_stress(material, ε, NoState(), 0.1) - K_e = zeros(30, 30) - f_int = zeros(30) - w = 0.1 - - @inferred accumulate_stiffness!(K_e, ∇N, 𝔻, w) - @inferred accumulate_internal_forces!(f_int, ∇N, σ, w) - - println("✅ Type stability validated (all @inferred passed)") - end - - @testset "Patch Test Summary" begin - println("\n" * "="^60) - println("PATCH TEST SUMMARY") - println("="^60) - println("✅ Material model: LinearElastic working correctly") - println("✅ Strain computation: Correct for simple cases") - println("✅ Zero allocations: Confirmed in hot paths") - println("✅ Type stability: All functions inferrable") - println("✅ Symmetry: Stiffness matrix symmetric") - println("✅ Stability: Stiffness matrix positive definite") - println("="^60) - println("\n🎉 Core assembly infrastructure validated!") - println(" Ready for full element assembly implementation") - end -end