From a469b2e6bba103f7074ab7e6c535ea00a66bbe55 Mon Sep 17 00:00:00 2001 From: JanM Date: Wed, 10 May 2023 11:28:13 +0200 Subject: [PATCH 1/2] introduced interpolation Lagrange{2,RefCube,3} --- src/interpolations.jl | 53 +++++++++++++++++++++++++++++++++++++ test/test_interpolations.jl | 4 ++- 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/src/interpolations.jl b/src/interpolations.jl index 49da85d4ba..c97d27c94b 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -13,6 +13,7 @@ The following interpolations are implemented: * `Lagrange{RefLine,2}` * `Lagrange{RefQuadrilateral,1}` * `Lagrange{RefQuadrilateral,2}` +* `Lagrange{RefQuadrilateral,3}` * `Lagrange{RefTriangle,1}` * `Lagrange{RefTriangle,2}` * `Lagrange{RefTriangle,3}` @@ -491,6 +492,58 @@ function shape_value(ip::Lagrange{RefQuadrilateral, 2}, ξ::Vec{2}, i::Int) throw(ArgumentError("no shape function $i for interpolation $ip")) end +##################################### +# Lagrange RefQuadrilateral order 3 # +##################################### +getnbasefunctions(::Lagrange{RefQuadrilateral, 3}) = 16 + +facedof_indices(::Lagrange{2,RefCube,3}) = ((1,2, 5,6), (2,3, 7,8), (3,4, 9,10), (4,1, 11,12)) +facedof_interior_indices(::Lagrange{2,RefCube,3}) = ((5,6), (7,8), (9,10), (11,12)) +celldof_interior_indices(::Lagrange{2,RefCube,3}) = (13,14,15,16) + +function reference_coordinates(::Lagrange{2,RefCube,3}) + return [Vec{2, Float64}((-1.0, -1.0)), + Vec{2, Float64}(( 1.0, -1.0)), + Vec{2, Float64}(( 1.0, 1.0)), + Vec{2, Float64}((-1.0, 1.0)), + Vec{2, Float64}((-1/3, -1.0)), + Vec{2, Float64}(( 1/3, -1.0)), + Vec{2, Float64}(( 1.0, -1/3)), + Vec{2, Float64}(( 1.0, 1/3)), + Vec{2, Float64}(( 1/3, 1.0)), + Vec{2, Float64}((-1/3, 1.0)), + Vec{2, Float64}((-1.0, 1/3)), + Vec{2, Float64}((-1.0, -1/3)), + Vec{2, Float64}((-1/3, -1/3)), + Vec{2, Float64}(( 1/3, -1/3)), + Vec{2, Float64}((-1/3, 1/3)), + Vec{2, Float64}(( 1/3, 1/3))] +end + +function shape_value(ip::Lagrange{RefQuadrilateral, 3}, ξ::Vec{2}, i::Int) + # See https://defelement.com/elements/examples/quadrilateral-Q-3.html + # Transform domain from [-1, 1] × [-1, 1] to [0, 1] × [0, 1] + ξ_x = ξ[1]*0.5 + 0.5 + ξ_y = ξ[2]*0.5 + 0.5 + i == 1 && return (81*ξ_x^3*ξ_y^3)/4 - (81*ξ_x^3*ξ_y^2)/2 + (99*ξ_x^3*ξ_y)/4 - (9*ξ_x^3)/2 - (81*ξ_x^2*ξ_y^3)/2 + (81*ξ_x^2*ξ_y^2) - (99*ξ_x^2*ξ_y)/2 + (9*ξ_x^2) + (99*ξ_x*ξ_y^3)/4 - (99*ξ_x*ξ_y^2)/2 + (121*ξ_x*ξ_y)/4 - (11*ξ_x)/2 - (9*ξ_y^3)/2 + 9*ξ_y^2 - (11*ξ_y)/2 + 1 + i == 2 && return (ξ_x*( - 81*ξ_x^2*ξ_y^3 + 162*ξ_x^2*ξ_y^2 - 99*ξ_x^2*ξ_y + 18*ξ_x^2 + 81*ξ_x*ξ_y^3 - 162*ξ_x*ξ_y^2 + 99*ξ_x*ξ_y - 18*ξ_x - 18*ξ_y^3 + 36*ξ_y^2 - 22*ξ_y + 4))/4 + i == 4 && return (ξ_y*( - 81*ξ_x^3*ξ_y^2 + 81*ξ_x^3*ξ_y - 18*ξ_x^3 + 162*ξ_x^2*ξ_y^2 - 162*ξ_x^2*ξ_y + 36*ξ_x^2 - 99*ξ_x*ξ_y^2 + 99*ξ_x*ξ_y - 22*ξ_x + 18*ξ_y^2 - 18*ξ_y + 4))/4 + i == 3 && return (ξ_x*ξ_y*(81*ξ_x^2*ξ_y^2 - 81*ξ_x^2*ξ_y + 18*ξ_x^2 - 81*ξ_x*ξ_y^2 + 81*ξ_x*ξ_y - 18*ξ_x + 18*ξ_y^2 - 18*ξ_y + 4))/4 + i == 5 && return (9*ξ_x*( - 27*ξ_x^2*ξ_y^3 + 54*ξ_x^2*ξ_y^2 - 33*ξ_x^2*ξ_y + 6*ξ_x^2 + 45*ξ_x*ξ_y^3 - 90*ξ_x*ξ_y^2 + 55*ξ_x*ξ_y - 10*ξ_x - 18*ξ_y^3 + 36*ξ_y^2 - 22*ξ_y + 4))/4 + i == 6 && return (9*ξ_x*(27*ξ_x^2*ξ_y^3 - 54*ξ_x^2*ξ_y^2 + 33*ξ_x^2*ξ_y - 6*ξ_x^2 - 36*ξ_x*ξ_y^3 + 72*ξ_x*ξ_y^2 - 44*ξ_x*ξ_y + 8*ξ_x + 9*ξ_y^3 - 18*ξ_y^2 + 11*ξ_y - 2))/4 + i == 12 && return (9*ξ_y*( - 27*ξ_x^3*ξ_y^2 + 45*ξ_x^3*ξ_y - 18*ξ_x^3 + 54*ξ_x^2*ξ_y^2 - 90*ξ_x^2*ξ_y + 36*ξ_x^2 - 33*ξ_x*ξ_y^2 + 55*ξ_x*ξ_y - 22*ξ_x + 6*ξ_y^2 - 10*ξ_y + 4))/4 + i == 11 && return (9*ξ_y*(27*ξ_x^3*ξ_y^2 - 36*ξ_x^3*ξ_y + 9*ξ_x^3 - 54*ξ_x^2*ξ_y^2 + 72*ξ_x^2*ξ_y - 18*ξ_x^2 + 33*ξ_x*ξ_y^2 - 44*ξ_x*ξ_y + 11*ξ_x - 6*ξ_y^2 + 8*ξ_y - 2))/4 + i == 7 && return (9*ξ_x*ξ_y*(27*ξ_x^2*ξ_y^2 - 45*ξ_x^2*ξ_y + 18*ξ_x^2 - 27*ξ_x*ξ_y^2 + 45*ξ_x*ξ_y - 18*ξ_x + 6*ξ_y^2 - 10*ξ_y + 4))/4 + i == 8 && return (9*ξ_x*ξ_y*( - 27*ξ_x^2*ξ_y^2 + 36*ξ_x^2*ξ_y - 9*ξ_x^2 + 27*ξ_x*ξ_y^2 - 36*ξ_x*ξ_y + 9*ξ_x - 6*ξ_y^2 + 8*ξ_y - 2))/4 + i == 10 && return (9*ξ_x*ξ_y*(27*ξ_x^2*ξ_y^2 - 27*ξ_x^2*ξ_y + 6*ξ_x^2 - 45*ξ_x*ξ_y^2 + 45*ξ_x*ξ_y - 10*ξ_x + 18*ξ_y^2 - 18*ξ_y + 4))/4 + i == 9 && return (9*ξ_x*ξ_y*( - 27*ξ_x^2*ξ_y^2 + 27*ξ_x^2*ξ_y - 6*ξ_x^2 + 36*ξ_x*ξ_y^2 - 36*ξ_x*ξ_y + 8*ξ_x - 9*ξ_y^2 + 9*ξ_y - 2))/4 + i == 13 && return (81*ξ_x*ξ_y*(9*ξ_x^2*ξ_y^2 - 15*ξ_x^2*ξ_y + 6*ξ_x^2 - 15*ξ_x*ξ_y^2 + 25*ξ_x*ξ_y - 10*ξ_x + 6*ξ_y^2 - 10*ξ_y + 4))/4 + i == 14 && return (81*ξ_x*ξ_y*( - 9*ξ_x^2*ξ_y^2 + 15*ξ_x^2*ξ_y - 6*ξ_x^2 + 12*ξ_x*ξ_y^2 - 20*ξ_x*ξ_y + 8*ξ_x - 3*ξ_y^2 + 5*ξ_y - 2))/4 + i == 15 && return (81*ξ_x*ξ_y*( - 9*ξ_x^2*ξ_y^2 + 12*ξ_x^2*ξ_y - 3*ξ_x^2 + 15*ξ_x*ξ_y^2 - 20*ξ_x*ξ_y + 5*ξ_x - 6*ξ_y^2 + 8*ξ_y - 2))/4 + i == 16 && return (81*ξ_x*ξ_y*(9*ξ_x^2*ξ_y^2 - 12*ξ_x^2*ξ_y + 3*ξ_x^2 - 12*ξ_x*ξ_y^2 + 16*ξ_x*ξ_y - 4*ξ_x + 3*ξ_y^2 - 4*ξ_y + 1))/4 + throw(ArgumentError("no shape function $i for interpolation $ip")) +end + ################################ # Lagrange RefTriangle order 1 # ################################ diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 08c9215ae5..c759513fb0 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -5,6 +5,7 @@ Lagrange{RefLine, 2}(), Lagrange{RefQuadrilateral, 1}(), Lagrange{RefQuadrilateral, 2}(), + Lagrange{RefQuadrilateral, 3}(), Lagrange{RefTriangle, 1}(), Lagrange{RefTriangle, 2}(), Lagrange{RefTriangle, 3}(), @@ -102,7 +103,8 @@ if k == dof @test N_dof ≈ 1.0 else - @test N_dof ≈ 0.0 atol=4eps(Float64) #broken=typeof(interpolation)==Lagrange{2, RefTetrahedron, 5}&&dof==4&&k==18 + factor = interpolation isa Lagrange{RefQuadrilateral, 3} ? 200 : 4 + @test N_dof ≈ 0.0 atol = factor * eps(Float64) end end end From dda353d1a92c11547a576595cf27396dbc837f1f Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Fri, 26 May 2023 10:12:44 +0200 Subject: [PATCH 2/2] this got lost --- src/interpolations.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/interpolations.jl b/src/interpolations.jl index c97d27c94b..c6d0606580 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -497,11 +497,11 @@ end ##################################### getnbasefunctions(::Lagrange{RefQuadrilateral, 3}) = 16 -facedof_indices(::Lagrange{2,RefCube,3}) = ((1,2, 5,6), (2,3, 7,8), (3,4, 9,10), (4,1, 11,12)) -facedof_interior_indices(::Lagrange{2,RefCube,3}) = ((5,6), (7,8), (9,10), (11,12)) -celldof_interior_indices(::Lagrange{2,RefCube,3}) = (13,14,15,16) +facedof_indices(::Lagrange{RefQuadrilateral, 3}) = ((1,2, 5,6), (2,3, 7,8), (3,4, 9,10), (4,1, 11,12)) +facedof_interior_indices(::Lagrange{RefQuadrilateral, 3}) = ((5,6), (7,8), (9,10), (11,12)) +celldof_interior_indices(::Lagrange{RefQuadrilateral, 3}) = (13,14,15,16) -function reference_coordinates(::Lagrange{2,RefCube,3}) +function reference_coordinates(::Lagrange{RefQuadrilateral, 3}) return [Vec{2, Float64}((-1.0, -1.0)), Vec{2, Float64}(( 1.0, -1.0)), Vec{2, Float64}(( 1.0, 1.0)),