Refactoring dof distrubtion to allow using generic local celldof orderings

docs/src/devdocs/dofhandler.md

@@ -24,6 +24,6 @@ Ferrite.get_grid
 Ferrite.find_field
 Ferrite._find_field
 Ferrite._close_subdofhandler!
-Ferrite._distribute_dofs_for_cell!
-Ferrite.permute_and_push!
+Ferrite._distribute_field_dofs_for_cell!
+Ferrite.permute_and_set!
 ```

docs/src/devdocs/interpolations.md

@@ -129,13 +129,16 @@ vertex. In this case, we only have a single dof per vertex,
 Ferrite.vertexdof_indices(::QTI) = ((1,), (2,), (3,))
 compare_test(Ferrite.vertexdof_indices) # hide
 ```
-Note that the dof index numbering must be in order, first all in the first
-vertex, then all in the next vertex, and so on. (If we had two indices per
-vertex, we would have `((1,2), (3,4), (5,6))`). Next, we number the edge dofs,
-following the same convention, matching how we defined it above,
+Note that the dofs are assigned in order of increasing codimension, followed by the
+index of the local geometry, consistent with its local orientation. E.g. first all
+dofs of the first vertex, then all dofs of the second vertex, and so on. The dof index
+can be arbitrarily assigned, as long as they are assigned consistently and between 1 and
+the number of basis functions (here 6).
 ```@example InterpolationExample
+#                                           e1    e2    e3
 Ferrite.edgedof_interior_indices(::QTI) = ((4,), (5,), (6,))
 compare_test(Ferrite.edgedof_interior_indices) # hide
+#                                  v1 v2 e1    v2 v3 e2    v3 v1 e3
 Ferrite.edgedof_indices(::QTI) = ((1, 2, 4,), (2, 3, 5,), (3, 1, 6,))
 compare_test(Ferrite.edgedof_indices) # hide
 ```

src/interpolations.jl

