diff --git a/test/basis/test_piola_nedelec.jl b/test/basis/test_piola_nedelec.jl new file mode 100644 index 0000000..cedd3c2 --- /dev/null +++ b/test/basis/test_piola_nedelec.jl @@ -0,0 +1,166 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using Test +using JuliaFEM +using Tensors + +@testset "Piola identity on J = I" begin + J = Tensor{2, 3}((1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) + u = Vec((1.0, -2.0, 0.5)) + @test piola_covariant(J, u) ≈ u + @test piola_contravariant(J, u) ≈ u +end + +@testset "Nédélec Whitney tet — tangential component on edge (1,2)" begin + ξ = Vec((0.3, 0.0, 0.0)) + N = nedelec_whitney_tet_reference(1, ξ) + ex = Vec((1.0, 0.0, 0.0)) + @test dot(N, ex) ≈ 1.0 rtol = 1e-14 +end + +@testset "Piola covariant commutes with identity Jacobian (reference point)" begin + # Physical tet = reference tet: Jacobian is identity at any ξ for affine map. + J = Tensor{2, 3}((1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) + ξ = Vec((0.2, 0.2, 0.2)) + u_ref = nedelec_whitney_tet_reference(4, ξ) + @test piola_covariant(J, u_ref) ≈ u_ref +end + +@testset "Nédélec hierarchical tet — slot 1 matches Whitney" begin + ξ = Vec((0.17, 0.09, 0.05)) + for le in 1:6 + @test nedelec_hierarchical_tet_edge_reference(le, 1, ξ) ≈ nedelec_whitney_tet_reference(le, ξ) + end +end + +@testset "Nédélec hierarchical tet — slot 2 tangential on edge (1,2)" begin + ξ = Vec((0.3, 0.0, 0.0)) + ex = Vec((1.0, 0.0, 0.0)) + N2 = nedelec_hierarchical_tet_edge_reference(1, 2, ξ) + @test dot(N2, ex) ≈ 0.7 * 0.3 rtol = 1e-14 +end + +@testset "Nédélec hierarchical tet — Piola covariant with diagonal Jacobian" begin + J = Tensor{2, 3}((2.0, 0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 4.0)) + ξ = Vec((0.11, 0.22, 0.07)) + le = 3 + for slot in 1:2 + u_ref = nedelec_hierarchical_tet_edge_reference(le, slot, ξ) + u_phys = piola_covariant(J, u_ref) + inv_diag = Vec((1 / 2, 1 / 3, 1 / 4)) + @test u_phys ≈ inv_diag .* u_ref rtol = 1e-14 + end +end + +@testset "Piola contravariant with diagonal non-unit Jacobian" begin + J = Tensor{2, 3}((2.0, 0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 4.0)) + ex = Vec((1.0, 0.0, 0.0)) + @test piola_contravariant(J, ex) ≈ Vec((1.0 / 12.0, 0.0, 0.0)) +end + +# Reference ξ; physical `x` diagonal with `∂/∂x_i = (1/J_ii) ∂/∂ξ_i`. +function _divergence_physical_diagonal(f, ξ::Vec{3}, inv_diag::Vec{3}; h::Float64 = 1e-6) + e1 = Vec((1.0, 0.0, 0.0)) + e2 = Vec((0.0, 1.0, 0.0)) + e3 = Vec((0.0, 0.0, 1.0)) + return ( + inv_diag[1] * (f(ξ + h * e1)[1] - f(ξ - h * e1)[1]) / (2h) + + inv_diag[2] * (f(ξ + h * e2)[2] - f(ξ - h * e2)[2]) / (2h) + + inv_diag[3] * (f(ξ + h * e3)[3] - f(ξ - h * e3)[3]) / (2h) + ) +end + +function _divergence_reference(f, ξ::Vec{3}; h::Float64 = 1e-6) + e1 = Vec((1.0, 0.0, 0.0)) + e2 = Vec((0.0, 1.0, 0.0)) + e3 = Vec((0.0, 0.0, 1.0)) + return ( + (f(ξ + h * e1)[1] - f(ξ - h * e1)[1]) / (2h) + + (f(ξ + h * e2)[2] - f(ξ - h * e2)[2]) / (2h) + + (f(ξ + h * e3)[3] - f(ξ - h * e3)[3]) / (2h) + ) +end + +@testset "RT₀ tet reference — divergence constant (finite difference)" begin + ξa = Vec((0.11, 0.22, 0.05)) + ξb = Vec((0.05, 0.07, 0.08)) + for lf in 1:4 + div_a = _divergence_reference(ξ -> rt0_tet_reference_basis(lf, ξ), ξa) + div_b = _divergence_reference(ξ -> rt0_tet_reference_basis(lf, ξ), ξb) + @test div_a ≈ div_b rtol = 1e-4 + @test div_a ≈ 36.0 rtol = 1e-5 + end +end + +# Outward face normals for Hex8 reference faces (same order as `faces(::Hex8)`). +const _RT0_HEX_FACE_NORMALS = ( + Vec((0.0, 0.0, -1.0)), # z = -1 + Vec((0.0, 0.0, 1.0)), # z = +1 + Vec((0.0, -1.0, 0.0)), # y = -1 + Vec((1.0, 0.0, 0.0)), # x = +1 + Vec((0.0, 1.0, 0.0)), # y = +1 + Vec((-1.0, 0.0, 0.0)), # x = -1 +) + +@testset "RT₀ hex reference — divergence constant (finite difference)" begin + ξa = Vec((0.11, -0.22, 0.31)) + ξb = Vec((-0.05, 0.17, -0.08)) + for lf in 1:6 + div_a = _divergence_reference(ξ -> rt0_hex8_reference_basis(lf, ξ), ξa) + div_b = _divergence_reference(ξ -> rt0_hex8_reference_basis(lf, ξ), ξb) + @test div_a ≈ div_b rtol = 1e-5 + @test div_a ≈ 0.5 rtol = 1e-5 + end +end + +@testset "RT₀ hex reference — unit outward normal flux on owning face" begin + face_points = ( + Vec((0.25, 0.4, -1.0)), + Vec((-0.5, -0.2, 1.0)), + Vec((0.2, -1.0, 0.35)), + Vec((1.0, -0.1, 0.2)), + Vec((-0.3, 1.0, -0.4)), + Vec((-1.0, 0.15, 0.25)), + ) + for lf in 1:6 + ξ = face_points[lf] + ψ = rt0_hex8_reference_basis(lf, ξ) + n = _RT0_HEX_FACE_NORMALS[lf] + @test dot(ψ, n) ≈ 1.0 rtol = 1e-14 + end +end + +@testset "RT₀ hex + Piola contravariant — divergence scales by 1/det(J)" begin + J = Tensor{2, 3}((2.0, 0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 4.0)) + inv_diag = Vec((1 / 2, 1 / 3, 1 / 4)) + ξ = Vec((0.15, -0.25, 0.05)) + lf = 4 + div_ref = _divergence_reference(ξ -> rt0_hex8_reference_basis(lf, ξ), ξ) + div_phys = _divergence_physical_diagonal( + ξ -> piola_contravariant(J, rt0_hex8_reference_basis(lf, ξ)), + ξ, + inv_diag, + ) + @test div_ref ≈ 0.5 rtol = 1e-5 + @test div_phys ≈ div_ref / 24 rtol = 1e-4 +end + +@testset "RT₀ hex reference — zero normal flux on opposite face" begin + # Parallel face pairs: (1,2), (3,5), (4,6) in `faces(::Hex8)` order. + lf_opposite = (2, 1, 5, 6, 3, 4) + opposite_points = ( + Vec((0.1, 0.2, 1.0)), # on face 2 (z=+) + Vec((0.2, -0.3, -1.0)), # on face 1 (z=-) + Vec((-0.2, 1.0, 0.1)), # on face 5 (y=+) + Vec((-1.0, 0.3, -0.1)), # on face 6 (x=-) + Vec((0.4, -1.0, 0.2)), # on face 3 (y=-) + Vec((1.0, -0.2, 0.3)), # on face 4 (x=+) + ) + for lf in 1:6 + ξ = opposite_points[lf] + ψ = rt0_hex8_reference_basis(lf, ξ) + n_opp = _RT0_HEX_FACE_NORMALS[lf_opposite[lf]] + @test dot(ψ, n_opp) ≈ 0.0 atol = 1e-14 + end +end