Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change parameterization of (Abstract)Cell #677

Closed
fredrikekre opened this issue Apr 11, 2023 · 26 comments · Fixed by #679
Closed

Change parameterization of (Abstract)Cell #677

fredrikekre opened this issue Apr 11, 2023 · 26 comments · Fixed by #679
Milestone

Comments

@fredrikekre
Copy link
Member

Currently we have

abstract type AbstractCell{dim,N,M} end
and
struct Cell{dim,N,M} <: AbstractCell{dim,N,M}
nodes::NTuple{N,Int}
end
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 Lots of stuff... #125 because at the time we stored an M-tuple with on-boundary info. This struct field was removed in store boundaries in a matrix #130, but not the type parameter. In Changes to MixedDofHandler #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:
    # Typealias for commonly used cells
    const implemented_celltypes = (
    (const Line = Cell{1,2,2}),
    (const Line2D = Cell{2,2,1}),
    (const Line3D = Cell{3,2,0}),
    (const QuadraticLine = Cell{1,3,2}),
    (const Triangle = Cell{2,3,3}),
    (const QuadraticTriangle = Cell{2,6,3}),
    (const Quadrilateral = Cell{2,4,4}),
    (const Quadrilateral3D = Cell{3,4,1}),
    (const QuadraticQuadrilateral = Cell{2,9,4}),
    (const Tetrahedron = Cell{3,4,4}),
    (const QuadraticTetrahedron = Cell{3,10,4}),
    (const Hexahedron = Cell{3,8,6}),
    (Cell{2,20,6}),
    (const Wedge = Cell{3,6,5})
    )
    . 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).

@fredrikekre fredrikekre added this to the v0.4.0 milestone Apr 11, 2023
@termi-official
Copy link
Member

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
Copy link
Member Author

  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
Copy link
Member

KnutAM commented Apr 11, 2023

  1. 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
Copy link
Member Author

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
Copy link
Member

KnutAM commented Apr 11, 2023

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
Copy link
Member Author

fredrikekre commented Apr 11, 2023

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

@termi-official
Copy link
Member

  1. Yes, the grid defines the spatial dimension and connectivity. From my current experiments in [BREAKING] Embedded elements #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
Copy link
Member

KnutAM commented Apr 11, 2023

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
Copy link
Member Author

fredrikekre commented Apr 11, 2023

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
Copy link
Member Author

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
Copy link
Member

KnutAM commented Apr 11, 2023

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
Copy link
Member

KnutAM commented Apr 11, 2023

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
Copy link
Member Author

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 fredrikekre changed the title Change parametrization of (Abstract)Cell Change parameterization of (Abstract)Cell Apr 12, 2023
@fredrikekre
Copy link
Member Author

fredrikekre commented Apr 12, 2023

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

Ferrite.jl/src/Grid/grid.jl

Lines 373 to 374 in aa26130

mutable struct Grid{dim,C<:AbstractCell,T<:Real} <: AbstractGrid{dim}
cells::Vector{C}
, 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
Copy link
Member

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 JuliaGeometry/GeometryBasics.jl#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
Copy link
Collaborator

lijas commented Apr 12, 2023

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
Copy link
Member

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

@fredrikekre
Copy link
Member Author

@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
Copy link
Member

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
Copy link
Member

KnutAM commented Apr 12, 2023

@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
Copy link
Member

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
Copy link
Member

KnutAM commented Apr 12, 2023

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
Copy link
Member

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
Copy link
Member

KnutAM commented Apr 12, 2023

"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
Copy link
Member Author

(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
Copy link
Member

KnutAM commented Apr 13, 2023

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)

