Extend Lagrange elements for 2D Triangles

CHANGELOG.md

@@ -6,6 +6,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+### Added
+ - New higher order function interpolations for triangles (`Lagrange{2,RefTetrahedron,3}`,
+   `Lagrange{2,RefTetrahedron,4}`, and `Lagrange{2,RefTetrahedron,5}`). ([#482][github-482],
+   [#512][github-512])
 ### Changed
  - The default components to constrain in `Dirichlet` and `PeriodicDirichlet` have changed
    from component 1 to all components of the field. For scalar problems this has no effect.
@@ -120,6 +124,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 [github-475]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/475
 [github-478]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/478
 [github-481]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/481
+[github-482]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/482
 [github-487]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/487
 [github-494]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/494
 [github-496]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/496
@@ -129,6 +134,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 [github-505]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/505
 [github-506]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/506
 [github-509]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/509
+[github-512]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/512
 
 [Unreleased]: https://github.com/Ferrite-FEM/Ferrite.jl/compare/v0.3.8...HEAD
 [0.3.8]: https://github.com/Ferrite-FEM/Ferrite.jl/compare/v0.3.7...v0.3.8

src/interpolations.jl

@@ -330,6 +330,97 @@ 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
+
+const Lagrange2Tri345 = Union{
+    Lagrange{2,RefTetrahedron,3},
+    Lagrange{2,RefTetrahedron,4},
+    Lagrange{2,RefTetrahedron,5},
+}
+
+function getnbasefunctions(ip::Lagrange2Tri345)
+    order = getorder(ip)
+    return (order + 1) * (order + 2) ÷ 2
+end
+
+nvertexdofs(::Lagrange2Tri345) = 1
+nedgedofs(ip::Lagrange2Tri345) = getorder(ip) - 1
+nfacedofs(ip::Lagrange2Tri345) = getorder(ip) - 1
+function ncelldofs(ip::Lagrange2Tri345)
+    order = getorder(ip)
+    return (order + 1) * (order + 2) ÷ 2 - 3 * order
+end
+
+# Permutation to switch numbering to Ferrite ordering
+const permdof2D = Dict{Int,Vector{Int}}(
+    1 => [1, 2, 3],
+    2 => [3, 6, 1, 5, 4, 2],
+    3 => [4, 10, 1, 7, 9, 8, 5, 2, 3, 6],
+    4 => [5, 15, 1, 9, 12, 14, 13, 10, 6, 2, 3, 4, 7, 8, 11],
+    5 => [6, 21, 1, 11, 15, 18, 20, 19, 16, 12, 7, 2, 3, 4, 5, 8, 9, 10, 13, 14, 17],
+)
+
+vertices(::Lagrange2Tri345) = (1, 2, 3)
+
+function faces(ip::Lagrange2Tri345)
+    order = getorder(ip)
+    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("unreachable")
+    end
+end
+
+function reference_coordinates(ip::Lagrange2Tri345)
+    order = getorder(ip)
+    coordpts = Vector{Vec{2, Float64}}()
+    for k = 0:order
+        for l = 0:(order - k)
+            push!(coordpts, Vec{2, Float64}((l / order, k / order)))
+        end
+    end
+    return permute!(coordpts, permdof2D[order])
+end
+
+function value(ip::Lagrange2Tri345, i::Int, ξ::Vec{2})
+    order = getorder(ip)
+    i = permdof2D[order][i]
+    ξ_x = ξ[1]
+    ξ_y = ξ[2]
+    γ = 1. - ξ_x - ξ_y
+    i1, i2, i3 = _numlin_basis2D(i, order)
+    val = one(γ)
+    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 #

test/test_cellvalues.jl

@@ -6,6 +6,9 @@ for (func_interpol, quad_rule) in  (
                                     (Lagrange{2, RefCube, 2}(), QuadratureRule{2, RefCube}(2)),
                                     (Lagrange{2, RefTetrahedron, 1}(), QuadratureRule{2, RefTetrahedron}(2)),
                                     (Lagrange{2, RefTetrahedron, 2}(), QuadratureRule{2, RefTetrahedron}(2)),
+                                    (Lagrange{2, RefTetrahedron, 3}(), QuadratureRule{2, RefTetrahedron}(2)),
+                                    (Lagrange{2, RefTetrahedron, 4}(), QuadratureRule{2, RefTetrahedron}(2)),
+                                    (Lagrange{2, RefTetrahedron, 5}(), QuadratureRule{2, RefTetrahedron}(2)),
                                     (Lagrange{3, RefCube, 1}(), QuadratureRule{3, RefCube}(2)),
                                     (Serendipity{2, RefCube, 2}(), QuadratureRule{2, RefCube}(2)),
                                     (Lagrange{3, RefTetrahedron, 1}(), QuadratureRule{3, RefTetrahedron}(2)),

test/test_grid_dofhandler_vtk.jl

@@ -425,3 +425,33 @@ end
     # unconnected subset
     test_coloring(generate_grid(Triangle, (10, 10)), union(Set(1:10), Set(70:80)))
 end
+
+@testset "DoF distribution" begin
+    # _________
+    # |\      |
+    # |  \  2 |
+    # | 1  \  |
+    # |______\|
+    grid = generate_grid(Triangle, (1, 1))
+
+    ## Lagrange{2,RefTetrahedron,3}
+    dh = DofHandler(grid)
+    push!(dh, :u, 1, Lagrange{2,RefTetrahedron,3}())
+    close!(dh)
+    @test celldofs(dh, 1) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
+    @test celldofs(dh, 2) == [2, 11, 3, 12, 13, 14, 15, 7, 6, 16]
+
+    ## Lagrange{2,RefTetrahedron,4}
+    dh = DofHandler(grid)
+    push!(dh, :u, 1, Lagrange{2,RefTetrahedron,4}())
+    close!(dh)
+    @test celldofs(dh, 1) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
+    @test celldofs(dh, 2) == [2, 16, 3, 17, 18, 19, 20, 21, 22, 9, 8, 7, 23, 24, 25]
+
+    ## Lagrange{2,RefTetrahedron,5}
+    dh = DofHandler(grid)
+    push!(dh, :u, 1, Lagrange{2,RefTetrahedron,5}())
+    close!(dh)
+    @test celldofs(dh, 1) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]
+    @test celldofs(dh, 2) == [2, 22, 3, 23, 24, 25, 26, 27, 28, 29, 30, 11, 10, 9, 8, 31, 32, 33, 34, 35, 36]
+end

test/test_interpolations.jl

@@ -1,11 +1,15 @@
 @testset "interpolations" begin
 
-for interpolation in (Lagrange{1, RefCube, 1}(),
+for interpolation in (
+                      Lagrange{1, RefCube, 1}(),
                       Lagrange{1, RefCube, 2}(),
                       Lagrange{2, RefCube, 1}(),
                       Lagrange{2, RefCube, 2}(),
                       Lagrange{2, RefTetrahedron, 1}(),
                       Lagrange{2, RefTetrahedron, 2}(),
+                      Lagrange{2, RefTetrahedron, 3}(),
+                      Lagrange{2, RefTetrahedron, 4}(),
+                      Lagrange{2, RefTetrahedron, 5}(),
                       Lagrange{3, RefCube, 1}(),
                       Lagrange{3, RefCube, 2}(),
                       Serendipity{2, RefCube, 2}(),
@@ -21,7 +25,8 @@ for interpolation in (Lagrange{1, RefCube, 1}(),
                       #
                       BubbleEnrichedLagrange{2,RefTetrahedron,1}(),
                       #
-                      CrouzeixRaviart{2,1}(),)
+                      CrouzeixRaviart{2,1}(),
+    )
 
     # Test of utility functions
     ndim = Ferrite.getdim(interpolation)
@@ -47,7 +52,7 @@ for interpolation in (Lagrange{1, RefCube, 1}(),
             if k == node
                 @test N_node[k] ≈ 1.0
             else
-                @test N_node[k] ≈ 0.0
+                @test N_node[k] ≈ 0.0 atol=4eps(Float64)
             end
         end
     end

test/test_utils.jl

@@ -126,6 +126,12 @@ function calculate_volume(::Lagrange{2, RefTetrahedron, 2}, x::Vector{Vec{dim, T
     return vol
 end
 
+# TODO: Only correct for linear sides
+function calculate_volume(::Lagrange{2, RefTetrahedron, O}, x::Vector{Vec{dim, T}}) where {T, dim, O}
+    vol = norm((x[1] - x[3]) × (x[2] - x[3])) * 0.5
+    return vol
+end
+
 function calculate_volume(::Lagrange{3, RefTetrahedron, order}, x::Vector{Vec{3, T}}) where {T, order}
     vol = norm((x[2] - x[1]) ⋅ ((x[3] - x[1]) × (x[4] - x[1]))) / 6.0
     return vol