resolve inconsistency between refcube and refshape interpolation in 3D

[koehlerson]

#520

[codecov-commenter]

Codecov Report

Base: 92.37% // Head: 92.37% // No change to project coverage 👍

Coverage data is based on head (6cb9ac4) compared to base (084a7c9).
Patch coverage: 0.00% of modified lines in pull request are covered.

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #523   +/-   ##
=======================================
  Coverage   92.37%   92.37%           
=======================================
  Files          22       22           
  Lines        3776     3776           
=======================================
  Hits         3488     3488           
  Misses        288      288           

Impacted Files	Coverage Δ
src/interpolations.jl	89.16% <0.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[termi-official]

Can we include some regression test to catch future problems with wrong dispatches on these functions?

[koehlerson]

second order Hexahedron has inconsistent edges as well :(

highorderedge = (1, 2, 9)
firstorderedge = (1, 2)
highorderedge = (2, 3, 10)
firstorderedge = (2, 3)
highorderedge = (3, 4, 11)
firstorderedge = (3, 4)
highorderedge = (4, 1, 12)
firstorderedge = (4, 1)
highorderedge = (5, 6, 13)
firstorderedge = (1, 5)
highorderedge = (6, 7, 14)
firstorderedge = (2, 6)
highorderedge = (7, 8, 15)
firstorderedge = (3, 7)
highorderedge = (8, 5, 16)
firstorderedge = (4, 8)
highorderedge = (1, 5, 17)
firstorderedge = (5, 6)
highorderedge = (2, 6, 18)
firstorderedge = (6, 7)
highorderedge = (3, 7, 19)
firstorderedge = (7, 8)
highorderedge = (4, 8, 20)
firstorderedge = (8, 5)

I guess the points match but are numbered differently

[koehlerson]

is there a better way to extract the type of the interpolation, e.g. Lagrange, Serendipity etc.?

[fredrikekre]

You can check interpolation isa Lagrange instead

[koehlerson]

ideally this should be checked for non-lagrangian interpolations, too. Since in the codebase DiscontinuousLagrange, CrouzeixRaviart, BubbleEnrichedLagrange is already introduced

[termi-official]

Also I will include H(div) and H(curl) interpolations soon, because our group need them.

[KnutAM]

is there a better way to extract the type of the interpolation, e.g. Lagrange, Serendipity etc.?

@koehlerson : I noted down this trick with ref to this thread (used it in testing as well), but recently I saw

@inline basetypeof(x::T) where T = basetypeof(T)
@generated function basetypeof(::Type{T}) where T
    if T isa Union
        T
    else
        getfield(parentmodule(T), nameof(T))
    end
end

from https://github.com/rafaqz/DimensionalData.jl/blob/3087e3e445ddb57e3f02630d6077cff585a132ac/src/LookupArrays/utils.jl#L73-L80

[koehlerson]

I needed to change the edge ordering and IIRC @lijas you mentioned something that the order is different in vtk. Do we need to change something additional for the vtk export?

[lijas]

I dont remember where I noticed that... but for hexahedron it looks like we use the same numbering
https://vtk.org/doc/nightly/html/classvtkQuadraticHexahedron.html

[koehlerson]

I dont remember where I noticed that... but for hexahedron it looks like we use the same numbering https://vtk.org/doc/nightly/html/classvtkQuadraticHexahedron.html

I had this conversation in mind: #353 (comment)

[lijas]

Actually, it seems like they order their edges for hexahedron differently:
VTK: (0,1), (1,2), (2,3), (3,0), (4,5), (5,6), (6,7), (7,4), (0,4), (1,5), (2,6), (3,7) (0-based indexing)
Ferrite. (1,2), (2,3), (3,4), (4,1), (1,5), (2,6), (3,7), (4,8), (5,6), (6,7), (7,8), (8,5)

So they number their vertical edges 9,10,11, and 12, and for us they are 5,6,7 and 8

[koehlerson]

can you tell me where we need to change code in the vtk.jl file? I'm also fine if you want to push the changes on this branch in case you find some time in the coming days :)

Edit: hm does it actually matter when we export it for normal stuff (H1 problems, no Hcurl or Hdiv, nodal interpolation ansatz)? Because nodes_to_vtkorder should still be identity mapping, shouldn't it?

[lijas]

It should matter if we output the results on a quadratic hexahedron mesh, right?

I guess we should change the ordering of the edges
here:

Ferrite.jl/src/Grid/grid.jl

Line 732 in f3057f3

edges(c::Union{Hexahedron,Cell{3,20,6}}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]), (c.nodes[5],c.nodes[6]), (c.nodes[6],c.nodes[7]), (c.nodes[7],c.nodes[8]), (c.nodes[8],c.nodes[5]), (c.nodes[1],c.nodes[5]), (c.nodes[2],c.nodes[6]), (c.nodes[3],c.nodes[7]), (c.nodes[4],c.nodes[8]))

and here:

Ferrite.jl/src/interpolations.jl

Line 501 in f3057f3

edges(::Lagrange{3,RefCube,1}) = ((1,2), (2,3), (3,4), (4,1), (1,5), (2,6), (3,7), (4,8), (5,6), (6,7), (7,8), (8,5))

to be consistent with that of VTK.

Or we can implement nodes_to_vtkorder for QuadraticHexahedrons.

[koehlerson]

ye but nodes_to_vtkorder only gives the nodes in the vtk order and not the edges. So, we can still pass a hexahedron with matching node order but with incorrect edge order? I'm also okay to revert the last commit and change the serendipity edge and linear hexahedron edge definition the be align with vtk. But of course that is somehow breaking if anyone uses super extensively linear hexahedron and some wild stuff w.r.t. their edges.

What do you think @fredrikekre @KristofferC @kimauth ?

[termi-official]

We really should decouple the vtk representation from our internal representation. If I find some time I will make a better VTK export based on interpolation of data to the nodes instead of some reordering magic (which does not work e.g. for CR elements).

[lijas]

So, we can still pass a hexahedron with matching node order but with incorrect edge order?

Yes, this is true, I did not think of that. So our edge ordering can be different I guess 🤔

[koehlerson]

So, we can still pass a hexahedron with matching node order but with incorrect edge order?

Yes, this is true, I did not think of that. So our edge ordering can be different I guess thinking

Actually it shouldn't matter that we order them differently, because apparently all of us use the whole time the linear hexahedron export and the linear hexahedron shares the "wrong" edge numbering (just without the higher order nodes). I just aligned the full quadratic interpolation (note that I don't mean Cell{3,20,6} or the 3D Serendipity but rather the Lagrange{3,RefCube,2} stuff) to the linear definition