fredrikekre added a commit that referenced this issue Apr 18, 2023
This patch reworks the `AbstractCell` interface (refer to #677 for
motivation and discussion for alternatives to this patch). Specifically,
this patch

 - introduces the reference dimension as a type parameter on `AbstractRefShape`,
 - introduces `RefInterval`, `RefTriangle`, `RefQuadrilateral`, and
   `RefHexahedron` as new subtypes of `AbstractRefShape{refdim}` (note
   that interpolations are not updated to this system yet since it can
   be done in a separate patch),
 - deprecates the old `Cell{refdim, nnodes, nfaces}` struct,
 - introduces new structs for all supported cell types, e.g.
   `struct Triangle <: AbstractCell{RefTriangle}`,
 - defines the unique vertex/edge/face identifiers for DoF-distribution
   in terms of the reference shape.

The new parametrization makes it possible to dispatch on the reference
shape instead of union of concrete implementations, i.e.
`::AbstractCell{RefTriangle}` instead of `::Union{Triangle,
QuadraticTriangle}`, and similarly for the reference dimension.

One drawback is that it is no longer possible to use "anonymous
celltypes" after this patch, i.e. using, e.g., `Cell{3, 20, 6}` and
defining methods for that instead of explicitly naming it
`SerendipityQuadraticHexahedron` or similar. Instead, you now have to
define a struct and some methods. Arguably this is a good thing --
nobody understands the meaning of `Cell{3, 20, 6}`.

For normal usage this patch is not breaking (see e.g. no required
changes for examples), unless one dispatches on the `Cell` struct
directly.

Specific things that remain to be decided:

 - `struct Triangle <: AbstractCell{2, RefTriangle}` or `struct Triangle
   <: AbstractRefCell{2, RefTriangle} <: AbstractCell`? I don't know if
   it is too restrictive to "force" `{refdim, refshape}` on all subtypes
   of `AbstractCell`. On the other hand, if one implements a
   non-refshape based element one can define a dummy reference shape.
   For a real-world example, looking at e.g.
   https://github.com/kimauth/FerriteCohesiveZones.jl/blob/main/src/cohesive_cell.jl
   I think that would simply be implemented as
   ```julia
   struct CohesiveQuadrilateral <: AbstractCell{2, RefQuadrilateral}
       nodes::NTuple{4, Int}
   end
   ```
   which I don't see a direct problem with, other than the fact that the
   cell will now have 4 edges (according to the fallback definitions).
   To solve this, one would have to implement a method for
   `faces(::CohesiveQuadrilateral)` that only return the upper and lower
   edge. I don't think this is too much to ask for custom elements like
   this. Alternatively, with #581 merged it is possible to distribute
   different number of dofs for different edges, so this can also be
   solved on the (default) interpolation level by distributing 0 dofs on
   the left and right edge.

 - `AbstractCell{refdim, refshape}` or just `AbstractCell{refshape}`?
   The dimension is implicit from the shape now, and you can dispatch on
   `AbstractCell{<:AbstractRefShape{refdim}} where refdim`. Personally I
   am in favor of including the dimension, even if redundant. It does
   not cost anything, and sometimes you want to dispatch on the
   dimension, sometimes on the refshape.

Fixes #677.
fredrikekre added a commit that referenced this issue Apr 18, 2023
This patch reworks the `AbstractCell` interface (refer to #677 for
motivation and discussion for alternatives to this patch). Specifically,
this patch

 - introduces the reference dimension as a type parameter on `AbstractRefShape`,
 - introduces `RefInterval`, `RefTriangle`, `RefQuadrilateral`, and
   `RefHexahedron` as new subtypes of `AbstractRefShape{refdim}` (note
   that interpolations are not updated to this system yet since it can
   be done in a separate patch),
 - deprecates the old `Cell{refdim, nnodes, nfaces}` struct,
 - introduces new structs for all supported cell types, e.g.
   `struct Triangle <: AbstractCell{2, RefTriangle}`,
 - defines the unique vertex/edge/face identifiers for DoF-distribution
   in terms of the reference shape.

The new parametrization makes it possible to dispatch on the reference
shape instead of union of concrete implementations, i.e.
`::AbstractCell{2, RefTriangle}` instead of `::Union{Triangle,
QuadraticTriangle}`, and similarly for the reference dimension.

