Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-25 11:35:30 +00:00
Code Issues Packages Projects Releases Wiki Activity

test(geometry): slim Jacobian suite to Triangle{3} Lagrange{1} smoke tests

Browse Source 
The previous file duplicated a long narrative and many scenarios that better live
in basis/quadrature coverage; keep the Jacobian helpers exercised with tuple vs
vector coordinates and physical derivative consistency.

- SPDX header; drop unused `LinearAlgebra` import.
- Fix basis calls to `Triangle{3}` / `Lagrange{1}` and verify scaling + PoU gradient sum.
This commit is contained in:
Jukka Aho
2026-05-09 18:46:03 +03:00
parent cc71c10033
commit 00e69548b9
1 changed files with 29 additions and 371 deletions
Download Patch File Download Diff File Expand all files Collapse all files
test/geometry/test_jacobian.jl
+29 -371
View File
@@ -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!")