@@ -88,21 +88,39 @@ this cache is of the same type no matter the interpolation: the purpose is to ma
 dof-distribution type-stable.
 """
 struct InterpolationInfo
-    nvertexdofs::Vector{Int}
-    nedgedofs::Vector{Int}
-    nfacedofs::Vector{Int}
-    nvolumedofs::Int
+    lvertexdofs::Vector{Int}
+    lvertexdofoffsets::Vector{Int}
+    ledgedofs::Vector{Int}
+    ledgedofoffsets::Vector{Int}
+    lfacedofs::Vector{Int}
+    lfacedofoffsets::Vector{Int}
+    lvolumedofs::Vector{Int}
     reference_dim::Int
     adjust_during_distribution::Bool
     n_copies::Int
     is_discontinuous::Bool
 end
 function InterpolationInfo(interpolation::Interpolation{shape}, n_copies) where {rdim, shape <: AbstractRefShape{rdim}}
+    lvertexdofs = Int[]
+    for ii in vertexdof_indices(interpolation)
+        append!(lvertexdofs, ii)
+    end
+    ledgedofs = Int[]
+    for ii in edgedof_interior_indices(interpolation)
+        append!(ledgedofs, ii)
+    end
+    lfacedofs = Int[]
+    for ii in facedof_interior_indices(interpolation)
+        append!(lfacedofs, ii)
+    end
     info = InterpolationInfo(
-        [length(i) for i in vertexdof_indices(interpolation)],
-        [length(i) for i in edgedof_interior_indices(interpolation)],
-        [length(i) for i in facedof_interior_indices(interpolation)],
-        length(volumedof_interior_indices(interpolation)),
+        lvertexdofs,
+        cumsum([1; [length(i) for i in vertexdof_indices(interpolation)]]),
+        ledgedofs,
+        cumsum([1; [length(i) for i in edgedof_interior_indices(interpolation)]]),
+        lfacedofs,
+        cumsum([1; [length(i) for i in facedof_interior_indices(interpolation)]]),
+        [volumedof_interior_indices(interpolation)...],
         rdim,
         adjust_dofs_during_distribution(interpolation),
         n_copies,

test/integration/test_simple_scalar_convergence.jl

@@ -14,6 +14,7 @@ module ConvergenceTestHelper
     get_geometry(::Ferrite.Interpolation{RefTetrahedron}) = Tetrahedron
     get_geometry(::Ferrite.Interpolation{RefPyramid}) = Pyramid
 
+    get_quadrature_order(::Main.TensorProductQ9TestInterpolation) = 3
     get_quadrature_order(::Lagrange{shape, order}) where {shape, order} = max(2 * order - 1, 2)
     get_quadrature_order(::Lagrange{RefTriangle, 5}) where {shape, order} = 8
     get_quadrature_order(::Lagrange{RefPrism, order}) where {order} = 2 * order # Don't know why
@@ -159,6 +160,7 @@ end # module ConvergenceTestHelper
             Lagrange{RefPyramid, 1}(),
             #
             Serendipity{RefQuadrilateral, 2}(),
+            TensorProductQ9TestInterpolation(),
             Serendipity{RefHexahedron, 2}(),
             #
             BubbleEnrichedLagrange{RefTriangle, 1}(),

test/test_grid_dofhandler_vtk.jl

@@ -723,6 +723,7 @@ end
         grid = generate_grid(Line, (2,))
 
         Ferrite.vertexdof_indices(::VectorLagrangeTest{RefLine, 1, 2}) = ((1, 2), (3, 4))
+        Ferrite.getnbasefunctions(::VectorLagrangeTest{RefLine, 1, 2}) = 4
         dh1 = DofHandler(grid)
         add!(dh1, :u, VectorLagrangeTest{RefLine, 1, 2}())
         close!(dh1)
@@ -733,6 +734,7 @@ end
         @test dh1.cell_dofs == dh2.cell_dofs
 
         Ferrite.vertexdof_indices(::VectorLagrangeTest{RefLine, 1, 3}) = ((1, 2, 3), (4, 5, 6))
+        Ferrite.getnbasefunctions(::VectorLagrangeTest{RefLine, 1, 3}) = 6
         dh1 = DofHandler(grid)
         add!(dh1, :u, VectorLagrangeTest{RefLine, 1, 3}())
         close!(dh1)
@@ -746,6 +748,7 @@ end
     @testset "2d" begin
         grid = generate_grid(Quadrilateral, (2, 2))
         Ferrite.vertexdof_indices(::VectorLagrangeTest{RefQuadrilateral, 1, 2}) = ((1, 2), (3, 4), (5, 6), (7, 8))
+        Ferrite.getnbasefunctions(::VectorLagrangeTest{RefQuadrilateral, 1, 2}) = 8
         dh1 = DofHandler(grid)
         add!(dh1, :u, VectorLagrangeTest{RefQuadrilateral, 1, 2}())
         close!(dh1)
@@ -755,7 +758,7 @@ end
         @test dh1.cell_dofs == dh2.cell_dofs
 
         Ferrite.vertexdof_indices(::VectorLagrangeTest{RefQuadrilateral, 1, 3}) = ((1, 2, 3), (4, 5, 6), (7, 8, 9), (10, 11, 12))
-        Ferrite.facedof_indices(::VectorLagrangeTest{RefQuadrilateral, 1, 3}) = ((1, 2, 3, 4, 5, 6), (4, 5, 6, 7, 8, 9), (7, 8, 9, 10, 11, 12), (10, 11, 12, 1, 2, 3))
+        Ferrite.getnbasefunctions(::VectorLagrangeTest{RefQuadrilateral, 1, 3}) = 12
         dh1 = DofHandler(grid)
         add!(dh1, :u, VectorLagrangeTest{RefQuadrilateral, 1, 3}())
         close!(dh1)

test/test_interpolations.jl

@@ -79,17 +79,15 @@ function _test_interpolation_properties(dofs::NamedTuple, rs::NamedTuple)
     all_dofs = Int[]
     # Vertices numbered first
     append!(all_dofs, collect_all_dofs(dofs.vert))
-    @test all(all_dofs .== 1:length(all_dofs))
+    @test all(1 .≤ all_dofs .≤ dofs.n)
+    @test length(all_dofs) == length(Set(all_dofs))
+
     if rs.rdim ≥ 1 # Test edges
-        # Edges numbered next, no gaps or missing numbers. Sorted by edge number.
+        # Edges have no numbers. Sorted by edge number.
         all_edofs_i = collect_all_dofs(dofs.edge_i)
-        @test all(all_edofs_i .== length(all_dofs) .+ (1:length(all_edofs_i)))
         # - all edge dofs include both vertexdofs and interior edegdofs, and nothing more.
         append!(all_dofs, all_edofs_i)
-        @test all(all_dofs .== 1:length(all_dofs))
         @test length(all_dofs) == length(collect_all_dofs(dofs.vert)) + length(all_edofs_i)
-        # Coarse check for C
-        @test Set(collect_all_dofs(dofs.edge)) == Set(1:length(all_dofs))
         # - test each edge individually (Detailed check for C)
         for (edge_nr, edge_vertices) in enumerate(rs.edges)
             vdofs_e = Int[] # dofs.vert for vertices belonging to the current edge
@@ -99,15 +97,14 @@ function _test_interpolation_properties(dofs::NamedTuple, rs::NamedTuple)
             @test Set(dofs.edge[edge_nr]) == Set(vcat(vdofs_e, collect(dofs.edge_i[edge_nr])))
         end
     end
+    @test all(1 .≤ all_dofs .≤ dofs.n)
+    @test length(all_dofs) == length(Set(all_dofs))
+
     if rs.rdim ≥ 2 # Test faces
         # Face numbered next, no gaps or missing numbers. Sorted by face number.
         all_fdofs_i = collect_all_dofs(dofs.face_i)
-        @test all(all_fdofs_i .== length(all_dofs) .+ (1:length(all_fdofs_i)))
         # - all dofs now include vertex dofs, edge dofs and face dofs, but not volume dofs.
         append!(all_dofs, all_fdofs_i)
-        @test all(all_dofs .== 1:length(all_dofs))
-        # Coarse check for C
-        @test Set(collect_all_dofs(dofs.face)) == Set(1:length(all_dofs))
         # - test each face individually (Detailed check for C)
         for (facenr, face_verts) in enumerate(rs.faces)
             vdofs_f = Int[]
@@ -124,12 +121,15 @@ function _test_interpolation_properties(dofs::NamedTuple, rs::NamedTuple)
             @test Set(dofs.face[facenr]) == Set(vcat(vdofs_f, edofs_f, collect(dofs.face_i[facenr])))
         end
     end
+    @test all(1 .≤ all_dofs .≤ dofs.n)
+    @test length(all_dofs) == length(Set(all_dofs))
+
     # Test volume
     # We always test this, since volumedofs are also used by lower-dimensional
     # discontinuous inteprolations to make them internal to the cell, e.g. DiscontinuousLagrange
     # Volumedofs numbered last
     append!(all_dofs, collect(dofs.vol_i))
-    @test all(all_dofs .== 1:length(all_dofs))        # Numbering convention
+    @test all(1 .≤ all_dofs .≤ dofs.n)
 
     # Test D: getnbasefunctions matching number of dof indices
     return @test length(all_dofs) == dofs.n
@@ -140,6 +140,7 @@ end
             Lagrange{RefLine, 2}(),
             Lagrange{RefQuadrilateral, 1}(),
             Lagrange{RefQuadrilateral, 2}(),
+            TensorProductQ9TestInterpolation(),
             Lagrange{RefQuadrilateral, 3}(),
             Lagrange{RefTriangle, 1}(),
             Lagrange{RefTriangle, 2}(),
@@ -254,7 +255,7 @@ end
 
         # regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/520
         interpolation_type = typeof(interpolation).name.wrapper
-        if func_order > 1 && interpolation_type != Ferrite.Serendipity
+        if func_order > 1 && interpolation_type ∉ (Ferrite.Serendipity, TensorProductQ9TestInterpolation)
             first_order = interpolation_type{ref_shape, 1}()
             for (highorderface, firstorderface) in zip(Ferrite.facedof_indices(interpolation), Ferrite.facedof_indices(first_order))
                 for (h_node, f_node) in zip(highorderface, firstorderface)
@@ -546,7 +547,7 @@ end
             # Test functionals associated with the edges
             for edge_nr in 1:Ferrite.nedges(RefShape)
                 edge_coords = getindex.((Ferrite.reference_coordinates(ipg),), Ferrite.reference_edges(RefShape)[edge_nr])
-                tangent = normalize(edge_coords[2] - edge_coords[1])
+                tangent = Tensors.normalize(edge_coords[2] - edge_coords[1])
                 dof_inds = Ferrite.edgedof_interior_indices(ip)[edge_nr]
 
                 ξ(s) = parameterize_edge(edge_coords, s)

test/test_utils.jl

@@ -2,6 +2,45 @@
 
 using Ferrite: reference_shape_value
 
+struct TensorProductQ9TestInterpolation <: Ferrite.ScalarInterpolation{RefQuadrilateral, 2} end
+Ferrite.adjust_dofs_during_distribution(::TensorProductQ9TestInterpolation) = false
+Ferrite.getnbasefunctions(::TensorProductQ9TestInterpolation) = 9
+Ferrite.vertexdof_indices(::TensorProductQ9TestInterpolation) = ((1,), (3,), (9,), (7,))
+Ferrite.edgedof_indices(::TensorProductQ9TestInterpolation) = ((1, 3, 2), (3, 9, 6), (9, 7, 8), (7, 1, 4))
+Ferrite.edgedof_interior_indices(::TensorProductQ9TestInterpolation) = ((2,), (6,), (8,), (4,))
+Ferrite.facedof_indices(ip::TensorProductQ9TestInterpolation) = ((1, 3, 9, 7, 2, 6, 8, 4, 5),)
+Ferrite.facedof_interior_indices(::TensorProductQ9TestInterpolation) = ((5,))
+
+function Ferrite.reference_coordinates(::TensorProductQ9TestInterpolation)
+    return [
+        Vec{2, Float64}((-1.0, -1.0)),
+        Vec{2, Float64}((0.0, -1.0)),
+        Vec{2, Float64}((1.0, -1.0)),
+        Vec{2, Float64}((-1.0, 0.0)),
+        Vec{2, Float64}((0.0, 0.0)),
+        Vec{2, Float64}((1.0, 0.0)),
+        Vec{2, Float64}((-1.0, 1.0)),
+        Vec{2, Float64}((0.0, 1.0)),
+        Vec{2, Float64}((1.0, 1.0)),
+    ]
+end
+
+function Ferrite.reference_shape_value(ip::TensorProductQ9TestInterpolation, ξ::Vec{2}, i::Int)
+    ξ_x = ξ[1]
+    ξ_y = ξ[2]
+    i == 1 && return (ξ_x^2 - ξ_x) * (ξ_y^2 - ξ_y) / 4
+    i == 3 && return (ξ_x^2 + ξ_x) * (ξ_y^2 - ξ_y) / 4
+    i == 9 && return (ξ_x^2 + ξ_x) * (ξ_y^2 + ξ_y) / 4
+    i == 7 && return (ξ_x^2 - ξ_x) * (ξ_y^2 + ξ_y) / 4
+    i == 2 && return (1 - ξ_x^2) * (ξ_y^2 - ξ_y) / 2
+    i == 6 && return (ξ_x^2 + ξ_x) * (1 - ξ_y^2) / 2
+    i == 8 && return (1 - ξ_x^2) * (ξ_y^2 + ξ_y) / 2
+    i == 4 && return (ξ_x^2 - ξ_x) * (1 - ξ_y^2) / 2
+    i == 5 && return (1 - ξ_x^2) * (1 - ξ_y^2)
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+
 #####################################
 # Volume for the reference elements #
 #####################################