Quadratic Serendipity Hex Interpolation

[koehlerson]

Added quadratic hex interpolation, although I'm not 100% sure about the faces dispatch, in particular I'm not sure if my numbering is correct with the new quadratic nodes.

Other than that, I think QuadraticHexahedron should be renamed, because its an incomplete 20 node version and the complete one has 27 nodes. If you agree to this, then I think the new interpolation should go under Serendipity.

Out of curiosity: is there any advantage of the 20 the node version compared to the 27 node one?

[codecov-commenter]

Codecov Report

Merging #353 (719efd1) into master (b73c61c) will increase coverage by 0.51%.
The diff coverage is 98.55%.

@@            Coverage Diff             @@
##           master     #353      +/-   ##
==========================================
+ Coverage   88.19%   88.70%   +0.51%     
==========================================
  Files          20       20              
  Lines        2219     2284      +65     
==========================================
+ Hits         1957     2026      +69     
+ Misses        262      258       -4     

Impacted Files	Coverage Δ
src/interpolations.jl	90.69% <96.96%> (+1.13%)	⬆️
src/Export/VTK.jl	75.71% <100.00%> (+0.35%)	⬆️
src/Grid/grid.jl	86.50% <100.00%> (+0.61%)	⬆️
src/Grid/grid_generators.jl	100.00% <100.00%> (ø)
src/Dofs/MixedDofHandler.jl	84.64% <0.00%> (+1.75%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b73c61c...719efd1. Read the comment docs.

[lijas]

Nice!

VTK calls there 20-node hexahedron "Quadratichexahedron" (https://vtk.org/doc/nightly/html/classvtkQuadraticHexahedron.html)
And their 27-node hexhedron is called "TriQuadraticHexahedron" (https://vtk.org/doc/release/5.10/html/classvtkTriQuadraticHexahedron.html)
(They also have a 24-node hexahedron which we cant handle in Ferrite since we assume the same number of nodes on all edges. Would be cool to add this functionality maybe.).

If I remember correctly, I think we have renumbered our edges differently compared to VTK (we should compare with the numbering in the first link), which would cause output issues in paraview/vtk (but I dont think this is a problem atm since we have no mesh_generator generating 20-node Hexahedrons).

Could be good to add test here aswell?:

Ferrite.jl/test/test_grid_dofhandler_vtk.jl

Line 12 in b73c61c

@testset "Grid, DofHandler, vtk" begin

[koehlerson]

VTK calls there 20-node hexahedron "Quadratichexahedron" (https://vtk.org/doc/nightly/html/classvtkQuadraticHexahedron.html)
And their 27-node hexhedron is called "TriQuadraticHexahedron" (https://vtk.org/doc/release/5.10/html/classvtkTriQuadraticHexahedron.html)

OK! Then let's stick to their names. I was just wondering, because I want to use it for some quadratic hex mesh and when I meshed with gmsh, I just noticed that their default seems to be the TriQuadratic one, while ours now defaults to the 20 node version.

If I remember correctly, I think we have renumbered our edges differently compared to VTK (we should compare with the numbering in the first link), which would cause output issues in paraview/vtk (but I dont think this is a problem atm since we have no mesh_generator generating 20-node Hexahedrons).

Could be good to add test here aswell?:

Ferrite.jl/test/test_grid_dofhandler_vtk.jl

Line 12 in b73c61c

@testset "Grid, DofHandler, vtk" begin

This would be ideal, but for this I'd need a grid_generator as well. I'm currently reading my quadratic hex mesh from gmsh in order to avoid writing the generator :D

W.r.t. the faces is there some place in Ferrite where I can look this up or is this function defining it? As of now e.g. (5,6,7,8,13,14,15,16) corresponds to

8-----15-----7
|            |
16          14
|            |
5-----13-----6

[koehlerson]

I'm asking because

Ferrite.jl/src/interpolations.jl

Line 351 in 8dc3be3

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

defines at least the first face clockwise and for what I've seen in the repo usually it's defined anti clockwise. Maybe it doesn't matter in the end, I'm just confused

[lijas]

The order of the face-dofs is not important anywhere in the code right now (i think). But ideally it should be ordered in the same way as we distribute dofs: Corners/vertices, edges, faces and inner celldofs. So the way you have ordered it now looks good to me.

Side note:
Searching for how VTK order their dofs, I found this: https://blog.kitware.com/wp-content/uploads/2020/03/Implementation-of-rational-Be%CC%81zier-cells-into-VTK-Report.pdf

For full compatibility with vtk, we should follow this, but since we never use such high-order interpolations we dont really have to think about it.

[kimauth]

For testing you could just generate a mesh with a few cells by manually giving the nodes and cells. The MixedDofHandler does something similar in its tests, e.g. here.

[koehlerson]

The order of the face-dofs is not important anywhere in the code right now (i think). But ideally it should be ordered in the same way as we distribute dofs: Corners/vertices, edges, faces and inner celldofs. So the way you have ordered it now looks good to me.

Side note:
Searching for how VTK order their dofs, I found this: https://blog.kitware.com/wp-content/uploads/2020/03/Implementation-of-rational-Be%CC%81zier-cells-into-VTK-Report.pdf

For full compatibility with vtk, we should follow this, but since we never use such high-order interpolations we dont really have to think about it.

Thanks for the nice ressource, I didn't get it in the beginning till I noticed it's cubic stuff. But I, too, think this is aligned with my ordering now

For testing you could just generate a mesh with a few cells by manually giving the nodes and cells. The MixedDofHandler does something similar in its tests, e.g. here.

In @lijas mentioned testset we are looping over the cells and invoke on each cell the grid_generator dispatch, I think we should stay with this testing structure, however, this obviously requires then a new QuadraticHexahedron grid_generator.

I'd say what's left now:

• grid_generator
• vtk export

I tested it locally with a small simulation, seems to work fine :)

[lijas]

Here is a grid generator that seems to work :)

I noticed that our edge order does not follow vtk.
We order them as: bottom horizontal edges -- vertical mid edges -- top horizontal edges
edges(c::Union{Hexahedron,QuadraticHexahedron}) = ((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[1],c.nodes[5]), (c.nodes[2],c.nodes[6]), (c.nodes[3],c.nodes[7]), (c.nodes[4],c.nodes[8]), (c.nodes[5],c.nodes[6]), (c.nodes[6],c.nodes[7]), (c.nodes[7],c.nodes[8]), (c.nodes[8],c.nodes[5]))

but vtk has: bottom horizontal edges -- top horizontal edges -- vertical mid edges
Ferrite.edges(c::Union{Hexahedron,QuadraticHexahedron}) = ((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]))

How is the order of the dofs for your interpolation (on the edges)?

function Ferrite.generate_grid(::Type{QuadraticHexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_z = nel[3]; nel_tot = nel_x*nel_y*nel_z
    nnode_x = 2nel_x + 1; nnode_y = 2nel_y + 1; nnode_z = 2nel_z + 1 #Note: not the actually number of nodes in x/y/z, just a temporary variables

    # Generate nodes
    coords_x = range(left[1], stop=right[1], length=nnode_x)
    coords_y = range(left[2], stop=right[2], length=nnode_y)
    coords_z = range(left[3], stop=right[3], length=nnode_z)
    nodes = Node{3,T}[]

    node_array = fill(0, (nnode_x,nnode_y,nnode_z))
    nodeid = 0
    for k in 1:nnode_z, j in 1:nnode_y, i in 1:nnode_x
        (iseven(i) && iseven(j)) && continue
        (iseven(i) && iseven(k)) && continue
        (iseven(k) && iseven(j)) && continue
        push!(nodes, Node((coords_x[i], coords_y[j], coords_z[k])))
        nodeid += 1
        node_array[i,j,k] = nodeid
    end


    # Generate cells
    cells = QuadraticHexahedron[]
    for k in 1:2:2nel_z, j in 1:2:2nel_y, i in 1:2:2nel_x     
        push!(cells, QuadraticHexahedron((
                node_array[i,j,k], node_array[i+2,j,k], node_array[i+2,j+2,k], node_array[i,j+2,k], # vertices bot 
                node_array[i,j,k+2], node_array[i+2,j,k+2], node_array[i+2,j+2,k+2], node_array[i,j+2,k+2], # vertices top
                node_array[i+1,j,k], node_array[i+2,j+1,k], node_array[i+1,j+2,k], node_array[i,j+1,k], # edges horizontal bottom
                node_array[i+1,j,k+2], node_array[i+2,j+1,k+2], node_array[i+1,j+2,k+2], node_array[i,j+1,k+2], # edges horizontal top
                node_array[i,j,k+1], node_array[i+2,j,k+1], node_array[i+2,j+2,k+1], node_array[i,j+2,k+1] ))  # edges vertical
            )
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y, nel_z))
    boundary = FaceIndex[[FaceIndex(cl, 1) for cl in cell_array[:,:,1][:]];
                              [FaceIndex(cl, 2) for cl in cell_array[:,1,:][:]];
                              [FaceIndex(cl, 3) for cl in cell_array[end,:,:][:]];
                              [FaceIndex(cl, 4) for cl in cell_array[:,end,:][:]];
                              [FaceIndex(cl, 5) for cl in cell_array[1,:,:][:]];
                              [FaceIndex(cl, 6) for cl in cell_array[:,:,end][:]]]

    boundary_matrix = Ferrite.boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String,Set{FaceIndex}}()
    facesets["bottom"] = Set{FaceIndex}(boundary[(1:length(cell_array[:,:,1][:]))   .+ offset]); offset += length(cell_array[:,:,1][:])
    facesets["front"]  = Set{FaceIndex}(boundary[(1:length(cell_array[:,1,:][:]))   .+ offset]); offset += length(cell_array[:,1,:][:])
    facesets["right"]  = Set{FaceIndex}(boundary[(1:length(cell_array[end,:,:][:])) .+ offset]); offset += length(cell_array[end,:,:][:])
    facesets["back"]   = Set{FaceIndex}(boundary[(1:length(cell_array[:,end,:][:])) .+ offset]); offset += length(cell_array[:,end,:][:])
    facesets["left"]   = Set{FaceIndex}(boundary[(1:length(cell_array[1,:,:][:]))   .+ offset]); offset += length(cell_array[1,:,:][:])
    facesets["top"]    = Set{FaceIndex}(boundary[(1:length(cell_array[:,:,end][:])) .+ offset]); offset += length(cell_array[:,:,end][:])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
    
end

#Check outpput
Ferrite.cell_to_vtkcell(::Type{QuadraticHexahedron}) = VTKCellTypes.VTK_QUADRATIC_HEXAHEDRON
grid = generate_grid(QuadraticHexahedron, (3,3,3))
vtk_grid("test.vtu", grid) do vtk

end

[koehlerson]

Here is a grid generator that seems to work :)

