Adding NgonFace like types for volumetric elements (e.g. polyhedra)

[Kevin-Mattheus-Moerman]

This is a feature request. I am interested in using GeometryBasics in combination with volumetric meshing for finite element analysis, which features tetra-, hexa-, pentahedral elements for instance (both linear and higher order variants). I know there is some volumetric element stuff there in GeometryBasics, however unless I am mistaken they all "come with the points" e.g. P here:

const TetrahedronP{T,P<:AbstractPoint{3,T}} = Simplex{3,T,4,P}

In addition, I think it would be useful to have a polyhedron description where the points and point indices are separate. It would be great to have something equivalent to NgonFace but for polyhedra. Therefore I propose NhedraElement * and perhaps it can be achieved by adding something like this (which I believe is fully in line with how NgonFace is defined:

abstract type AbstractElement{N,T} <: StaticVector{N,T} end
abstract type AbstractNhedraElement{N,T} <: AbstractElement{N,T} end

@fixed_vector Nhedra = AbstractNhedraElement

Next we can define useful aliases, e.g. similar to TriangleFace or QuadFace. In the same vein one could propose TetrahedronElement and HexahedronElement. However I am proposing the shorter tet4 and hex8 (although I could also just implement them as aliases):

const tet4{T} = Nhedra{4,T} # A 4-noded tetrahedron
const hex8{T} = Nhedra{8,T} # An 8-noded hexahedron
const penta6{T} = Nhedra{6,T} # A 6 noded pentahedron (aka wedge element)

I am partial to, but not hung up on, these names. However I choose them in anticipation of support for higher-order finite elements (which may also feature mid-edge and mid-face points) e.g. something like:

const tet10{T} = Nhedra{10,T} # A 10-noded quadratic tetrahedron
const hex20{T} = Nhedra{20,T} # An 20-noded quadratic hexahedron

Let me know if this is of interest, I'm relatively new to Julia but I'd be happy to try add it and submit a PR.

*In an NgonFace the N should formally stand for the number of edges but it happens to be equal also to the number of points, certainly, functionally the code treats it as the latter. For the NhedraElement I propose it would be for the number of points (as the number of faces can be and often is different). The mathematical purist may feel upset though that N does not stand for the number of faces, however functionally, like with NgonFace the N really relates to the number of points we have to account for, so I'd be happy with this naming convention.

I'd also be interested to see what @fverdugo @santiagobadia @termi-official @fredrikekre @KristofferC feel about the above. They may have "baked their own cookies" over at their respective finite element related packages, but may still be able to offer insight in what they would ideally like to see on this in GeometryBasics.

[termi-official]

Things are unfortunately a bit more complicated. We had something similar in an old Ferrite version and we moved away from this idea (see https://github.com/Ferrite-FEM/Ferrite.jl/pull/679 and related discussions). Basically we had an element parametrized by the triplet (dimension, number of vertices, number of faces). Turns out this is already ambiguous. I think the simplest construction here is a quadrilateral vs bubble triangle. Both have 4 nodes. There are more subtle ones when it comes to nonlinear geometry, e.g. uniform spaces quadrilateral vs Lobatto quadrilateral. Our solution for now is basically to have one type per element. I also tried now multiple times to define something useful (see e.g. https://github.com/JuliaGeometry/GeometryBasics.jl/issues/161 and https://github.com/JuliaGeometry/GeometryBasics.jl/pull/204) but could not get something satisfactory for now.

[Kevin-Mattheus-Moerman]

@termi-official thanks for the comment. I understand what you're saying. In what I propose Nhedra does not contain information on the nature of the polyhedron other than the number of nodes. Instead using tet4 vs quad4 distinguishes the 4-noded trilinear tetrahedron from the 4-noded trilinear quadrilateral. Although for surface types one should use Ngon Rather than Nhedra so these two at least should not be confused. So I think at the very least the dimensionality (or whether something is a surface or a volume element) can be agreed upon by all I think. There are, like you say, of course polyhedra with the same number of nodes, e.g. the pentahedron and the octahedron are both (volumetric) polyhedra and both have 6 nodes. However, would it be okay to let them be distinguished by their type alias, i.e. penta6 and octa6? Is this what you mean when you say: "Our solution for now is basically to have one type per element."? It may be tricky though to have GeometryBasics have an exhaustive set of all possible polyhedra (and nonlinear ones too), and with naming conventions the community agrees upon... So perhaps we only add something like the Nhedra functionality and let users bake their own things with them?

[termi-official]

So I think at the very least the dimensionality (or whether something is a surface or a volume element) can be agreed upon by all I think.

I also think so, but I also cannot say with 100% confindence what the others think.

There are, like you say, of course polyhedra with the same number of nodes, e.g. the pentahedron and the octahedron are both (volumetric) polyhedra and both have 6 nodes. However, would it be okay to let them be distinguished by their type alias, i.e. penta6 and octa6? Is this what you mean when you say: "Our solution for now is basically to have one type per element."?

In my opinion this design is problematic, because users might think that they can write a function which specifically dispatches on e.g. penta6, but it would dispatch for both, penta6 and octa6. This is part the issue we tried to solve in PR 679 linked above. Our next thought then was: How would a generic interface for a geometry used in finite elements look like? The only common nominator we could identify here was the reference geometry (e.g. triangle, hexahedron,...). Furthermore, as pointed our above, you can also easily construct cases for e.g. nonlinear geometries where the only difference between the elements is the relative placement of the nodes. I am also not aware of a combination of characteristic parameters which can, in general, identify a polyhedron.

It may be tricky though to have GeometryBasics have an exhaustive set of all possible polyhedra (and nonlinear ones too), and with naming conventions the community agrees upon... So perhaps we only add something like the Nhedra functionality and let users bake their own things with them?

Totally agree. My point is that GeometryBasics.jl could rather provide an abstract interface which users can extend with their own geometry types rather than providing an exhaustive interface. But since you distinguish between polyhedron and n-hedron, what is the difference in definition for you?

[Kevin-Mattheus-Moerman]

@termi-official

My point is that GeometryBasics.jl could rather provide an abstract interface which users can extend with their own geometry types rather than providing an exhaustive interface.

And I was wondering if this abstraction could be something like what I propose e.g. Nhedra. And that then I would use that personally to create the tet4, hex8, penta6, etc. types.

But since you distinguish between polyhedron and n-hedron, what is the difference in definition for you?

Actually I do not distinguish between the two. An Nhedron or a polyhedron would be the same for me. Not sure if those names would be better then my Nhedra? I do think the type should relate to the number of points not the number of faces as the former is much more practical. I was only proposing Nhedra along the same lines as the current Ngon in GeometryBasics. And also because Nhedra is shorter than polyhedra. Using Nhedron or polyhedron is fine too.

So could a type be added to GeometryBasics that is for polyhedra (any volumetric element) with N-points? And to let that be the abstraction that users build their things on? I would say that the type should say nothing about the face connectivity (that all depends as we've highlighted...) and in my opinion should not ship with "per element points" so to leave the actual coordinates out. I think this type is currently missing if I'm not mistaken, hence my Nhedra proposal. Given the lack of documentation on GeometryBasics I found it hard to figure out exactly if this is missing or not. I see Simplex and NSimplex but I think these ship with the points if I'm not mistaken.

Thanks for your replies here!

[Kevin-Mattheus-Moerman]

@termi-official just updated my comment. I misread polyhedron as polygon, but fixed my comment now. I think we're on the same page anyway. :smile:

[Kevin-Mattheus-Moerman]

@termi-official if you are in agreement I can try to add this as a PR. I think the naming Nhedron is perhaps best since that is singular like Ngon.

Would it be to add something like this? :

abstract type AbstractPolyhedron{N,T} <: StaticVector{N,T} end
abstract type AbstractNhedron{N,T} <: AbstractPolyhedron{N,T} end

@fixed_vector Nhedron = AbstractNhedron