Unite Face subtypes and Polygon type

Hi, I was experimenting with the library and saw that

julia> Quadrangle <: Polygon
false
julia> Triangle <: Polygon
false

which I found a bit odd at first. Then I realized that Polygon is a concrete type. I might be wrong, but wouldn't any 2-face be a polygon? I also understand that the current Polygon struct serves to construct any kind of polygon, simple or not, including triangles and quadrangles (geometrically, not Triangle nor Quadrangle instances).

Maybe would it be better to define a Polygon abstract type, and have the current Polygon concrete type be something like a PolygonArea or GenericPolygon?

I guess it could also be possible to use a new Polygon abstract type in a trait-based approach instead, which would apply to Polytope objects.

I was also thinking about this a lot, and ended up naming it a Polygon to match what most people would expect in the GIS world. However, I agree that this is confusing, and we should probably rename it to Ngon or some more specialized name.

We have Quadrangle <: Polytope and that works as expected. If we rename Polygon -> Ngon I think it will be more difficult to expect people trying Quadrangle < Ngon.

Would you like to submit a PR with this renaming?

I have been thinking some more about this, and I believe the confusion may come from the definition of the abstract types Polytope and Face. In the following I assume that the current Polytope type refers to a geometric polytope, and not to an abstract polytope.

You can classify a k-polytope by its "dimension" or "number" k. You can have a polygon or a polyhedron for example for k=2 and k=3. But a polygon can be defined in 2D, 3D or higher, same as a polyhedron can be defined in 4D. In addition to geometry parameters (embedding space dimension Dim and element type T), the polytope should probably contain this k parameter to better describe it and e.g. differentiate between polygons and polyhedra.

Now, I believe the abstract polytope theory to be useful in bringing another level of abstraction. According to Wikipedia, an abstract polytope does not depend on the dimension of the Euclidian space. Only the geometric polytope does:

An ordinary geometric polytope is said to be a realization in some real N-dimensional space, typically Euclidean, of the corresponding abstract polytope.

Therefore, an abstract polytope type should only possess a parameter which refers to its dimension, but not to the space of its geometric realization. From there, a geometric polytope type could be defined, parametrized by its embedding space and the abstract polytope it implements.

I am also quite unsure about the use of the Face type. A face is defined as an element that constitutes a (abstract) polytope. So, a polyhedron can be declined into (2-)faces, which in turn can be declined into edges (1-face), and finally into vertices (0-face). A 2-face of a polyhedron is a polygon, which itself is a polytope. However, that is about it: a 2-face can be any kind of polytope, and therefore should not require any special treatment, unless one wants to link it with its parent polytope. This works for any n-face.

The current Polygon is merely a special representation which allows to build any kind of polygons. And any polygon could in theory be used for finite element analysis, although having uniform and simple polygon types is usually required for most algorithms.

I am proposing to implement something like this:

abstract type AbstractPolytope{N} end

const AbstractPolygon = AbstractPolytope{2}
const AbstractPolyhedron = AbstractPolytope{3}

# geometric polytope
abstract type Polytope{N, Dim, T} <: Geometry{Dim, T} end

const Polygon = Polytope{2}
const Polyhedron = Polytope{3}

# not used for now, but may be useful for dispatching on abstract polytopes
# (since a Polytope has already a Geometry as supertype, and a struct can only have one supertype)
implements(::Type{Polytope{N}}) where {N} = AbstractPolytope{N}

# V should contain 3 elements
struct Triangle{Dim, T, V<:AbstractVector{Point{Dim,T}}} <: Polygon{Dim, T}
    vertices::V
end

# V should contain 4 elements
struct Quadrangle{Dim, T, V<:AbstractVector{Point{Dim,T}}} <: Polygon{Dim, T}
    vertices::V
end

struct BoundedPolygon{Dim,T,C<:Chain{Dim,T}} <: Polygon{Dim,T}
  outer::C
  inners::Vector{C}

  function BoundedPolygon{Dim,T,C}(outer, inners) where {Dim,T,C}
    @assert isclosed(outer) "invalid outer ring"
    @assert all(isclosed.(inners)) "invalid inner rings"
    new(outer, inners)
  end
end

I think this would be clearer and closer to the mathematical definition. That would require some refactoring other than a renaming because with this reasoning Face is not relevant anymore and should be deleted.