Very nice! Thank you!

I noticed that our edge order does not follow vtk.
We order them as: bottom horizontal edges -- vertical mid edges -- top horizontal edges

So, should we follow the vtk definition?

How is the order of the dofs for your interpolation (on the edges)?

The ordering follows this graphic

[lijas]

So, should we follow the vtk definition?

yeah, that way we dont have to modify anything for the outputting, just add Ferrite.cell_to_vtkcell(::Type{QuadraticHexahedron}) = VTKCellTypes.VTK_QUADRATIC_HEXAHEDRON

[koehlerson]

ok, done. I included the test you've mentioned earlier which checks some sha sum. I have overwritten the cheksums.sha1 file now, since the QuadraticHexahedron wasn't there. I just set the overwrite variable to true and afterwards again to false. Was this the correct way? Other than that I'm going to check one last time with a computation if the export is correct

[koehlerson]

looks good :)

[lijas]

Nice :)

Why did you change FaceIndex to Tuple{Int,Int} in the gridgenerator? All other generators uses FaceIndex.

[koehlerson]

short version: derp

long version: I branched out before pulling, so I was on a relatively old master commit and therefore I wondered, why FaceIndex is not working and why you would use it in the first place 😄 . Anyhow, I think it's fine now

[koehlerson]

Let me know if it's OK to squash and merge or if further things need to be changed :)

