From 98b126563f25e4054863148579ab5145dd2f3415 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 15 Dec 2025 08:16:26 +0200 Subject: [PATCH] test(geometry): add Jacobian computation test MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New 384-line test file for Jacobian matrix computation: - Tests J = ∂x/∂ξ mapping between parametric and physical coordinates - Tests identity, scaled, and rotated element transformations - Tests 2D triangles and 3D tetrahedra - Tests physical derivatives: dN/dx = J⁻¹ · dN/dξ - Validates element quality detection (det(J) > 0) - Tests consistency with manual calculations - Validates type stability and zero-allocation performance - Tests both NTuple and Vector interfaces Comprehensive test for Jacobian computation essential for coordinate transformation and integration in FEM assembly. --- test/geometry/test_jacobian.jl | 384 +++++++++++++++++++++++++++++++++ 1 file changed, 384 insertions(+) create mode 100644 test/geometry/test_jacobian.jl diff --git a/test/geometry/test_jacobian.jl b/test/geometry/test_jacobian.jl new file mode 100644 index 0000000..89d41ef --- /dev/null +++ b/test/geometry/test_jacobian.jl @@ -0,0 +1,384 @@ +""" +# 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! +""" + +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 + 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 + 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)) + ) + + 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 + 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 + 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!")