I am willing to submit a PR to implement this new model, if you agree that this would be an improvement.

You can classify a k-polytope by its "dimension" or "number" k. You can have a polygon or a polyhedron for example for k=2 and k=3. But a polygon can be defined in 2D, 3D or higher, same as a polyhedron can be defined in 4D. In addition to geometry parameters (embedding space dimension Dim and element type T), the polytope should probably contain this k parameter to better describe it and e.g. differentiate between polygons and polyhedra.

It is already there. The embedding dimension is embeddim and the k in your explanation I called paramdim for parametric dimension.

I am also quite unsure about the use of the Face type. A face is defined as an element that constitutes a (abstract) polytope. So, a polyhedron can be declined into (2-)faces, which in turn can be declined into edges (1-face), and finally into vertices (0-face). A 2-face of a polyhedron is a polygon, which itself is a polytope. However, that is about it: a 2-face can be any kind of polytope, and therefore should not require any special treatment, unless one wants to link it with its parent polytope. This works for any n-face.

This is also already there. I am having the impression that you are just describing the current version of the source code. If you check each Face type, there is a facets function that returns the polytopes of lower dimension.

Regarding your proposal, I spent some time thinking about how to organize the types so that they are useful for multiple dispatch. There is not much use in splitting Polyhedron from Polygon for example, they are all Polytopes. If someone is interested in querying only polygons, the paramdim is known at compile time just for that.

The only thing that makes sense to me, is maybe rename the Polygon to Ngon or some more specialized name to avoid the confusion.

Regarding the need of the Face type, it is there to make it explicit what is acceptable as a finite element in a mesh. Not all polytopes are valid finite elements in meshing algorithms. So it is nice to be dispatch on Face only, instead of all kinds of polygons.

It is already there. The embedding dimension is embeddim and the k in your explanation I called paramdim for parametric dimension.

What benefit is there in not using a type parameter to describe this parametric dimension? The same could be done as in embeddim: the type is parametrized by Dim, and the function embeddim is a convenience function for querying Dim. You can define a paramdim in the exact same way. Since a polygon is just a polytope with 2 as its parametric dimension, having this number as a type parameter would allow the definition of a const alias. I am not saying that it is important to separate polygon or polyhedron elements, as you say it describes the same thing in the end (up to a parameter). It would just be clearer.

By the way, paramdim throw an error for any Polygon instance, while it is supposed to be a polytope of parametric dimension 2. Having this struct be a subtype of the right alias would fix the problem automatically.

The only thing that makes sense to me, is maybe rename the Polygon to Ngon or some more specialized name to avoid the confusion.

Ngon and Polygon are the same thing so I see little point in doing so. It is not about avoiding people to notice that it is confusing. It's about fixing the confusion. It would need another name.

This is also already there. I am having the impression that you are just describing the current version of the source code.

The code aims at implementing this mathematical background, and it is very well organized. I looked at most of the source code already. I am not pointing out any missing functionality, it's just that to me the application of the mathematical concepts is semantically wrong in the few places I mentioned. I believe clarity is important, which is why I'd like to fix these blurry corners.

If you check each Face type, there is a facets function that returns the polytopes of lower dimension.

Well, that's not really what I meant. Saying that some element is a face would only make sense if you consider the polytope it is a face of. Also, any polytope has facets, not just the implementations of what you identified as faces. In its current state, it brings confusion.

Regarding the need of the Face type, it is there to make it explicit what is acceptable as a finite element in a mesh. Not all polytopes are valid finite elements in meshing algorithms. So it is nice to be dispatch on Face only, instead of all kinds of polygons.

How do you define what is acceptable as a finite element? This looks somewhat arbitrary. If you want to dispatch on only some polytopes, you could add a trait for that, right? Unless there would be a noticeable loss of performance - I haven't used traits much. It just feels wrong to me to call these "acceptable" elements "faces". They are faces, yes, of the whole model when applying finite element methods, but sometimes you want just a Quadrangle without considering FEM. Well, that's the case for me in 2D rendering. I know I could use the current Polygon struct instead, but it just want it to be made of 4 vertices without inner chains. So, a quadrangle. Again, I know there is no loss of functionality, because you can easily convert and such, but it feels wrong to specifically label these as "elements", if you don't intend this package to be opinionated.

Let me say that I like the feedback. It is good polish these types as much as possible as we progress. Below are some quick comments before I call the day here. I can continue the discussion tomorrow.

