Use generated functions for edgedof_indices and facedof_indices

[KnutAM]

This removes the need to manually implement these.

While tests pass, the following changes occur for facedof_indices:

# Nedelec{RefTetrahedron, 1, 3}
((1, 2, 3), (1, 4, 5), (2, 5, 6), (3, 4, 6)) # master
((3, 2, 1), (1, 5, 4), (2, 6, 5), (4, 6, 3)) # pr
# Nedelec{RefHexahedron, 1, 3}
((1, 2, 3, 4), (1, 5, 9, 10), (2, 6, 10, 11), (3, 7, 11, 12), (4, 8, 9, 12), (5, 6, 7, 8)) # master
((4, 3, 2, 1), (1, 10, 5, 9), (2, 11, 6, 10), (3, 12, 7, 11), (9, 8, 12, 4), (5, 6, 7, 8)) # pr
# Lagrange{RefPyramid, 2}
((1, 3, 4, 2, 7, 11, 9, 6, 14), (1, 2, 5, 6, 10, 8), (1, 5, 3, 7, 12, 8), (2, 4, 5, 9, 13, 10), (3, 5, 4, 12, 13, 11)) # master
((1, 3, 4, 2, 7, 11, 9, 6, 14), (1, 2, 5, 6, 10, 8), (1, 5, 3, 8, 12, 7), (2, 4, 5, 9, 13, 10), (3, 5, 4, 12, 13, 11)) #pr

The numbering scheme as not been fully clear for facedof_indices, and with this PR we also ensure that it is consistent between interpolations.

test script

using Ferrite
using Serialization

function store_results(filename)
    Hcurl_interpolations = [
        Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(), Nedelec{RefQuadrilateral, 1}(),
        Nedelec{RefTetrahedron, 1}(), Nedelec{RefHexahedron, 1}(),
        ]
    Hdiv_interpolations = [
        RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), RaviartThomas{RefQuadrilateral, 1}(),
        RaviartThomas{RefTetrahedron, 1}(), RaviartThomas{RefHexahedron, 1}(),
        BrezziDouglasMarini{RefTriangle, 1}(),
        ]
    regular_interpolations = Any[Lagrange{RefLine, 1}(),
                Lagrange{RefLine, 2}(),
                Lagrange{RefQuadrilateral, 1}(),
                Lagrange{RefQuadrilateral, 2}(),
                Lagrange{RefQuadrilateral, 3}(),
                Lagrange{RefTriangle, 1}(),
                Lagrange{RefTriangle, 2}(),
                Lagrange{RefTriangle, 3}(),
                Lagrange{RefTriangle, 4}(),
                Lagrange{RefTriangle, 5}(),
                Lagrange{RefHexahedron, 1}(),
                Lagrange{RefHexahedron, 2}(),
                Serendipity{RefQuadrilateral, 2}(),
                Serendipity{RefHexahedron, 2}(),
                Lagrange{RefTetrahedron, 1}(),
                Lagrange{RefTetrahedron, 2}(),
                Lagrange{RefPrism, 1}(),
                Lagrange{RefPrism, 2}(),
                Lagrange{RefPyramid, 1}(),
                Lagrange{RefPyramid, 2}(),
                #
                DiscontinuousLagrange{RefLine, 0}(),
                DiscontinuousLagrange{RefQuadrilateral, 0}(),
                DiscontinuousLagrange{RefHexahedron, 0}(),
                DiscontinuousLagrange{RefTriangle, 0}(),
                DiscontinuousLagrange{RefTetrahedron, 0}(),
                DiscontinuousLagrange{RefLine, 1}(),
                DiscontinuousLagrange{RefQuadrilateral, 1}(),
                DiscontinuousLagrange{RefHexahedron, 1}(),
                DiscontinuousLagrange{RefTriangle, 1}(),
                DiscontinuousLagrange{RefTetrahedron, 1}(),
                DiscontinuousLagrange{RefPrism, 1}(),
                DiscontinuousLagrange{RefPyramid, 1}(),
                #
                BubbleEnrichedLagrange{RefTriangle, 1}(),
                #
                CrouzeixRaviart{RefTriangle, 1}(),
                CrouzeixRaviart{RefTetrahedron, 1}(),
                RannacherTurek{RefQuadrilateral, 1}(),
                RannacherTurek{RefHexahedron, 1}(),
            ]

    ips = append!(regular_interpolations, Hdiv_interpolations, Hcurl_interpolations)

    println("edgedof_indices")
    edi = Dict{String,Any}()
    fdi = Dict{String,Any}()
    for ip in ips
        key = string(typeof(ip))
        edi[key] = Ferrite.edgedof_indices(ip)
        fdi[key] = Ferrite.facedof_indices(ip)
    end
    serialize(filename, (edi, fdi))