One drawback is that it is no longer possible to use "anonymous
celltypes" after this patch, i.e. using, e.g., `Cell{3, 20, 6}` and
defining methods for that instead of explicitly naming it
`SerendipityQuadraticHexahedron` or similar. Instead, you now have to
define a struct and some methods. Arguably this is a good thing --
nobody understands the meaning of `Cell{3, 20, 6}`.

For normal usage this patch is not breaking (see e.g. no required
changes for examples), unless one dispatches on the `Cell` struct
directly.

Fixes #677.
fredrikekre added a commit that referenced this issue Apr 20, 2023
This patch reworks the `AbstractCell` interface (refer to #677 for
motivation and discussion for alternatives to this patch). Specifically,
this patch

 - introduces the reference dimension as a type parameter on `AbstractRefShape`,
 - introduces `RefLine`, `RefTriangle`, `RefQuadrilateral`, and
   `RefHexahedron` as new subtypes of `AbstractRefShape{refdim}` (note
   that interpolations are not updated to this system yet since it can
   be done in a separate patch),
 - deprecates the old `Cell{refdim, nnodes, nfaces}` struct,
 - introduces new individual structs for all supported cell types, e.g.
   `struct Triangle <: AbstractCell{2, RefTriangle}`,
 - defines the unique vertex/edge/face identifiers for DoF-distribution
   in terms of dispatch on the reference shape.

The new parametrization makes it possible to dispatch on the reference
shape instead of union of concrete implementations, i.e.
`::AbstractCell{2, RefTriangle}` instead of `::Union{Triangle,
QuadraticTriangle}`, and similarly for the reference dimension.

One drawback is that it is no longer possible to use "anonymous
celltypes" after this patch, i.e. using, e.g., `Cell{3, 20, 6}` and
defining methods for that instead of explicitly naming it
`SerendipityQuadraticHexahedron` or similar. Instead, you now have to
define a struct and some methods. Arguably this is a good thing --
nobody understands the meaning of `Cell{3, 20, 6}`.

For normal usage this patch is not breaking (see e.g. no required
changes for examples), unless one dispatches on the `Cell` struct
directly.

Fixes #677.
fredrikekre added a commit that referenced this issue Apr 20, 2023
This patch reworks the `AbstractCell` interface (refer to #677 for
motivation and discussion for alternatives to this patch). Specifically,
this patch

 - introduces the reference dimension as a type parameter on `AbstractRefShape`,
 - introduces `RefLine`, `RefTriangle`, `RefQuadrilateral`, and
   `RefHexahedron` as new subtypes of `AbstractRefShape{refdim}` (note
   that interpolations are not updated to this system yet since it can
   be done in a separate patch),
 - deprecates the old `Cell{refdim, nnodes, nfaces}` struct,
 - introduces new individual structs for all supported cell types, e.g.
   `struct Triangle <: AbstractCell{2, RefTriangle}`,
 - defines the unique vertex/edge/face identifiers for DoF-distribution
   in terms of dispatch on the reference shape.

The new parametrization makes it possible to dispatch on the reference
shape instead of union of concrete implementations, i.e.
`::AbstractCell{2, RefTriangle}` instead of `::Union{Triangle,
QuadraticTriangle}`, and similarly for the reference dimension.

One drawback is that it is no longer possible to use "anonymous
celltypes" after this patch, i.e. using, e.g., `Cell{3, 20, 6}` and
defining methods for that instead of explicitly naming it
`SerendipityQuadraticHexahedron` or similar. Instead, you now have to
define a struct and some methods. Arguably this is a good thing --
nobody understands the meaning of `Cell{3, 20, 6}`.

For normal usage this patch is not breaking (see e.g. no required
changes for examples), unless one dispatches on the `Cell` struct
directly.

Fixes #677.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

4 participants