What benefit is there in not using a type parameter to describe this parametric dimension? The same could be done as in embeddim: the type is parametrized by Dim, and the function embeddim is a convenience function for querying Dim. You can define a paramdim in the exact same way.

That is a good point, and I remember iterating over this design some weeks ago. The problem I faced specifically with the parametric dimension is that it is a property of the type of the polytope, not the instance. So for example, even if we introduce a type parameter Polytope{Dim,T,P} corresponding to the convenience function paramdim, we will still need to specialize paramdim for every "face" type otherwise paramdim(Triangle) won't work. Why do I need paramdim(Triangle) and not paramdim(Triangle{Dim,T}) you may ask? It is because the meshes are lazy views on a point set. We never construct faces nor elements until the last minute and they are just "pointers" to an array of vertices. The connectivity lists of the mesh don't contain the necessary information to instantiate the type. We also need a list of actual points in an embedding space to materialize the polytopes and the points contain the missing Dim and T information.

Additionally, the more type parameters the more time we spend in compilation of the package. We add type parameters when we feel that they will be used a lot in dispatch. However, so far I couldn't find good use cases for an additional type parameter. We may change mind in the future, specially if you have some application in mind that requires dispatching specifically on 2-polytopes embedded in N > 2 spaces.

By the way, paramdim throw an error for any Polygon instance, while it is supposed to be a polytope of parametric dimension 2. Having this struct be a subtype of the right alias would fix the problem automatically.

The fix is simple as you know, we just need to specialize paramdim. The Polygon type is the one that I am working currently, so it is still work in progress.

Ngon and Polygon are the same thing so I see little point in doing so. It is not about avoiding people to notice that it is confusing. It's about fixing the confusion. It would need another name.

I partially agree. Ngon can also be interpreted as an arbitrary "N" gon without a particular name or structure, like the USA map. It would be nice to find a better name, and I am been thinking about it really hard, and couldn't find one that is good and short yet.

The code aims at implementing this mathematical background, and it is very well organized. I looked at most of the source code already. I am not pointing out any missing functionality, it's just that to me the application of the mathematical concepts is semantically wrong in the few places I mentioned. I believe clarity is important, which is why I'd like to fix these blurry corners.

That is very welcome, and I really appreciate you taking the time to review the code. I am 100% with you that we need to polish things as much as possible. It is a huge effort to get the names and type hierarchy right from the beginning. We should continue with discussions like this as we progress.

Well, that's not really what I meant. Saying that some element is a face would only make sense if you consider the polytope it is a face of. Also, any polytope has facets, not just the implementations of what you identified as faces. In its current state, it brings confusion.

Agree, and this is because I was thinking about meshes when I introduced this type. Perhaps there is a better name to describe this concept, which is the concept of a geometry that can be used to discretize a volume of given dimension. Maybe Element is a better name to avoid confusion.

How do you define what is acceptable as a finite element? This looks somewhat arbitrary. If you want to dispatch on only some polytopes, you could add a trait for that, right?

Finite elements have been widely documented in FEM literature. They are specific polytopes with a pre-defined order of vertices that facilitate the definition of interpolation functions of high-order. I am adopting the conventions of the GMSH project: https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering

They are faces, yes, of the whole model when applying finite element methods, but sometimes you want just a Quadrangle without considering FEM.

How is the current Quadrangle restricting your use case? Apart from inheriting from Face, it is also considered a Polytope in the current type hierarchy. We can brainstorm this more tomorrow or during the week.

That is a good point, and I remember iterating over this design some weeks ago. The problem I faced specifically with the parametric dimension is that it is a property of the type of the polytope, not the instance. So for example, even if we introduce a type parameter Polytope{Dim,T,P} corresponding to the convenience function paramdim, we will still need to specialize paramdim for every "face" type otherwise paramdim(Triangle) won't work.

If we put the parametric dimension in the Polytope type itself, you could construct other types with a fixed parametric dimension by "instantiating" a type from it. I believe this could be enough. The intention is not to use a general polytope with an unknown parametric dimension for future supertypes, but a Polytope{2, Dim, T} where {Dim, T} for example. From the snippet above

const Polygon = Polytope{2}

struct Triangle{Dim, T, V<:AbstractVector{Point{Dim,T}}} <: Polygon{Dim, T}
# skipped

if we define the paramdim function as