end

# store_results("pr.bin")
# store_results("master.bin")

pr_e, pr_f = deserialize("pr.bin")
ma_e, ma_f = deserialize("master.bin")

for k in keys(pr_e)
    if pr_e[k] !== ma_e[k]
        println("edges for $k don't match (master/pr):")
        println(ma_e[k])
        println(pr_e[k])
    end
    if pr_f[k] !== ma_f[k]
        println("faces for $k don't match (master/pr):")
        println(ma_f[k])
        println(pr_f[k])
    end
end

#=
faces for Nedelec{RefTetrahedron, 1, 3} don't match (master/pr):
((1, 2, 3), (1, 4, 5), (2, 5, 6), (3, 4, 6))
((3, 2, 1), (1, 5, 4), (2, 6, 5), (4, 6, 3))
faces for Nedelec{RefHexahedron, 1, 3} don't match (master/pr):
((1, 2, 3, 4), (1, 5, 9, 10), (2, 6, 10, 11), (3, 7, 11, 12), (4, 8, 9, 12), (5, 6, 7, 8))
((4, 3, 2, 1), (1, 10, 5, 9), (2, 11, 6, 10), (3, 12, 7, 11), (9, 8, 12, 4), (5, 6, 7, 8))
faces for Lagrange{RefPyramid, 2} don't match (master/pr):
((1, 3, 4, 2, 7, 11, 9, 6, 14), (1, 2, 5, 6, 10, 8), (1, 5, 3, 7, 12, 8), (2, 4, 5, 9, 13, 10), (3, 5, 4, 12, 13, 11))
((1, 3, 4, 2, 7, 11, 9, 6, 14), (1, 2, 5, 6, 10, 8), (1, 5, 3, 8, 12, 7), (2, 4, 5, 9, 13, 10), (3, 5, 4, 12, 13, 11))
=#

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 94.30%. Comparing base (efdf74d) to head (7d6a514).
⚠️ Report is 8 commits behind head on master.

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1293   +/-   ##
=======================================
  Coverage   94.30%   94.30%           
=======================================
  Files          40       40           
  Lines        6762     6748   -14     
=======================================
- Hits         6377     6364   -13     
+ Misses        385      384    -1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[lijas]

Nice! Do we want the order of the indices in such a way that cross(dxdxi, dxdeta) (on the face) is a outward pointing normal? It seems like we test something like that here.

Also, I guess if we re-implement n_indices_on_face and(::Interpolation) n_indices_on_edge(::Interpolation) which we once had but deleted, we can also generate the functions facedofs_interior_indeices and edgedofs_interior_indices.

[KnutAM]

Nice! Do we want the order of the indices in such a way that cross(dxdxi, dxdeta) (on the face) is a outward pointing normal? It seems like we test something like that here.

Currently, we don't have the mapping (AFAIK) for reducing from a facet of a dim cell/interpolation to the dim-1 cell/interpolation on that facet. However, the ordering of the vertices defining the reference cell fullfils this and is are defined by the reference_faces here.

Ferrite.jl/src/Grid/grid.jl

Line 214 in 236eb50

reference_faces(::Type{RefTetrahedron}) = ((1, 3, 2), (1, 2, 4), (2, 3, 4), (1, 4, 3)) # f1 ... f4

