diff --git a/test/materials/test_hyperelasticity_new_api.jl b/test/materials/test_hyperelasticity_new_api.jl
new file mode 100644
index 0000000..46141f6
--- /dev/null
+++ b/test/materials/test_hyperelasticity_new_api.jl
@@ -0,0 +1,693 @@
+"""
+# Hyperelasticity - NEW API (Test-Driven Development)
+
+**What:** Shows how hyperelastic materials SHOULD work with the NEW API
+
+**Why:**
+- **Finite strain** - Large deformations (rubber, soft tissue, biomechanics)
+- **Energy-based** - Strain energy function Ψ(F)
+- **Frame-invariant** - Material objectivity (rotation independence)
+- **Multiple models** - Neo-Hookean, Mooney-Rivlin, Ogden, etc.
+- **Incompressibility** - Nearly incompressible (ν ≈ 0.5)
+
+**NEW API Concepts:**
+1. **Hyperelastic material types** - Neo-Hookean, Mooney-Rivlin, Ogden
+2. **Strain energy function** - Ψ(F) and derivatives
+3. **Push-forward stress** - σ = (1/J) P F^T (Cauchy from 1st Piola-Kirchhoff)
+4. **Tangent moduli** - C_ijkl = ∂²Ψ/∂F_ij∂F_kl
+5. **Incompressibility constraint** - det(F) ≈ 1
+
+**Test Problems:**
+
+## Test 1: Neo-Hookean (Simplest Hyperelastic)
+- Ψ = μ/2 (I₁ - 3) - μ ln(J) + λ/2 (ln J)²
+- Validates stress computation from energy
+- Tests incompressibility limit
+
+## Test 2: Mooney-Rivlin (Two-Parameter)
+- Ψ = C₁(I₁ - 3) + C₂(I₂ - 3) + κ/2 (J - 1)²
+- Better for rubber than Neo-Hookean
+- Validates second invariant I₂
+
+## Test 3: Ogden Model (Multi-Term)
+- Ψ = Σᵢ μᵢ/αᵢ (λ₁^αᵢ + λ₂^αᵢ + λ₃^αᵢ - 3)
+- Uses principal stretches λᵢ
+- Most accurate for rubber
+
+## Test 4: Uniaxial Tension Test
+- Compare to experimental data
+- Validates material parameters
+- Tests large strain (λ > 2)
+
+**Expected Behavior (when implemented):**
+✅ Stress computed from ∂Ψ/∂F correctly
+✅ Tangent moduli symmetric and positive-definite
+✅ Incompressibility enforced (J ≈ 1)
+✅ Frame-invariant (rotations don't change Ψ)
+✅ Matches experimental stress-strain curves
+✅ Works with Newton-Krylov solver
+
+**Status:** 🚧 VISIONARY TEST - Implementation in progress
+"""
+
+using Test
+using JuliaFEM
+using Tensors
+using LinearAlgebra
+using Statistics
+
+@testset "Hyperelasticity - NEW API (TDD)" begin
+
+    # =============================================================================
+    # NEO-HOOKEAN MODEL (SIMPLEST)
+    # =============================================================================
+
+    @testset "Neo-Hookean Material (Visionary)" begin
+        @test_skip begin  # Skip until implemented
+
+            # Material parameters
+            μ = 1000.0    # Shear modulus
+            λ = 2000.0    # Lame parameter (nearly incompressible)
+
+            # NEW: Hyperelastic material type
+            material = NeoHookean(
+                μ=μ,
+                λ=λ,
+                formulation=:compressible  # or :incompressible
+            )
+
+            # Test deformation gradient
+            F = Tensor{2,3}((
+                1.2, 0.1, 0.0,
+                0.0, 0.9, 0.0,
+                0.0, 0.0, 1.0
+            ))
+
+            # Compute strain energy
+            Ψ = strain_energy(material, F)
+
+            # Neo-Hookean energy:
+            # Ψ = μ/2 (I₁ - 3) - μ ln(J) + λ/2 (ln J)²
+            C = tdot(F)  # Right Cauchy-Green: C = F^T F
+            I₁ = tr(C)    # First invariant
+            J = det(F)    # Volume ratio
+
+            Ψ_analytical = μ / 2 * (I₁ - 3) - μ * log(J) + λ / 2 * log(J)^2
+
+            @test isapprox(Ψ, Ψ_analytical, rtol=1e-10)
+
+            # First Piola-Kirchhoff stress: P = ∂Ψ/∂F
+            P = first_piola_kirchhoff_stress(material, F)
+
+            # Analytical P for Neo-Hookean
+            F_inv = inv(F)
+            P_analytical = μ * (F - tdot(F_inv)) + λ * log(J) * tdot(F_inv)
+
+            @test P ≈ P_analytical rtol = 1e-10
+
+            # Cauchy stress: σ = (1/J) P F^T
+            σ = cauchy_stress(material, F)
+            σ_from_P = (1 / J) * P ⊡ transpose(F)
+
+            @test σ ≈ σ_from_P rtol = 1e-10
+
+        end
+    end
+
+    # =============================================================================
+    # INCOMPRESSIBILITY CONSTRAINT
+    # =============================================================================
+
+    @testset "Incompressibility (Nearly) (Visionary)" begin
+        @test_skip begin
+
+            # Nearly incompressible (ν → 0.5)
+            E = 1000.0
+            ν = 0.499  # Nearly incompressible
+
+            μ = E / (2 * (1 + ν))
+            λ = E * ν / ((1 + ν) * (1 - 2ν))  # Very large!
+
+            material = NeoHookean(
+                μ=μ,
+                λ=λ,
+                formulation=:compressible
+            )
+
+            # Deformation (should preserve volume)
+            F = Tensor{2,3}((
+                1.5, 0.0, 0.0,
+                0.0, 1 / sqrt(1.5), 0.0,
+                0.0, 0.0, 1 / sqrt(1.5)
+            ))
+
+            J = det(F)
+
+            # For incompressible: J = 1 (volume preserving)
+            @test isapprox(J, 1.0, atol=1e-3)
+
+            # Hydrostatic pressure enforces incompressibility
+            σ = cauchy_stress(material, F)
+            p = -tr(σ) / 3  # Hydrostatic pressure
+
+            # For nearly incompressible, pressure should be large
+            @test abs(p) > 1000.0  # Significant pressure
+
+        end
+    end
+
+    # =============================================================================
+    # MOONEY-RIVLIN MODEL
+    # =============================================================================
+
+    @testset "Mooney-Rivlin Material (Visionary)" begin
+        @test_skip begin
+
+            # Material parameters (typical for rubber)
+            C₁ = 0.5  # MPa
+            C₂ = 0.1  # MPa
+            κ = 100.0  # Bulk modulus (incompressibility)
+
+            # NEW: Mooney-Rivlin material
+            material = MooneyRivlin(
+                C1=C₁,
+                C2=C₂,
+                bulk_modulus=κ
+            )
+
+            # Test deformation
+            F = Tensor{2,3}((
+                1.5, 0.2, 0.0,
+                0.1, 0.8, 0.0,
+                0.0, 0.0, 1.0
+            ))
+
+            # Strain energy: Ψ = C₁(I₁ - 3) + C₂(I₂ - 3) + κ/2 (J - 1)²
+            C = tdot(F)
+            I₁ = tr(C)
+            I₂ = 0.5 * (tr(C)^2 - tr(C ⊡ C))  # Second invariant
+            J = det(F)
+
+            Ψ = strain_energy(material, F)
+            Ψ_analytical = C₁ * (I₁ - 3) + C₂ * (I₂ - 3) + κ / 2 * (J - 1)^2
+
+            @test isapprox(Ψ, Ψ_analytical, rtol=1e-10)
+
+            # Stress computation
+            P = first_piola_kirchhoff_stress(material, F)
+            σ = cauchy_stress(material, F)
+
+            # Validate symmetry of Cauchy stress
+            @test isapprox(σ, transpose(σ), atol=1e-10)
+
+        end
+    end
+
+    # =============================================================================
+    # OGDEN MODEL (PRINCIPAL STRETCHES)
+    # =============================================================================
+
+    @testset "Ogden Material (Visionary)" begin
+        @test_skip begin
+
+            # Ogden parameters (multi-term)
+            μ_terms = [1.0, 0.5, 0.2]    # Shear moduli
+            α_terms = [2.0, 3.0, -2.0]   # Exponents
+            κ = 100.0
+
+            # NEW: Ogden material
+            material = Ogden(
+                mu=μ_terms,
+                alpha=α_terms,
+                bulk_modulus=κ
+            )
+
+            # Deformation
+            F = Tensor{2,3}((
+                2.0, 0.0, 0.0,
+                0.0, 0.6, 0.0,
+                0.0, 0.0, 0.8
+            ))
+
+            # Compute principal stretches
+            C = tdot(F)
+            eigenvalues_C = eigvals(C)
+            λ = sqrt.(eigenvalues_C)  # Principal stretches
+
+            # Strain energy: Ψ = Σᵢ μᵢ/αᵢ (λ₁^αᵢ + λ₂^αᵢ + λ₃^αᵢ - 3)
+            Ψ = strain_energy(material, F)
+
+            J = det(F)
+            Ψ_analytical = sum(
+                μ_terms[i] / α_terms[i] * (sum(λ .^ α_terms[i]) - 3)
+                for i in 1:3
+            ) + κ / 2 * (J - 1)^2
+
+            @test isapprox(Ψ, Ψ_analytical, rtol=1e-10)
+
+        end
+    end
+
+    # =============================================================================
+    # UNIAXIAL TENSION TEST
+    # =============================================================================
+
+    @testset "Uniaxial Tension (Large Strain) (Visionary)" begin
+        @test_skip begin
+
+            # Geometry: Unit cube under tension
+            mesh = generate_mesh(
+                geometry=UnitCube(),
+                element_type=Hex8,
+                n_elements=(4, 4, 4)
+            )
+
+            # Material (Neo-Hookean rubber)
+            material = NeoHookean(
+                μ=1.0,    # MPa
+                λ=10.0,   # Nearly incompressible
+                density=1000.0
+            )
+
+            # Physics
+            elastic_physics = ContinuumPhysics{Displacement}(
+                material=material,
+                formulation=FullThreeD(),
+                finite_strain=true  # CRITICAL!
+            )
+
+            domain = Domain(
+                name="RUBBER",
+                elements=mesh.elements,
+                physics=elastic_physics
+            )
+
+            # Boundary conditions
+            bc_fixed = DirichletBC(
+                nodes=mesh.node_sets["LEFT"],
+                dof=:displacement,
+                values=[0.0, 0.0, 0.0]
+            )
+
+            # Applied stretch (λ = 2.0, 100% strain!)
+            u_applied = 1.0  # Stretch from 1.0 to 2.0
+
+            bc_stretch = DirichletBC(
+                nodes=mesh.node_sets["RIGHT"],
+                dof=:displacement,
+                component=:x,
+                value=u_applied
+            )
+
+            # Nonlinear problem (finite strain)
+            problem = NonlinearProblem(
+                domains=[domain],
+                boundary_conditions=[bc_fixed, bc_stretch]
+            )
+
+            # Newton-Krylov solver
+            solver = NewtonKrylovSolver(
+                max_iterations=20,
+                convergence_tol=1e-6,
+                krylov_solver=GMRES(restart=30),
+                line_search=BacktrackingLineSearch()
+            )
+
+            solution = solve!(problem, solver)
+
+            # Extract stress-stretch curve
+            λ = 1.0 + u_applied  # Stretch ratio
+
+            # Compute engineering stress
+            F_avg = compute_average_deformation_gradient(solution)
+            P_avg = first_piola_kirchhoff_stress(material, F_avg)
+
+            # Engineering stress: σ_eng = P_11 (1st Piola-Kirchhoff in x)
+            σ_eng = P_avg[1, 1]
+
+            # For Neo-Hookean uniaxial:
+            # σ_eng = μ(λ - 1/λ²)
+            σ_analytical = material.μ * (λ - 1 / λ^2)
+
+            @test isapprox(σ_eng, σ_analytical, rtol=0.05)  # 5% tolerance
+
+        end
+    end
+
+    # =============================================================================
+    # FRAME INVARIANCE (OBJECTIVITY)
+    # =============================================================================
+
+    @testset "Frame Invariance (Visionary)" begin
+        @test_skip begin
+
+            material = NeoHookean(μ=1000.0, λ=2000.0)
+
+            # Deformation gradient
+            F = Tensor{2,3}((
+                1.2, 0.1, 0.0,
+                0.0, 0.9, 0.0,
+                0.0, 0.0, 1.0
+            ))
+
+            # Rotation tensor (90° about z-axis)
+            θ = π / 2
+            Q = Tensor{2,3}((
+                cos(θ), -sin(θ), 0.0,
+                sin(θ), cos(θ), 0.0,
+                0.0, 0.0, 1.0
+            ))
+
+            # Rotated deformation: F' = Q F
+            F_rotated = Q ⊡ F
+
+            # Strain energy should be INVARIANT
+            Ψ = strain_energy(material, F)
+            Ψ_rotated = strain_energy(material, F_rotated)
+
+            @test isapprox(Ψ, Ψ_rotated, rtol=1e-10)
+
+            # Cauchy stress should transform: σ' = Q σ Q^T
+            σ = cauchy_stress(material, F)
+            σ_rotated = cauchy_stress(material, F_rotated)
+
+            σ_transformed = Q ⊡ σ ⊡ transpose(Q)
+
+            @test σ_rotated ≈ σ_transformed rtol = 1e-10
+
+        end
+    end
+
+    # =============================================================================
+    # TANGENT MODULI (FOR NEWTON)
+    # =============================================================================
+
+    @testset "Tangent Moduli (Visionary)" begin
+        @test_skip begin
+
+            material = NeoHookean(μ=1000.0, λ=2000.0)
+
+            F = Tensor{2,3}((
+                1.2, 0.1, 0.0,
+                0.0, 0.9, 0.0,
+                0.0, 0.0, 1.0
+            ))
+
+            # Material tangent: C_ijkl = ∂²Ψ/∂F_ij∂F_kl
+            C = material_tangent(material, F)
+
+            # Validate major symmetry: C_ijkl = C_klij
+            # (Minor symmetries don't hold for finite strain)
+            for i in 1:3, j in 1:3, k in 1:3, l in 1:3
+                @test isapprox(C[i, j, k, l], C[k, l, i, j], atol=1e-10)
+            end
+
+            # Validate positive-definiteness
+            # For small perturbation δF, δ²Ψ = C_ijkl δF_ij δF_kl > 0
+            δF = 0.01 * rand(Tensor{2,3})
+            δ²Ψ = dcontract(dcontract(C, δF), δF)
+
+            @test δ²Ψ > 0  # Positive-definite
+
+        end
+    end
+
+    # =============================================================================
+    # PSEUDO-CODE: HYPERELASTIC ASSEMBLY
+    # =============================================================================
+
+    @testset "Hyperelastic Assembly Pattern (Visionary)" begin
+        # Pseudo-code showing finite strain assembly
+
+        println("\n" * "="^70)
+        println("HYPERELASTIC ASSEMBLY (FINITE STRAIN)")
+        println("="^70)
+
+        assembly_pseudo = """
+        # For hyperelastic materials, assembly uses current configuration
+
+        function compute_residual_hyperelastic(u, material, elements)
+            r = zeros(length(u))
+            
+            for elem in elements
+                # Current deformation gradient: F = I + ∇u
+                X = reference_coordinates(elem)  # Undeformed
+                x = X + u[elem.nodes]            # Deformed
+                
+                for ip in integration_points(elem)
+                    # Jacobian in reference config
+                    J₀, dN_dX = jacobian(elem, ip, X)
+                    
+                    # Deformation gradient: F = ∂x/∂X
+                    F = compute_deformation_gradient(x, dN_dX)
+                    
+                    # 1st Piola-Kirchhoff stress: P = ∂Ψ/∂F
+                    P = first_piola_kirchhoff_stress(material, F)
+                    
+                    # Residual: r = ∫_Ω₀ P : ∇_X(δu) dΩ₀
+                    # = ∫_Ω₀ P_iJ (dN_I/dX_J) dΩ₀
+                    for I in 1:n_nodes
+                        for i in 1:3
+                            for J in 1:3
+                                r[3*(I-1)+i] += P[i,J] * dN_dX[I,J] * J₀ * ip.weight
+                            end
+                        end
+                    end
+                end
+            end
+            
+            return r
+        end
+
+        # Tangent stiffness (for Newton):
+        # K_t = ∫_Ω₀ (dN_I/dX_K) C_iJkL (dN_J/dX_L) dΩ₀
+        # where C_iJkL = ∂²Ψ/∂F_iJ∂F_kL
+        """
+
+        println(assembly_pseudo)
+        println("="^70)
+        println("✓ Uses reference configuration Ω₀ (not current!)")
+        println("✓ Deformation gradient F = ∂x/∂X")
+        println("✓ 1st Piola-Kirchhoff stress P = ∂Ψ/∂F")
+        println("✓ Tangent moduli C = ∂²Ψ/∂F²")
+        println("✓ Works with nodal assembly (same pattern!)")
+        println("="^70)
+    end
+
+    # =============================================================================
+    # KEY ARCHITECTURAL INSIGHTS
+    # =============================================================================
+
+    println("\n" * "="^70)
+    println("HYPERELASTICITY ARCHITECTURE INSIGHTS (NEW API)")
+    println("="^70)
+    println("✓ Hyperelastic materials: NeoHookean, MooneyRivlin, Ogden")
+    println("✓ Strain energy function Ψ(F) is fundamental")
+    println("✓ Stress from energy: P = ∂Ψ/∂F, σ = (1/J) P F^T")
+    println("✓ Tangent from energy: C = ∂²Ψ/∂F²")
+    println("✓ Frame-invariant: Rotations don't change Ψ")
+    println("✓ Incompressibility: det(F) ≈ 1 for rubber")
+    println("✓ Works with Newton-Krylov (unsymmetric OK)")
+    println("✓ Assembly in reference config (not current!)")
+    println("✓ Large strains: λ > 2 (100%+ strain)")
+    println("="^70)
+
+end
+
+"""
+# IMPLEMENTATION NOTES
+
+## Hyperelastic Material Models
+
+### Neo-Hookean (Simplest)
+
+**Strain energy:**
+Ψ = μ/2 (I₁ - 3) - μ ln(J) + λ/2 (ln J)²
+
+where:
+- I₁ = tr(C) = tr(F^T F) (first invariant)
+- J = det(F) (volume ratio)
+- μ, λ = Lame parameters
+
+**1st Piola-Kirchhoff stress:**
+P = ∂Ψ/∂F = μ(F - F^{-T}) + λ ln(J) F^{-T}
+
+where F^{-T} = (F^{-1})^T.
+
+**Cauchy stress:**
+σ = (1/J) P F^T = μ/J (B - I) + λ/J ln(J) I
+
+where B = F F^T (left Cauchy-Green).
+
+**Use case:** Rubber at moderate strains (< 100%).
+
+### Mooney-Rivlin (Two-Parameter)
+
+**Strain energy:**
+Ψ = C₁(I₁ - 3) + C₂(I₂ - 3) + κ/2 (J - 1)²
+
+where:
+- I₁ = tr(C)
+- I₂ = 1/2 [(tr C)² - tr(C²)]
+- J = det(F)
+- C₁, C₂ = material parameters
+- κ = bulk modulus
+
+**Better fit for rubber** than Neo-Hookean.
+
+**Relation to Neo-Hookean:** C₂ = 0 → Neo-Hookean.
+
+### Ogden Model (Multi-Term)
+
+**Strain energy:**
+Ψ = Σᵢ μᵢ/αᵢ (λ₁^{αᵢ} + λ₂^{αᵢ} + λ₃^{αᵢ} - 3) + κ/2 (J - 1)²
+
+where:
+- λ₁, λ₂, λ₃ = principal stretches (eigenvalues of F)
+- μᵢ, αᵢ = material parameters (typically 3-6 terms)
+
+**Most accurate for rubber** (fits experimental data well).
+
+**Implementation:** Requires spectral decomposition of C.
+
+## Stress Measures
+
+### First Piola-Kirchhoff (P)
+
+**Definition:** P = ∂Ψ/∂F
+
+**Properties:**
+- Non-symmetric
+- Force per reference area
+- Work-conjugate to F
+
+**Use:** Weak form in reference config.
+
+### Cauchy Stress (σ)
+
+**Definition:** σ = (1/J) P F^T
+
+**Properties:**
+- Symmetric
+- True stress (force per current area)
+- What we measure
+
+**Use:** Post-processing, failure criteria.
+
+### Second Piola-Kirchhoff (S)
+
+**Definition:** S = F^{-1} P = J F^{-1} σ F^{-T}
+
+**Properties:**
+- Symmetric
+- Work-conjugate to E (Green-Lagrange strain)
+- Energy-conjugate
+
+**Use:** Theoretical derivations.
+
+## Incompressibility
+
+**Constraint:** det(F) = J = 1 (volume preserving)
+
+**Nearly incompressible:** ν → 0.5, λ → ∞
+
+**Enforcement:**
+
+1. **Penalty:** Add κ/2 (J - 1)² to Ψ (large κ)
+2. **Lagrange multiplier:** Introduce pressure p
+3. **Mixed formulation:** (u, p) unknowns
+
+**JuliaFEM approach:** Penalty for compressible materials, mixed for truly incompressible.
+
+## Frame Invariance (Objectivity)
+
+**Definition:** Strain energy invariant under rigid rotations.
+
+**Mathematical:** Ψ(Q F) = Ψ(F) for all rotations Q.
+
+**Why:** Material doesn't "know" about global rotations.
+
+**Implementation:** Use invariants (I₁, I₂, I₃) or principal stretches (λᵢ).
+
+**Validation:**
+```julia
+Q = rotation_matrix(θ)
+F_rotated = Q * F
+@test strain_energy(F_rotated) ≈ strain_energy(F)
+```
+
+## Tangent Moduli
+
+**Material tangent:**
+C_{iJkL} = ∂²Ψ/∂F_{iJ}∂F_{kL}
+
+**Spatial tangent:**
+c_{ijkl} = (1/J) F_{iI} F_{jJ} F_{kK} F_{lL} C_{IJKL}
+
+**Symmetries:**
+- Major: C_{iJkL} = C_{kLiJ} (from Ψ)
+- Minor: Generally NOT symmetric in hyperelasticity
+
+**Use in Newton:**
+δP = C : δF
+
+## Assembly (Finite Strain)
+
+**Residual (weak form):**
+r = ∫_{Ω₀} P : ∇_X(δu) dΩ₀ - f_ext
+
+**In components:**
+r_I^i = ∫_{Ω₀} P_{iJ} ∂N_I/∂X_J dΩ₀ - f_I^i
+
+**Tangent stiffness:**
+K_{IJ}^{ik} = ∫_{Ω₀} ∂N_I/∂X_K C_{iJkL} ∂N_J/∂X_L dΩ₀
+
+**Key differences from small strain:**
+1. Integrate over Ω₀ (reference), not Ω (current)
+2. Use ∂N/∂X (reference gradients), not ∂N/∂x
+3. P (1st Piola-Kirchhoff), not σ (Cauchy)
+
+## Nodal Assembly (Hyperelastic)
+
+```julia
+function tangent_matvec_hyperelastic!(w, v, u_current, material, elements, node_to_elements)
+    Threads.@threads for node_i in 1:n_nodes
+        w_local = zero(Vec{3})
+        
+        for elem in node_to_elements[node_i]
+            # Current deformation
+            F = deformation_gradient(elem, u_current)
+            
+            # Material tangent
+            C = material_tangent(material, F)
+            
+            for node_j in elem.nodes
+                # Tangent block: K_t_ij
+                K_t_ij = compute_hyperelastic_tangent_block(elem, node_i, node_j, C, F)
+                
+                v_j = Vec{3}(v[3*(node_j-1)+1:3*node_j])
+                w_local += K_t_ij ⊡ v_j
+            end
+        end
+        
+        w[3*(node_i-1)+1:3*node_i] = w_local
+    end
+end
+```
+
+**Same pattern as linear elasticity!** Only the tangent computation changes.
+
+## Next Steps
+
+1. Implement `NeoHookean` material type
+2. Implement `strain_energy` function
+3. Implement `first_piola_kirchhoff_stress`
+4. Implement `cauchy_stress` (push-forward)
+5. Implement `material_tangent` (C_ijkl)
+6. Implement `MooneyRivlin` material
+7. Implement `Ogden` material (spectral decomposition)
+8. Validate against analytical solutions
+9. Validate against experimental data
+10. Performance benchmarks
+
+"""