Change parameterization of (Abstract)Cell

[fredrikekre]

Currently we have https://github.com/Ferrite-FEM/Ferrite.jl/blob/aa261308181a5d856742385eafaa88e87481ad16/src/Grid/grid.jl#L25 and https://github.com/Ferrite-FEM/Ferrite.jl/blob/aa261308181a5d856742385eafaa88e87481ad16/src/Grid/grid.jl#L37-L39 which ... is a bit of a mess for several reasons, IMO, in particular:

1. The number of faces, M, shouldn't be a type parameter. It is not required for the storage (i.e. an NTuple or similar) and we never use it for dispatch. It was introduced in #125 because at the time we stored an M-tuple with on-boundary info. This struct field was removed in #130, but not the type parameter. In #279 the AbstractCell was introduced, and that propagated M also to the abstract type.
2. All concrete cells implemented in Ferrite are squeezed in under the same struct: https://github.com/Ferrite-FEM/Ferrite.jl/blob/aa261308181a5d856742385eafaa88e87481ad16/src/Grid/grid.jl#L85-L106. While it is a bit "cute" perhaps, it does lead to ambiguities.

Here is a suggested re-parametrization:

abstract type AbstractCell{dim,spacedim,refshape} end

struct Line{spacedim=1} <: AbstractCell{1,spacedim,RefCube}
    nodes::NTuple{2,Int}
end

struct Triangle{spacedim=2} <: AbstractCell{2,spacedim,RefTetrahedron}
    nodes::NTuple{3,Int}
end

...

Line{spacedim} can then replace Line, Line2D and Line3D, for example.

For the grid, all cells would have the same spacedim (i.e. the dim of the grid). In the grid:

struct Grid{dim=2,CT}
    cells::Vector{CT}
    # ...
end

when mixing cell types, for example Triangles and Quadrilaterals, we would then go from

CT = Union{Cell{2,3,3}, Cell{2,4,4}}

to

CT = Union{Triangle{2}, Quadrilateral{2}}

which shouldn't matter to the compiler (but let's benchmark).

[termi-official]

I agree that the parametrization is problematic. I have some experiments regarding a reparametrization here https://github.com/Ferrite-FEM/Ferrite.jl/pull/519/files#diff-aad722b9b915949ccae0c318858ae2725ae47adb50073f8b8d30092bffa1ac9f where I came to a similar conclusion. Here I have dropped the spatial dimension, but still have the number of nodes as type parameter. Dropping the number nodes is also an option. Still, I have a few questions here:

1. Why does the cell need to know about its spatial dimension? Currently we only use this as a hotfix for shell elements.
2. How should we deal with higher order geometry? Custom type QuadraticTriangle and friends or another type parameter as in the interpolation?
3. How do we couple the cell with its corresponding geometric interpolation?

[fredrikekre]

1. I guess it doesn't need to, since it will always be implied by the Grid, you mean?
2. My initial thought was to keep QuadraticTriangle just like we have now, but the order could be a parameter I suppose. Although, then the tuple storage type also depends on this so would have to be something like

   struct Triangle{order, N}
      nodes::NTuple{N,Int}
   end

   which might be a bit annoying.

3. Not sure what you mean, but in my suggestion the geometric interpolation is implicit given the reference shape parameter?

[KnutAM]

3. Coupled with the suggestion: What would be the disadvantage of coupling the cell parameterization to the geometric interpolation? I.e. struct Cell{IP<:Interpolation, N} Shouldn't that remove the mentioned type ambiguities, while still not having to implement a lot almost equal structs?

[fredrikekre]

Then you could not use e.g. linear geometric interpolation with quadratic grid though? And you couldn't reinit! CellValues with just the coordinates. Could duplicate the geometric interpolation though, but then it feels a bit strange to parametrize the cell on the interpolation (just to have it forward the same type parameters).

[KnutAM]

Then you could not use e.g. linear geometric interpolation with quadratic grid though?

Wouldn't this break current assumptions as well, e.g. reinit!(cv, getcoordinates(grid, cellnr)) fails because length(x)>getngeobasefunctions(cv).

[fredrikekre]

Okay, I thought that worked since all currently supported interpolations are hierarchical like that.

[termi-official]

1. Yes, the grid defines the spatial dimension and connectivity. From my current experiments in https://github.com/Ferrite-FEM/Ferrite.jl/pull/651 I do not find any operation where we need the spatial dimension of the element without its actual nodes.

2. I see. I also could not find a good angle here due to this dependence (and no arithmetic allowed).

3. It is only implicitly given, but the interpolation does make things unique (when we just have one Cell type). It is just another idea I had to resolve the existing ambiguity in a different way. I always thought that it is only possible to use default_interpolation for the geometry.

[KnutAM]