Also, I guess if we re-implement n_indices_on_face and(::Interpolation) n_indices_on_edge(::Interpolation) which we once had but deleted, we can also generate the functions facedofs_interior_indeices and edgedofs_interior_indices.

That would assume that we follow the current ordering scheme which would contradict e.g. #1188 I think (CC: @termi-official) since this has to be aligned with how the reference_shape_value is defined.

[termi-official]

Do we want the order of the indices in such a way that cross(dxdxi, dxdeta) (on the face) is a outward pointing normal?

I initially ordered them in a way that we can map between codimensional entities back and forth. So if your hex' faces are ordered in such that your normal is outwards pointing, then the faces in terms of assigned quads preserve this property.

[...] we can also generate the functions facedofs_interior_indeices and edgedofs_interior_indices.

How?

[KnutAM]

[...] we can also generate the functions facedofs_interior_indeices and edgedofs_interior_indices.

How?

I guess something like

function edgedof_interior_indices(ip)
    n = sum(vdofs -> sum(vdofs; init=0), vertexdof_indices(ip))
    dofs = Any[]
    for i = 1:nedges(ip)
        push!(dofs, ntuple(i -> n + i, n_indices_on_edge(ip)))
        n += n_indices_on_edge(ip)
    end
    return dofs
end

but as generated ? But this assumes dof-numbering vertex-edge-face-volume and I think it is a bit problematic that reference_shape_value would have to follow this implicitly then. On the other hand, this could help for implementation to make it all consistent and the same, and with the suitable reference_coordinates implementation we can test that we have the correct numbering.

Probably be shouldn't go overboard with the generation, but we could of course probably define the shape functions on a per-entity basis following the vertexdof_indices, edgedof_interior_indices, facedof_interior_indices, and volumedof_interior_indices and based on this generate the shape_reference_value function, but that is a significant refactor that risk obscuring the code a lot.

[termi-official]

Not sure if I would do that tbh. I really like the addition to always generate these two functions, because it is just error prone to sync with the functions uses during dof distribution. Especially since we want some very specific convention to be follows.

Instead of adding back n_indices_on_X I would rather suggest that we move forward, as also suggested by you here #1255 (comment) , and just utils to generate the full element with everything from the functionals on some ref shape (and support nodes), e.g. with some macro. For me it does not make much sense why we should add this, but leave the remaining part manual then.

[KnutAM]

Instead of adding back n_indices_on_X I would rather suggest that we move forward, as also suggested by you here #1255 (comment) , and just utils ...

Yes, I'm not suggesting that we should do that here. This PR is only to add the generated function for the non-interior indices on edges and faces, since these are pure duplication of information already provided, just in a different format.

Is it the current content of this PR that you don't like, or is it the addition in the comment above?

[termi-official]

Is it the current content of this PR that you don't like, or is it the addition in the comment above?

The latter, because it is quite arbitrary.

[termi-official]

@kimauth @DRollin any idea why this PR causes FerriteInterfaceElements to fail? I do not see these functions being used in there.

[KnutAM]

@kimauth @DRollin any idea why this PR causes FerriteInterfaceElements to fail? I do not see these functions being used in there.

Have tracked this down: Ferrite-FEM/FerriteInterfaceElements.jl#32

[termi-official]

As I do not think this is an urgent feature for someone to have immediately, let us wait until downstream has a fix and then merge.

[KnutAM]

let us wait until downstream has a fix and then merge.

Downstream fixed

[termi-official]

The Lagrange and Nedelec implementations were genuine wrong by not respecting the local orderings correctly. Other than that this should be a nice and helpful simplification!

Let us postpone any discussion on the convention for these orderings until we have cases where these orderings actually matter. Then we can also add suitable tests.

[KnutAM]

The Lagrange and Nedelec implementations were genuine wrong by not respecting the local orderings correctly

Where is that specified (asking since I obviously missed that when adding the Nedelec)?