[lijas]

Looks good to me :) Idk know if anyone else wants to take a look...

[termi-official]

Question on this: Does anyone know what the relation between the different VTK hexahedrons is?
I have never seen these 20 node hexahedrons and after a quick search I cannot find the theory behind it, other than the name CHX60 element. For completenes I also want to leave here that VTK provices the usual Lagrangian hexahedrons. See also the corresponding blog post. However, as lijas already noted, we cannot handle arbitrary order Lagrange bases yet.

A major point which I see as problematic is the naming, as we suggest to provice a "second order Lagrange hexahedron" but deliver some denegerate basis with order < 2. I would suggest to leave the "Lagrange dim 3 RefCube order 2" as the usual tensor product element of second order Lagrange basis functions. Other than the naming I like the implementation.

[koehlerson]

ye so that points to my initial doubts about the naming

In Zienkiewicz and Taylor they call three dimensional incomplete ones Serendipity:

[termi-official]

Aw yes, sorry Maxi, I derped here. You are right, the 20 node hex is a quadratic serendipity element. :) Not sure about the 24 node thingy tho.

Some infrastructure is there. However now we run into some mesh issue again. I think it makes sense to introduce some "serendipity geometry" like RefSerendipityCube? This way we can also handle serendipity elements on regular RefCubes. Alternatively there is also the possibility to fill in the missing nodes when importing a 20 node hex, reconstructing the quadratic cube geometry, i.e. RefCube.

