update WIP PR extend Lagrange dim 2

src/interpolations.jl

@@ -330,6 +330,86 @@ function value(ip::Lagrange{2,RefTetrahedron,2}, i::Int, ξ::Vec{2})
     i == 6 && return 4ξ_x * γ
     throw(ArgumentError("no shape function $i for interpolation $ip"))
 end
+###############################################
+# Lagrange dim 2 RefTetrahedron order 3, 4, 5 #
+###############################################
+# see https://getfem.readthedocs.io/en/latest/userdoc/appendixA.html
+
+function Ferrite.getnbasefunctions(::Lagrange{2,RefTetrahedron,order}) where order 
+    return (order + 1) * (order + 2) ÷ 2
+end
+
+Ferrite.nvertexdofs(::Lagrange{2,RefTetrahedron,order}) where order = 1
+Ferrite.nedgedofs(::Lagrange{2,RefTetrahedron,order}) where order = order - 1
+Ferrite.nfacedofs(::Lagrange{2,RefTetrahedron,order}) where order = order - 1
+Ferrite.ncelldofs(::Lagrange{2,RefTetrahedron,order}) where order = (order + 1) * (order + 2) ÷ 2 - 3 * order
+
+# permutation to switch to Ferrite ordering
+const permdof2D = Dict{Int, Vector{Int}}()
+permdof2D[1] = [1, 2, 3]
+permdof2D[2] = [3, 6, 1, 5, 4, 2]
+permdof2D[3] = [4, 10, 1, 7, 9, 8, 5, 2, 3, 6] 
+permdof2D[4] = [5, 15, 1, 9, 12, 14, 13, 10, 6, 2, 3, 4, 7, 8, 11] 
+permdof2D[5] = [6, 21, 1, 11, 15, 18, 20, 19, 16, 12, 7, 2, 3, 4, 5, 8, 9, 10,
+                13, 14, 17] 
+
+function Ferrite.vertices(::Lagrange{2,RefTetrahedron,order}) where order 
+    return (1, 2, 3)
+end
+
+function Ferrite.faces(::Lagrange{2,RefTetrahedron,order}) where order     
+    if order == 1
+        return ((1,2), (2,3), (3,1))
+    elseif order == 2 
+        return ((1,2,4), (2,3,5), (3,1,6))
+    elseif order == 3
+        return ((1,2,4,5), (2,3,6,7), (3,1,8,9))
+    elseif order == 4
+        return ((1,2,4,5,6), (2,3,7,8,9), (3,1,10,11,12))
+    elseif order == 5
+        return ((1,2,4,5,6,7), (2,3,8,9,10,11), (3,1,12,13,14,15))
+    else
+        error("Lagrange element of order > 6 not implemented yet")
+    end
+end
+
+function Ferrite.reference_coordinates(::Lagrange{2,RefTetrahedron,order}) where order
+    coordpts = Vector{Ferrite.Vec{2, Float64}}()
+    for k = 0:order
+        for l = 0:(order - k)
+            push!(coordpts, Ferrite.Vec{2, Float64}((l / order, k / order)))
+        end
+    end
+    return coordpts[permdof2D[order], :]
+end
+function Ferrite.value(ip::Lagrange{2,RefTetrahedron,order}, i::Int, 
+                           ξ::Ferrite.Vec{2}) where order
+    i = permdof2D[order][i]
+    ξ_x = ξ[1]
+    ξ_y = ξ[2]
+    γ = 1. - ξ_x - ξ_y
+    i1, i2, i3 = _numlin_basis2D(i, order) 
+    val = one(typeof(γ))
+    i1 ≥ 1 && (val *= prod((order * γ - j) / (j + 1) for j = 0:(i1 - 1))) 
+    i2 ≥ 1 && (val *= prod((order * ξ_x - j) / (j + 1) for j = 0:(i2 - 1)))
+    i3 ≥ 1 && (val *= prod((order * ξ_y - j) / (j + 1) for j = 0:(i3 - 1)))
+    return val
+end
+
+function _numlin_basis2D(i, order) 
+    c, j1, j2, j3 = 0, 0, 0, 0
+    for k = 0:order
+        if i <= c + (order + 1 - k)
+            j2 = i - c - 1
+            break
+        else
+            j3 += 1
+            c += order + 1 - k
+        end
+    end
+    j1 = order - j2 -j3
+    return j1, j2, j3
+end
 
 #########################################
 # Lagrange dim 3 RefTetrahedron order 1 #
@@ -704,7 +784,3 @@ function value(ip::CrouzeixRaviart{2,1}, i::Int, ξ::Vec{2})
     throw(ArgumentError("no shape function $i for interpolation $ip"))
 end
 
-#-------------------------------------------------------------------------------
-# tentative to implement Lagrange element of order k ≥ 3
-include("extend_interpolations.jl")
-#include("extend_interpolations_Ferriteorder.jl")