diff --git a/test/elements/test_celement_real_basis.jl b/test/elements/test_celement_real_basis.jl new file mode 100644 index 0000000..664b833 --- /dev/null +++ b/test/elements/test_celement_real_basis.jl @@ -0,0 +1,285 @@ +using Test +using JuliaFEM +using JuliaFEM: CElement, ScalarDOF, VectorDOF, interpolate, gradient +using Tensors + +""" +Test suite for CElement with REAL basis functions (not stubs). + +This test file verifies that CElement correctly uses the actual Lagrange +basis functions from basis_generated.jl, not the stub implementations. +""" + +@testset "CElement Real Basis Functions" begin + + @testset "Triangle Tri3 Linear Interpolation" begin + # Create a linear triangle element + # Nodes at vertices: (0,0), (1,0), (0,1) + nodes = [ + Vec{2}((0.0, 0.0)), # Node 1 + Vec{2}((1.0, 0.0)), # Node 2 + Vec{2}((0.0, 1.0)) # Node 3 + ] + + # Simple mesh structure + mesh = ( + connectivity = [[1, 2, 3]], + nodes = nodes + ) + + # Create element with scalar DOF + elem = CElement{Tri3, Lagrange{1}, ScalarDOF}(1, (1, 2, 3)) + + # Test nodal interpolation: At nodes, basis functions should be 1 at that node, 0 elsewhere + # At node 1: ξ = (0, 0) → should get u[1] exactly + u_global = [10.0, 20.0, 30.0] # Values at nodes 1, 2, 3 + + # At parametric coord (0, 0) = node 1 + val1 = interpolate(elem, mesh, u_global, Vec{2}((0.0, 0.0))) + @test val1 ≈ 10.0 atol=1e-12 + + # At parametric coord (1, 0) = node 2 + val2 = interpolate(elem, mesh, u_global, Vec{2}((1.0, 0.0))) + @test val2 ≈ 20.0 atol=1e-12 + + # At parametric coord (0, 1) = node 3 + val3 = interpolate(elem, mesh, u_global, Vec{2}((0.0, 1.0))) + @test val3 ≈ 30.0 atol=1e-12 + + # Test centroid: ξ = (1/3, 1/3) should give average + centroid = interpolate(elem, mesh, u_global, Vec{2}(1.0/3.0, 1.0/3.0)) + expected_centroid = (10.0 + 20.0 + 30.0) / 3.0 + @test centroid ≈ expected_centroid atol=1e-12 + + println("✓ Tri3 linear interpolation: Exact at nodes, correct at centroid") + end + + @testset "Triangle Tri3 Linear Gradient" begin + # Same triangle as above + nodes = [ + Vec{2}((0.0, 0.0)), + Vec{2}((1.0, 0.0)), + Vec{2}((0.0, 1.0)) + ] + + mesh = ( + connectivity = [[1, 2, 3]], + nodes = nodes + ) + + elem = CElement{Tri3, Lagrange{1}, ScalarDOF}(1, (1, 2, 3)) + + # Linear field: u = x + 2y (gradient should be [1, 2]) + # At node 1 (0,0): u = 0 + # At node 2 (1,0): u = 1 + # At node 3 (0,1): u = 2 + u_global = [0.0, 1.0, 2.0] + + # Gradient should be constant [1, 2] everywhere for linear element + grad_center = gradient(elem, mesh, u_global, Vec{2}(1.0/3.0, 1.0/3.0)) + @test grad_center[1] ≈ 1.0 atol=1e-10 + @test grad_center[2] ≈ 2.0 atol=1e-10 + + # Gradient should be same at different points + grad_node1 = gradient(elem, mesh, u_global, Vec{2}((0.0, 0.0))) + @test grad_node1[1] ≈ 1.0 atol=1e-10 + @test grad_node1[2] ≈ 2.0 atol=1e-10 + + println("✓ Tri3 gradient: Constant for linear field") + end + + @testset "Tetrahedron Tet4 Linear Interpolation" begin + # Create a linear tetrahedron + # Standard reference tet: (0,0,0), (1,0,0), (0,1,0), (0,0,1) + nodes = [ + Vec{3}((0.0, 0.0, 0.0)), # Node 1 + Vec{3}((1.0, 0.0, 0.0)), # Node 2 + Vec{3}((0.0, 1.0, 0.0)), # Node 3 + Vec{3}((0.0, 0.0, 1.0)) # Node 4 + ] + + mesh = ( + connectivity = [[1, 2, 3, 4]], + nodes = nodes + ) + + elem = CElement{Tet4, Lagrange{1}, ScalarDOF}(1, (1, 2, 3, 4)) + + # Test nodal interpolation + u_global = [5.0, 15.0, 25.0, 35.0] + + # At node 1: ξ = (0, 0, 0) + val1 = interpolate(elem, mesh, u_global, Vec{3}((0.0, 0.0, 0.0))) + @test val1 ≈ 5.0 atol=1e-12 + + # At node 2: ξ = (1, 0, 0) + val2 = interpolate(elem, mesh, u_global, Vec{3}((1.0, 0.0, 0.0))) + @test val2 ≈ 15.0 atol=1e-12 + + # At node 3: ξ = (0, 1, 0) + val3 = interpolate(elem, mesh, u_global, Vec{3}((0.0, 1.0, 0.0))) + @test val3 ≈ 25.0 atol=1e-12 + + # At node 4: ξ = (0, 0, 1) + val4 = interpolate(elem, mesh, u_global, Vec{3}((0.0, 0.0, 1.0))) + @test val4 ≈ 35.0 atol=1e-12 + + # Test centroid: ξ = (1/4, 1/4, 1/4) + centroid = interpolate(elem, mesh, u_global, Vec{3}((0.25, 0.25, 0.25))) + expected = (5.0 + 15.0 + 25.0 + 35.0) / 4.0 + @test centroid ≈ expected atol=1e-12 + + println("✓ Tet4 linear interpolation: Exact at nodes, correct at centroid") + end + + @testset "Tetrahedron Tet4 Linear Gradient" begin + # Same tet as above + nodes = [ + 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)) + ] + + mesh = ( + connectivity = [[1, 2, 3, 4]], + nodes = nodes + ) + + elem = CElement{Tet4, Lagrange{1}, ScalarDOF}(1, (1, 2, 3, 4)) + + # Linear field: u = 2x + 3y + 4z → gradient = [2, 3, 4] + # At nodes: u = [0, 2, 3, 4] + u_global = [0.0, 2.0, 3.0, 4.0] + + # Gradient should be constant [2, 3, 4] everywhere + grad_center = gradient(elem, mesh, u_global, Vec{3}((0.25, 0.25, 0.25))) + @test grad_center[1] ≈ 2.0 atol=1e-10 + @test grad_center[2] ≈ 3.0 atol=1e-10 + @test grad_center[3] ≈ 4.0 atol=1e-10 + + println("✓ Tet4 gradient: Constant for linear field") + end + + @testset "Vector DOF Deformation Gradient" begin + # Test 2D vector DOF (displacement field) + nodes = [ + Vec{2}((0.0, 0.0)), + Vec{2}((1.0, 0.0)), + Vec{2}((0.0, 1.0)) + ] + + mesh = ( + connectivity = [[1, 2, 3]], + nodes = nodes + ) + + # Element with 2D vector DOF + elem = CElement{Tri3, Lagrange{1}, VectorDOF{2}}(1, (1, 2, 3, 4, 5, 6)) + + # Displacement field: u = [x, 2y] → F = [∂u₁/∂x ∂u₁/∂y; ∂u₂/∂x ∂u₂/∂y] = [1 0; 0 2] + # At nodes: u = [[0,0], [1,0], [0,2]] + u_global = [0.0, 0.0, # Node 1: (ux, uy) + 1.0, 0.0, # Node 2 + 0.0, 2.0] # Node 3 + + # Deformation gradient + F = gradient(elem, mesh, u_global, Vec{2}(1.0/3.0, 1.0/3.0)) + + # Should be [1 0; 0 2] for linear displacement + @test F[1,1] ≈ 1.0 atol=1e-10 # ∂u₁/∂x + @test F[1,2] ≈ 0.0 atol=1e-10 # ∂u₁/∂y + @test F[2,1] ≈ 0.0 atol=1e-10 # ∂u₂/∂x + @test F[2,2] ≈ 2.0 atol=1e-10 # ∂u₂/∂y + + println("✓ VectorDOF deformation gradient: Correct for linear displacement") + end + + @testset "Partition of Unity (Completeness)" begin + # Basis functions should sum to 1 at any point + nodes = [ + Vec{2}((0.0, 0.0)), + Vec{2}((1.0, 0.0)), + Vec{2}((0.0, 1.0)) + ] + + mesh = ( + connectivity = [[1, 2, 3]], + nodes = nodes + ) + + elem = CElement{Tri3, Lagrange{1}, ScalarDOF}(1, (1, 2, 3)) + + # Test at several points + test_points = [ + Vec{2}((0.0, 0.0)), + Vec{2}((0.5, 0.0)), + Vec{2}((0.0, 0.5)), + Vec{2}((0.5, 0.5)), + Vec{2}((0.33, 0.33)), + Vec{2}((0.1, 0.7)) + ] + + u_ones = [1.0, 1.0, 1.0] # If all nodal values = 1, result should be 1 + + for ξ in test_points + # Skip if outside element (u + v > 1) + if ξ[1] + ξ[2] > 1.0 + continue + end + + val = interpolate(elem, mesh, u_ones, ξ) + @test val ≈ 1.0 atol=1e-12 + end + + println("✓ Partition of unity: Sum of basis = 1 at all points") + end + + @testset "Linear Reproduction (Consistency)" begin + # Linear functions should be reproduced exactly + nodes = [ + Vec{2}((0.0, 0.0)), + Vec{2}((2.0, 0.0)), + Vec{2}((0.0, 3.0)) + ] + + mesh = ( + connectivity = [[1, 2, 3]], + nodes = nodes + ) + + elem = CElement{Tri3, Lagrange{1}, ScalarDOF}(1, (1, 2, 3)) + + # Linear function: u(x,y) = 5 + 2x + 3y + u_nodal = [ + 5.0, # At (0,0): 5 + 5.0 + 2.0*2.0, # At (2,0): 9 + 5.0 + 3.0*3.0 # At (0,3): 14 + ] + + # Test at arbitrary physical points + test_points = [ + (Vec{2}((0.5, 0.0)), Vec{2}((1.0, 0.0))), # (ξ, physical) + (Vec{2}((0.0, 0.5)), Vec{2}((0.0, 1.5))), + (Vec{2}((0.25, 0.25)), Vec{2}((0.5, 0.75))) + ] + + for (ξ, x_phys) in test_points + val = interpolate(elem, mesh, u_nodal, ξ) + expected = 5.0 + 2.0*x_phys[1] + 3.0*x_phys[2] + @test val ≈ expected atol=1e-10 + end + + println("✓ Linear reproduction: Linear functions reproduced exactly") + end + +end # @testset "CElement Real Basis Functions" + +println("\n" * "="^70) +println("CElement Real Basis Test Summary") +println("="^70) +println("✅ All tests verify that CElement uses REAL Lagrange basis functions") +println("✅ Interpolation: Exact at nodes, correct partition of unity") +println("✅ Gradient: Constant for linear elements, correct deformation gradient") +println("✅ Math properties: Completeness and consistency verified") +println("="^70)