paramdim(::Type{<: Polytope{N}}) where {N} = N

then the parametric dimension of the Triangle is automatically specified when inheriting from Polygon. There is no need to specify it for each "face" type manually. It is implicitly defined from the polytope you inherit from. The Triangle is the same as before in terms of parameters. Except that we now have paramdim(Triangle) = 2 with no effort.

I partially agree. Ngon can also be interpreted as an arbitrary "N" gon without a particular name or structure, like the USA map. It would be nice to find a better name, and I am been thinking about it really hard, and couldn't find one that is good and short yet.

A suitable name is indeed hard to find... From a computer graphics perspective, a Ngon usually refers to a planar polygon with N sides, without holes. I have been thinking of PolyArea since the struct implements a polygon in terms of outer and inner boundaries to form a polygonal area. I just realized though that with this implementation there might (?) be a problem for some methods, e.g. plotting recipes, when the embedding space is in 3D, if the chains do not lie in the same plane... but that is probably out of scope for this discussion.

Agree, and this is because I was thinking about meshes when I introduced this type. Perhaps there is a better name to describe this concept, which is the concept of a geometry that can be used to discretize a volume of given dimension. Maybe Element is a better name to avoid confusion.

Element seems on the right track to me. Although, it is quite broad unless you specifically consider it from within the FE field, where the meaning is clearly unambiguous. Something like FElement maybe? I feel that one should not have to know about finite elements at all to go through the mathematical model. What I liked from a trait-based approach is that you separate the mathematical structures from their field-related use, but traits may not be the way to go. I would see for example a FEM package using Meshes, more than Meshes defining finite elements. Of course, if you want to keep everything in one place for your use cases, it could at least be somewhat isolated from the mathematical model, for clarity. Having a separate FEM module, or just put all FEM-related stuff inside a separate folder. With the current Face (or Element or whatsoever) type, a Triangle is defined as a finite element. But maybe it should be the finite element logic that externally exploits properties of simple polytopes, and inject FE meaning into it only at that stage.

Finite elements have been widely documented in FEM literature. They are specific polytopes with a pre-defined order of vertices that facilitate the definition of interpolation functions of high-order. I am adopting the conventions of the GMSH project: https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering

This could reside inside a specific finite element logic, since most of these nodes are not necessarily useful for defining what a Triangle or Quadrangle is for example (maybe Quadrilateral may be clearer when not talking about FEM quads, as we're here?). I see two different consequences of following this convention for these types. First, it allows the definition of multiple elements with the same shape, but more than the minimum required number of vertices. I mean, you need only 3 and 4 vertices to represent a triangle and quadrangle respectively, not 9 - until you apply FEM. Second, it automatically assigns k-faces, k > 0 from a set of vertices, from an arbitrary order which forms the convention. There is a loss of generality with this, because you can't use these polytopes for algorithms that respect a different convention, unless you define new types for semantically the same thing.

Let me present that in a different way. If we come back to their mathematical definition, polytopes are partially ordered sets. For a polyhedron for example, you need vertices, but also edges and "faces". You can get edges and then vertices from a set of "faces", but not vice-versa. The GMSH conventions act as a mapping (connectivity function) between the vertices and the faces (in the mathematical sense) of identified objects, which allow one to only use vertices to describe these objects using a "hard-coded" logic. This reasoning can be of course applied to any polytope, including this "Ngon" representation, where the logic is to follow clockwise and anticlockwise orderings of vertices as seen from a point inside a given chain. It would be really great to make such mappings implementable by the user, and possibly ship the GMSH convention by default. I hope this can be done with no runtime cost and minimal precompilation cost. It would be awesome to have some benchmark suites as a guide/safeguard, by the way.

Now, with what I just said, stripping connectivity would render the current "Element" types deprived of all interest. I would instead recommend the following:

add an ordering type parameter to the common polytopes used for FEM (triangles, quadrangles...)
define a ConnectivityMap (or other name) abstract type, with one subtype GMSH
specialize the facets function with methods that take a ::Type{<: Triangle{GMSH}} argument (replace Triangle for other specializations with the relevant type)

I hope all this makes sense... the more I dive into it the more it does, for me anyway.

I forgot to answer your question:

How is the current Quadrangle restricting your use case? Apart from inheriting from Face, it is also considered a Polytope in the current type hierarchy. We can brainstorm this more tomorrow or during the week.

For now it does what I want in that it stores vertices and connects them. But if for some reason it gets a different treatment because it is considered a finite element, then I'd have to convert it to something else. An example might be triangulation. I don't know if triangulating FEM quads would make sense, while it probably would for geometric quads if you consider any purpose. For quads it's easy, but it could be gain relevance for more complex (yet common) shapes. I know, it's not hard to just add a triangulate method for faces, but something looks off to me because of this FEM connotation. Although it's not really seen from the exterior of the package for now, it is when you look at the code.

Really good suggestions!

If we put the parametric dimension in the Polytope type itself, you could construct other types with a fixed parametric dimension by "instantiating" a type from it. I believe this could be enough. The intention is not to use a general polytope with an unknown parametric dimension for future supertypes, but a Polytope{2, Dim, T} where {Dim, T} for example. From the snippet above
const Polygon = Polytope{2}

struct Triangle{Dim, T, V<:AbstractVector{Point{Dim,T}}} <: Polygon{Dim, T}
# skipped
if we define the paramdim function as
paramdim(::Type{<: Polytope{N}}) where {N} = N
then the parametric dimension of the Triangle is automatically specified when inheriting from Polygon. There is no need to specify it for each "face" type manually. It is implicitly defined from the polytope you inherit from. The Triangle is the same as before in terms of parameters. Except that we now have paramdim(Triangle) = 2 with no effort.

Nice! Would you mind submitting a PR with the proposed change? We could then feel how it works out with what is already implemented in the test suite.

A suitable name is indeed hard to find... From a computer graphics perspective, a Ngon usually refers to a planar polygon with N sides, without holes. I have been thinking of PolyArea since the struct implements a polygon in terms of outer and inner boundaries to form a polygonal area.

That one deserves a PR too! What do you think of PolySurface instead to encompass the case in which the polygon is embedded in 3D?

I just realized though that with this implementation there might (?) be a problem for some methods, e.g. plotting recipes, when the embedding space is in 3D, if the chains do not lie in the same plane... but that is probably out of scope for this discussion.

I think this shouldn't be an issue if the plot recipe is aware of this possibility. I am working on a pure Julia triangulation, and we will soon be able to test it.

I feel that one should not have to know about finite elements at all to go through the mathematical model.

I like this idea too. The fact that some of these polytopes can become mesh elements is a detail that we need to separate well in the codebase. We need a type to dispatch on, but there is room for improvement regarding the conventions as you mentioned. Perhaps we can address this part of the issue after the new Polytope type is in place.

Let me present that in a different way. If we come back to their mathematical definition, polytopes are partially ordered sets. For a polyhedron for example, you need vertices, but also edges and "faces". You can get edges and then vertices from a set of "faces", but not vice-versa. The GMSH conventions act as a mapping (connectivity function) between the vertices and the faces (in the mathematical sense) of identified objects, which allow one to only use vertices to describe these objects using a "hard-coded" logic. This reasoning can be of course applied to any polytope, including this "Ngon" representation, where the logic is to follow clockwise and anticlockwise orderings of vertices as seen from a point inside a given chain. It would be really great to make such mappings implementable by the user, and possibly ship the GMSH convention by default. I hope this can be done with no runtime cost and minimal precompilation cost. It would be awesome to have some benchmark suites as a guide/safeguard, by the way.

Amazing. I really love this idea of splitting the FE logic from the mathematical concepts, and providing convention types with actual orderings. I suggest we first address the simpler issues above, and then open a separate issue to discuss just that.

Please let me know if you would like to work on the first two PRs I mentioned above. I can give it a try over the week if you are busy with other things. Then, we can open a new issue to start polishing the ordering conventions and other concepts related to FE.

Good! I'll make the two PRs.

Some renaming didn't get through. My mistake - the whole word option was toggled for search/replace and filenames are not affected by the search so some got through as well. I'll add a PR with these fixes.

I have made a series of experiments that implement my vision of the discussion above regarding the separation of FEM elements. I'll open a new issue to this end and close this one, since the original confusion is resolved via #7 and #8. I have converged towards a working prototype and I am quite satisfied with the result. I will submit it in several PR chunks, because there are changes that are not specifically related to FEM stuff but that do simplify, that way we can discuss each of them.

JuliaGeometry / Meshes.jl

Unite Face subtypes and Polygon type #6