I think, for example, that the implicit could theoretically cause ambiguity for some case of Serendipity vs Lagrange for higher order interpolations (but I haven't checked). But if not, an alternative approach would be to have struct Cell{RefShape,order,N}

[fredrikekre]

I do not find any operation where we need the spatial dimension of the element without its actual nodes.

True, pretty much the only thing it is used for is to extract the coordinates.

I also could not find a good angle here due to this dependence (and no arithmetic allowed).

We can keep linear and quadratic as special cases, and then if there is a need, implement a HigherOrderTriangle which have nodes::Vector{Int} instead (and then you don't need the type parameter at all).

[fredrikekre]

The simplicity of

abstract type AbstractCell{refshape} end

struct Line <: AbstractCell{RefLine}
    nodes::NTuple{2,Int}
end

struct Triangle <: AbstractCell{RefTriangle}
    nodes::NTuple{3,Int}
end

is tempting... Perhaps we can insert LagrangeCell in the hierarchy to avoid ambiguities?

abstract type AbstractCell{refshape} end
abstract type LagrangeCell{refshape} <: AbstractCell{refshape} end

struct Line <: LagrangeCell{RefLine}
    nodes::NTuple{2,Int}
end

struct Triangle <: LagrangeCell{RefTriangle}
    nodes::NTuple{3,Int}
end

?

[KnutAM]

I'm missing something, but I still don't follow the advantage of splitting up into multiple cell types, Instead of either parameterizing with the reference shape (especially if we introduce a reference shape for each dimension) or the geometric interpolation?

[KnutAM]

Using IP (or RefShape + N) as type parameter would solve e.g. the following issue with the latest version:

struct QuadraticQuadrilateral <: LagrangeCell{RefQuad} or SerendipityCell{RefQuad}
    nodes::NTuple{9 or 8,Int}
end

(Which currently is Lagrange, but for e.g. importing Abaqus meshes, the standard quadratic quad is a Serendipity)

[fredrikekre]

don't follow the advantage of splitting up into multiple cell types

The tuple size dont have to be a parameter, and it felt simpler to to squeeze it all into the same type.

the following issue with the latest version:

Here I assumed that the name QuadraticQuadrilateral would be the Lagrange version, just like in the proposal of having the interpolation as a type param, where it would be defined similarly:

const QuadraticQuadrilateral = Cell{Lagrange{2,RefCube,2}, 9}

[fredrikekre]

Perhaps AbstractCell should not have any parameters at all. I don't think we have any methods dispatching on e.g. ::AbstractCell{dim} etc. It is really only useful for https://github.com/Ferrite-FEM/Ferrite.jl/blob/aa261308181a5d856742385eafaa88e87481ad16/src/Grid/grid.jl#L373-L374, but also here it could be untyped and let people put any struct in the cells vector...

Here is my new suggestion then:

No abstract types. The concrete built-in implementation Cell store the interpolation. All existing interpolations are singletons so they have no struct fields. I think storing an instance would not cost anything, and it is more general and future proof if there will ever be interpolations with data. It is also a bit easier to work with instances instead of types, IMO.

# Abstract interface
abstract type AbstractCell end

get_geo_interpolation(::AbstractCell) = error("required method")
get_reference_dim(c::AbstractCell) = get_reference_dim(get_geo_interpolation(c))
get_reference_shape(c::AbstractCell) = get_reference_shape(get_geo_interpolation(c))
get_nodes(::AbstractCell) = error("required method")

# Concrete, built-in, implementation
struct Cell{IP <: Interpolation, N} <: AbstractCell
    interpolation::IP
    nodes::NTuple{N,Int}
end

get_geo_interpolation(c::Cell) = c.interpolation
get_nodes(c::Cell) = c.nodes

# Type aliases
const Line          = Cell{Lagrange{1, RefCube,        1}, 2}
const QuadraticLine = Cell{Lagrange{2, RefTetrahedron, 2}, 3}
# etc

[termi-official]

I think we start to converge to a solution here. I like that it allows us to fix one specific geometric realization via an interpolation. A bit unfortunate that we still need the number of nodes as an extra parameter, although the information should already be contained in the interpolation implicitly.

xref https://github.com/JuliaGeometry/GeometryBasics.jl/issues/161 for more discussion regarding a general geometric interface for FE grids. I will try to bring this one forward together with @koehlerson on MakieCon next week.

[lijas]

Sorry for being a bit late in to the discussion, but I dont really see the benefit with adding the interpolation as a type parameter and/or field variable. What is being improved with that approach, which we cant do with our current method i.e. default_interpolation(::Cell)?

I agree that our cell types can be improved in some way though (like getting rid of some ambiguities). But in that case I like the proposal with different struct for each celltype, e.g.

struct Line <: AbstractCell{...}
    nodes::NTuple{2,Int}
end

[termi-official]

E.g. if your interpolation has some extra information, then default_interpolation fails.

[fredrikekre]

@termi-official do you have an example of what type of information would be needed? We can probably replace default_interpolation(::Type{<: AbstractCell}) with e.g. get_geo_interpolation(grid::Grid, cell::AbstractCell) so that you can construct your interpolation with known grid perhaps:

# Entrypoint: Hook in here for AbstractCell implementation that needs information from the grid
get_geo_interpolation(grid::Grid, c::AbstractCell) = get_geo_interpolation(c) 
# Hook in here for implementations that needs the nodes or other struct fields
get_geo_interpolation(c::AbstractCell) = get_geo_interpolation(typeof(c))
# Default stuff just like we have now
get_geo_interpolation(::Type{Line}) = Lagrange{1,RefCube,1}()

[termi-official]

A (rather stupid) example would be high order Lagrange basis functions (thing e.g. about P adaptivity). We can automatically construct these from a given node set. In this case some of this nodal information is stored in the interpolation. I mean, we might still be able to fix this with more type parameters and reconstruction the interpolation every time we call get_geo_interpolation or by manually constructing the interpolation as we do right now for the required geometric cases.

[KnutAM]

@termi-official I think you mentioned Hermite polynomials vs Lagrange as a potential source of ambiguity?

I think all discussed options solve these issues. Certainly, a specific type for each cell cannot go wrong wrt. generality, but the biggest con I see is that there will be a lot of code duplication. For example, all types would have to implement get_nodes. (And we have* to come up with unique names for each, ending up in, e.g. QuadraticSerendipityQuadrilateral.)

At least for Lagrange and Serendipity interpolations, the structure would look the same. I'm not sure how that would be with Hermite though? If we restrict Cell to have only standard vector-node-based interpolations (Lagrange and Serendipity), I think Cell{RefShape,N} could be nice (if RefShape includes dimension) because we reuse the previously defined types. Most element codes I've seen are built from "Shape + num_nodes". And for cells with different geometric interpolation, one should define a new subtype of AbstractCell (which would be nicer once there are no required parameters)

*get_node_ids would probably be a better name as it wouldn't return the actual Node

[termi-official]

Hermite is at least an obvious one, this is why I like to take it as an example (c.f. https://defelement.com/elements/examples/triangle-Hermite-3.html and e.g. https://defelement.com/elements/examples/triangle-bubble-enriched-Lagrange-1.html which have both 4 nodes. Or even on segments. Regarding the transformations, I think we require H2 conformity, which I am not sure how to enforce correctly. Technically this can be even currently resolved by "just" defining a new HermiteCell<:AbstractCell type. Not sure what the best angle is.

[KnutAM]

The Hermite-3 is a good example of challenging the question of whether it suffices to have only nodes as the cell information!

Spontaneously, I would guess that a new type of node (containing at least one position and two rotations (probably including redundant direction information for easier implementation)) would be required for the vertices. If I understand correctly, we have 4 points and 6 rotation angles. So from that point of view, it couldn't even be solved locally on the cell level?

[termi-official]

We probably want to follow existing theory for now. I think [1] is a good start.

[1] Kirby, Robert C. "A general approach to transforming finite elements." The SMAI journal of computational mathematics 4 (2018): 197-224.

[KnutAM]

"In contrast to the Lagrange element, its [Hermite. my note] node set includes function values and derivatives at the nodes, as well as an interior function value" [1]

Yes, wouldn't that be the case? (2d-gradient ~ 2 rotation angles)? (But of course a bit more challenging to implement in practice:))

[fredrikekre]

(I started writing this yesterday so you have already discussed some of these things)

For example, all types would have to implement get_nodes*

All built in cell types could have the same nodes struct field so one method could cover them all.

I think Cell{RefShape,N} could be nice (if RefShape includes dimension) because we reuse the previously defined types.

I don't like this because how do you easily connect that e.g. Cell{RefLine,3} means quadratic Lagrange, and Cell{RefLine,4} means degree 3 Hermite? I think either explicitly including the interpolation, or define multiple structs is the way to go, and then it only boils down to whether we write

struct Line ... end

or

const Line = Cell{Lagrange{1,RefCube,1}, 2}

and that doesn't seem to matter much IMO. I think I would favor the explicit Line struct since it doesn't have the tuple-size as a parameter.

Using Cell{<:Interpolation} wouldn't "just work" for any interpolation anyway since you would still have to define all of the required methods:

vertices(c::Cell{Serendipity{2,RefCube,2}, 8}) = # ...
edges(c::Cell{Serendipity{2,RefCube,2}, 8})    = # ...
faces(c::Cell{Serendipity{2,RefCube,2}, 8})    = # ...

and I am not sure writing Cell{Serendipity{2,RefCube,2}, 8} is much more convenient compared to SerendipityQuadraticQuadrilateral or something.

[KnutAM]

Using Cell{<:Interpolation} wouldn't "just work" for any interpolation

OK, I thought one could do, e.g. faces(c::Cell{<:IP}) where IP = faces(c.nodes, getrefshape(IP)) or something. But if no such generalizations are possible, perhaps separate structs are the way to go!

Edit: But I guess that would still be possible with separate structs, just by calling get_geo_interpolation(c) or get_refshape(c)