[lijas]

This way we can also handle serendipity elements on regular RefCubes.

I think this is already possible, in the same way as we allow quadratic interpolation on linear elemets. No need for RefSerendipityCube.

[lijas]

But yeah, maybe call this interpolation serendipity instead, since there seems to be literature that calls it that :)

[koehlerson]

Some infrastructure is there. However now we run into some mesh issue again. I think it makes sense to introduce some "serendipity geometry" like RefSerendipityCube? This way we can also handle serendipity elements on regular RefCubes. Alternatively there is also the possibility to fill in the missing nodes when importing a 20 node hex, reconstructing the quadratic cube geometry, i.e. RefCube.

ye there is no need to introduce something new

But yeah, maybe call this interpolation serendipity instead, since there seems to be literature that calls it that :)

and the cell stays QuadraticHexahedron or should we name it differently?

[lijas]

and the cell stays QuadraticHexahedron or should we name it differently?

Hmm... QuadraticHexahedron seems like it should be reserved for the 27-node hexahedron?, but I dont know what we should call this 20-node hexa...

Maybe we can keep this un-aliased, and use Cell{3,20,6} for now. Unless someone comes up with a better name.

[koehlerson]

Ferrite.jl/src/interpolations.jl

Line 428 in 01cf5e1

getlowerdim(::Serendipity{3,RefCube,2}) = Serendipity{2,RefCube,2}()

Is this correct? Should have been wrong before I guess, but didn't test with face integrals where it probably matters?

Hmm... QuadraticHexahedron seems like it should be reserved for the 27-node hexahedron?, but I dont know what we should call this 20-node hexa...

Yes I think full quadratic hexas should get the name. Un-aliased seems like the best option for now

[lijas]

I dont think getlowerdim is used anywhere :P

[koehlerson]

cool, even better :D Still, I think this should be the correct dispatch in case anything needs it in the future

[termi-official]

@koehlerson This should be correct. The faces miss the central dof, which corresponds to the quadratic quadrilateral (serendipity element).
@lijas Handling serendipity interpolations on RefCubes is possible, I think we all agree here. However, I was refereing to a geometry which has just the serendipity nodes. But I am also totally fine with not having this as a feature, as it really is quite special. :)

[koehlerson]

So, I removed the alias and renamed the interpolation. I think this is again ready for review

[KristofferC]

Looks good to me. Always good to do some manual test as well, for example try solve the heat equation we have in the example and see that the result is reasonable.

[koehlerson]

So, this is how it looked with all facesets 0 Dirichlets before renaming

This is with "front" and "back" not set in the current commit:

[koehlerson]

Any other comments/suggestions or can I merge this?

[termi-official]

LGTM

Edit: Maybe we should change the title to Quadratic Serendipity Hex Interpolation for cleanness.