diff --git a/test/geometry/test_jacobian.jl b/test/geometry/test_jacobian.jl index 89d41ef..b5e2059 100644 --- a/test/geometry/test_jacobian.jl +++ b/test/geometry/test_jacobian.jl @@ -1,384 +1,42 @@ -""" -# Jacobian Computation Tests (test/geometry/) - -## What -Tests Jacobian matrix computation for mapping between parametric (reference) and -physical (mesh) coordinates. The Jacobian **J = ∂x/∂ξ** relates shape function -derivatives in parametric space to physical space. - -## Why -The Jacobian is **THE MOST CRITICAL** geometric computation in FEM: -- **Coordinate transformation**: Maps ∂/∂ξ → ∂/∂x via J⁻¹ -- **Differential volume**: det(J) provides dV = det(J) dξ for integration -- **Element quality**: det(J) > 0 required (negative → inverted element) -- **Strain computation**: ε = f(∇u), and ∇ requires Jacobian transformation - -Without correct Jacobian: -- Stiffness matrices are wrong -- Internal forces are wrong -- EVERYTHING is wrong! - -This test validates: -- **Correctness**: Known geometric transformations produce expected J -- **Orientation**: det(J) > 0 for well-shaped elements -- **Quality detection**: det(J) ≈ 0 for degenerate elements -- **Physical derivatives**: dN/dx = J⁻¹ · dN/dξ correct -- **Type stability**: Returns Tensor{2,D} (zero-allocation) -- **Zero allocations**: Hot path allocates nothing - -## How -**Test Cases:** - -**1. Identity Mapping** (reference → reference): -- Triangle: Vertices at (0,0), (1,0), (0,1) -- Expected: J = I (identity matrix), det(J) = 1.0 - -**2. Scaled Elements**: -- Triangle scaled 2× in x, 1.5× in y -- Expected: J = diag(2.0, 1.5), det(J) = 3.0 (area scaling) -- Tetrahedron scaled 2×, 3×, 4× -- Expected: det(J) = 24.0 (volume scaling) - -**3. Rotated Elements**: -- 90° rotation of triangle -- Expected: det(J) = 1.0 (area preserved), J contains rotation matrix - -**4. Physical Derivatives**: -- Computes dN/dx = J⁻¹ · dN/dξ -- Validates constant strain condition: ∑ᵢ dNᵢ/dx = 0 (partition of unity derivative) - -**5. Element Quality**: -- **Well-shaped**: det(J) > 0.1 (properly oriented, well-conditioned) -- **Degenerate**: det(J) < 1e-10 (collapsed to line/plane, unusable) - -**6. Type Stability**: -- `@inferred compute_jacobian(X, dN_dξ)` → Tensor{2,D} -- `@inferred physical_derivatives(J, dN_dξ)` → Tuple -- `@allocated` checks confirm zero allocations - -**7. Manual Verification**: -- Triangle with vertices (1,2), (4,3), (2,6) -- Hand-calculated J = [3 1; 1 4], det(J) = 11 -- Verifies implementation matches theory - -**8. Vector vs Tuple Interface**: -- Tests both NTuple{N,Vec{D}} (preferred) and Vector{Vec{D}} (legacy) -- Both produce same results (tuple faster due to stack allocation) - -## Expected Results -- ✅ **Identity mapping**: J = I, det(J) = 1.0 -- ✅ **Scaled mapping**: J diagonal with scale factors, det(J) = product of scales -- ✅ **Rotated mapping**: det(J) = 1.0 (area/volume preserved) -- ✅ **Physical derivatives**: ∑ᵢ dNᵢ/dx = 0 (constant strain condition) -- ✅ **Quality detection**: det(J) > 0 for valid, ≈ 0 for degenerate -- ✅ **Type stability**: All @inferred checks pass -- ✅ **Zero allocations**: All @allocated checks return 0 -- ✅ **Consistency**: Matches hand calculations - -## Mathematical Background -**Jacobian Matrix:** -``` -J_ij = ∂xᵢ/∂ξⱼ = ∑ₖ Xₖ,ᵢ · dNₖ/dξⱼ -``` - -Where: -- Xₖ = physical coordinates of node k (Vec{D}) -- dNₖ/dξⱼ = parametric derivative of shape function k w.r.t. ξⱼ -- D = physical dimension (2 or 3) -- d = parametric dimension (1, 2, or 3) - -**Physical Derivatives:** -``` -∂Nₖ/∂x = J⁻¹ · ∂Nₖ/∂ξ -``` - -**Integration:** -``` -∫_Ω f(x) dV = ∫_Ω_ref f(ξ) |det(J)| dξ -``` - -## Architecture Principle -**Tensors.jl for ALL Math** - -All geometric quantities use Tensors.jl types: -- Coordinates: Vec{D,Float64} -- Jacobian: Tensor{2,D} (or Tensor{2,D,Float64,M} for D×d) -- Derivatives: Vec{D,Float64} - -This provides: -- Natural mathematical notation (J[i,j]) -- Automatic differentiation compatible -- Zero-allocation operations -- GPU transferable - -## Critical for Assembly -Every integration point requires: -1. Evaluate dN/dξ at quadrature point -2. Compute J = X · dN/dξ -3. Compute dN/dx = J⁻¹ · dN/dξ -4. Compute det(J) for integration weight - -If ANY of these allocate, assembly is slow. These tests ensure they don't! -""" +# SPDX-FileCopyrightText: 2015-2026 Jukka Aho +# SPDX-License-Identifier: MIT using Test using JuliaFEM using Tensors -using LinearAlgebra -@testset "Jacobian Computation" begin - - @testset "2D Triangle - Identity Element" begin - # Reference triangle mapped to itself (identity transformation) - X = ( - Vec{2}((0.0, 0.0)), - Vec{2}((1.0, 0.0)), - Vec{2}((0.0, 1.0)) - ) - - # Evaluate at center +@testset "Jacobian helpers" begin + @testset "compute_jacobian tuple (triangle)" begin + X = (Vec{2}((0.0, 0.0)), Vec{2}((2.0, 0.0)), Vec{2}((0.0, 1.5))) xi = Vec{2}((1 / 3, 1 / 3)) - dN_dξ = get_basis_derivatives(Triangle(), Lagrange{Triangle,1}(), xi) - - # Compute Jacobian - J = compute_jacobian(X, dN_dξ) - - # For identity mapping, J should be identity matrix - @test J ≈ Tensor{2,2}((1.0, 0.0, 0.0, 1.0)) - @test det(J) ≈ 1.0 + dN = Tuple(get_basis_derivatives(Triangle{3}(), Lagrange{1}(), xi)) + J = compute_jacobian(X, dN) + @test J ≈ Tensor{2,2}((2.0, 0.0, 0.0, 1.5)) + @test det(J) ≈ 3.0 end - @testset "2D Triangle - Scaled Element" begin - # Triangle scaled by 2 in x and 1.5 in y - X = ( - Vec{2}((0.0, 0.0)), - Vec{2}((2.0, 0.0)), - Vec{2}((0.0, 1.5)) - ) + @testset "compute_jacobian AbstractVector" begin + Xv = [Vec{2}((0.0, 0.0)), Vec{2}((1.0, 0.0)), Vec{2}((0.0, 1.0))] + xi = Vec{2}((0.2, 0.2)) + dNsv = get_basis_derivatives(Triangle{3}(), Lagrange{1}(), xi) + dNv = collect(dNsv) + J1 = compute_jacobian(Xv, dNv) + J2 = compute_jacobian((Xv...,), Tuple(dNsv)) + @test isapprox(J1, J2; rtol=1e-14) + end + @testset "physical_derivatives" begin + X = (Vec{2}((0.0, 0.0)), Vec{2}((2.0, 0.0)), Vec{2}((0.0, 1.5))) xi = Vec{2}((1 / 3, 1 / 3)) - dN_dξ = get_basis_derivatives(Triangle(), Lagrange{Triangle,1}(), xi) - - J = compute_jacobian(X, dN_dξ) - - # Jacobian should reflect scaling - @test J[1, 1] ≈ 2.0 # ∂x/∂ξ - @test J[1, 2] ≈ 0.0 # ∂x/∂η - @test J[2, 1] ≈ 0.0 # ∂y/∂ξ - @test J[2, 2] ≈ 1.5 # ∂y/∂η - @test det(J) ≈ 3.0 # Area scaling = 2 × 1.5 - end - - @testset "2D Triangle - Rotated Element" begin - # 90° counter-clockwise rotation - θ = π / 2 - R = [cos(θ) -sin(θ); sin(θ) cos(θ)] - - # Original nodes - X_orig = [0.0 1.0 0.0; 0.0 0.0 1.0] - - # Rotate - X_rot = R * X_orig - X = ( - Vec{2}((X_rot[1, 1], X_rot[2, 1])), - Vec{2}((X_rot[1, 2], X_rot[2, 2])), - Vec{2}((X_rot[1, 3], X_rot[2, 3])) - ) - - xi = Vec{2}((1 / 3, 1 / 3)) - dN_dξ = get_basis_derivatives(Triangle(), Lagrange{Triangle,1}(), xi) - - J = compute_jacobian(X, dN_dξ) - - # Jacobian should contain rotation - @test det(J) ≈ 1.0 # Area preserved under rotation - @test norm(J) > 0 # Well-conditioned - end - - @testset "3D Tetrahedron - Identity Element" begin - # Reference tetrahedron mapped to itself - X = ( - Vec{3}((0.0, 0.0, 0.0)), - Vec{3}((1.0, 0.0, 0.0)), - Vec{3}((0.0, 1.0, 0.0)), - Vec{3}((0.0, 0.0, 1.0)) - ) - - xi = Vec{3}((0.25, 0.25, 0.25)) - dN_dξ = get_basis_derivatives(Tetrahedron(), Lagrange{Tetrahedron,1}(), xi) - - J = compute_jacobian(X, dN_dξ) - - # Identity mapping - @test J ≈ Tensor{2,3}((1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) - @test det(J) ≈ 1.0 - end - - @testset "3D Tetrahedron - Scaled Element" begin - # Tetrahedron scaled differently in each direction - X = ( - Vec{3}((0.0, 0.0, 0.0)), - Vec{3}((2.0, 0.0, 0.0)), - Vec{3}((0.0, 3.0, 0.0)), - Vec{3}((0.0, 0.0, 4.0)) - ) - - xi = Vec{3}((0.25, 0.25, 0.25)) - dN_dξ = get_basis_derivatives(Tetrahedron(), Lagrange{Tetrahedron,1}(), xi) - - J = compute_jacobian(X, dN_dξ) - - # Diagonal Jacobian (aligned with axes) - @test J[1, 1] ≈ 2.0 - @test J[2, 2] ≈ 3.0 - @test J[3, 3] ≈ 4.0 - @test det(J) ≈ 24.0 # Volume scaling = 2 × 3 × 4 - end - - @testset "Physical Derivatives - 2D Triangle" begin - X = ( - Vec{2}((0.0, 0.0)), - Vec{2}((2.0, 0.0)), - Vec{2}((0.0, 1.5)) - ) - - xi = Vec{2}((1 / 3, 1 / 3)) - dN_dξ = get_basis_derivatives(Triangle(), Lagrange{Triangle,1}(), xi) - - J = compute_jacobian(X, dN_dξ) - dN_dx = physical_derivatives(J, dN_dξ) - - # Verify constant strain condition: ∑ᵢ dNᵢ/dx = 0 - sum_dN_dx = sum(dN_dx) - @test norm(sum_dN_dx) < 1e-10 - - # Verify partition of unity holds - # (Not directly, but derivatives should be consistent) - @test length(dN_dx) == 3 - end - - @testset "Physical Derivatives - 3D Tetrahedron" begin - X = ( - Vec{3}((0.0, 0.0, 0.0)), - Vec{3}((1.0, 0.0, 0.0)), - Vec{3}((0.0, 1.0, 0.0)), - Vec{3}((0.0, 0.0, 1.0)) - ) - - xi = Vec{3}((0.25, 0.25, 0.25)) - dN_dξ = get_basis_derivatives(Tetrahedron(), Lagrange{Tetrahedron,1}(), xi) - - J = compute_jacobian(X, dN_dξ) - dN_dx = physical_derivatives(J, dN_dξ) - - # Constant strain condition - sum_dN_dx = sum(dN_dx) - @test norm(sum_dN_dx) < 1e-10 - - # Check each derivative is a 3D vector - for dN in dN_dx - @test length(dN) == 3 + dN = Tuple(get_basis_derivatives(Triangle{3}(), Lagrange{1}(), xi)) + J = compute_jacobian(X, dN) + dNdx_t = physical_derivatives(J, dN) + dNdx_v = physical_derivatives(J, collect(dN)) + @test length(dNdx_t) == 3 + @test length(dNdx_v) == 3 + for i in 1:3 + @test dNdx_t[i] ≈ dNdx_v[i] end - end - - @testset "Jacobian Determinant - Element Quality" begin - # Well-shaped triangle - X_good = ( - Vec{2}((0.0, 0.0)), - Vec{2}((1.0, 0.0)), - Vec{2}((0.0, 1.0)) - ) - - xi = Vec{2}((1 / 3, 1 / 3)) - dN_dξ = get_basis_derivatives(Triangle(), Lagrange{Triangle,1}(), xi) - J_good = compute_jacobian(X_good, dN_dξ) - - @test det(J_good) > 0 # Positive (properly oriented) - @test abs(det(J_good)) > 0.1 # Well-conditioned - - # Degenerate triangle (collapsed to line) - X_bad = ( - Vec{2}((0.0, 0.0)), - Vec{2}((1.0, 0.0)), - Vec{2}((2.0, 0.0)) # Collinear! - ) - - J_bad = compute_jacobian(X_bad, dN_dξ) - @test abs(det(J_bad)) < 1e-10 # Nearly zero (degenerate) - end - - @testset "Type Stability and Zero Allocation" begin - X = ( - Vec{2}((0.0, 0.0)), - Vec{2}((1.0, 0.0)), - Vec{2}((0.0, 1.0)) - ) - - xi = Vec{2}((1 / 3, 1 / 3)) - dN_dξ = get_basis_derivatives(Triangle(), Lagrange{Triangle,1}(), xi) - - # Type stability - J = @inferred compute_jacobian(X, dN_dξ) - @test J isa Tensor{2,2} - - dN_dx = @inferred physical_derivatives(J, dN_dξ) - @test dN_dx isa Tuple - - # Zero allocation (run twice to avoid compilation) - compute_jacobian(X, dN_dξ) - allocs = @allocated compute_jacobian(X, dN_dξ) - @test allocs == 0 - - physical_derivatives(J, dN_dξ) - allocs = @allocated physical_derivatives(J, dN_dξ) - @test allocs == 0 - end - - @testset "Consistency with Manual Calculation" begin - # Triangle with known Jacobian - X = ( - Vec{2}((1.0, 2.0)), - Vec{2}((4.0, 3.0)), - Vec{2}((2.0, 6.0)) - ) - - xi = Vec{2}((0.5, 0.25)) - dN_dξ = get_basis_derivatives(Triangle(), Lagrange{Triangle,1}(), xi) - # dN_dξ = (Vec(-1, -1), Vec(1, 0), Vec(0, 1)) - - J = compute_jacobian(X, dN_dξ) - - # Manual calculation: - # J = X2 - X1 in first column, X3 - X1 in second column - # J = [4-1 2-1] = [3 1] - # [3-2 6-2] [1 4] - - @test J[1, 1] ≈ 3.0 - @test J[1, 2] ≈ 1.0 - @test J[2, 1] ≈ 1.0 - @test J[2, 2] ≈ 4.0 - @test det(J) ≈ 11.0 # 3*4 - 1*1 = 11 + @test sum(dNdx_t) ≈ Vec{2}((0.0, 0.0)) end end - -@testset "Jacobian - AbstractVector Interface" begin - # Test that Vector interface also works (less efficient) - X_vec = [Vec{2}((0.0, 0.0)), Vec{2}((1.0, 0.0)), Vec{2}((0.0, 1.0))] - - xi = Vec{2}((1 / 3, 1 / 3)) - dN_dξ_tuple = get_basis_derivatives(Triangle(), Lagrange{Triangle,1}(), xi) - dN_dξ_vec = collect(dN_dξ_tuple) - - J_tuple = compute_jacobian(tuple(X_vec...), dN_dξ_tuple) - J_vec = compute_jacobian(X_vec, dN_dξ_vec) - - @test J_tuple ≈ J_vec - - # Physical derivatives - dN_dx_tuple = physical_derivatives(J_tuple, dN_dξ_tuple) - dN_dx_vec = physical_derivatives(J_vec, dN_dξ_vec) - - @test all(dN_dx_tuple[i] ≈ dN_dx_vec[i] for i in 1:3) -end - -println("✅ All Jacobian tests passed!")