Suggestion to change faces(::Quad) to edges(::Quad) #789
Conversation
|
Thanks for the suggestion and working on this. A few comments on this. Deal.ii uses faces for lines in 2D. MFEM for example does not have a consistent definition for this. I also asked a bit around and the response whether people would call the boundary of a 2D element an edge or face is actually pretty split. Saying an edge is a topological 1D entity further raises the question whether the cell for a line element is also an edge, and e.g. how do we handle this edge case in the code. Regarding the removal of the internal hacks I thought about having an dimension-independent layer in the Ferrite core, which is independent of the actual definition of face/edge/... . |
|
Would be really nice if this could simplify the internals! IMO any such change must be done without removing the possibility of straight-forward writing dimension-independent code. Would this also turn |
|
@termi-official
I think this is true, a line element would also be an edge, but what would be the problem with this? What edge-cases (pun intended? :P ) needs to be handled. @KnutAM Regarding FaceValues, I think it would be nice to introduce FacetValues which returns FaceValues for 3d-surfaces and EdgeValues for Lines. |
|
#780 I have there a lot of |
| @@ -70,9 +70,14 @@ mutable struct ExclusiveTopology <: AbstractTopology | |||
| edge_edge_neighbor::Matrix{EntityNeighborhood{EdgeIndex}} | |||
| # lazy constructed face topology | |||
| face_skeleton::Union{Vector{FaceIndex}, Nothing} | |||
| edge_skeleton::Union{Vector{EdgeIndex}, Nothing} | |||
There was a problem hiding this comment.
This seems only valid for 2D and should rather be facet_skeleton ?
There was a problem hiding this comment.
Hmm, I dont think so. My idea is that you generate the edge_skeleton, face_skeleton etc the same for both a 2d and 3d grid (because topologicaly they are definied in the same way now), and then if you want to loop over the dim-1 entity, you would dispatch on the dim-parameter:
function facetskeleton(top::ExclusiveTopology, grid::Grid{2})
return top.edge_skeleton
end
function facetskeleton(top::ExclusiveTopology, grid::Grid{3})
return top.face_skeleton
end
There was a problem hiding this comment.
I think I'm wondering in general about the dispatch. If I'm not mistaken, then you define a facet as a dim-1 entity and you can store directly the facetskeleton without having the extra field in the struct
function facetskeleton(top::ExclusiveTopology, grid)
return top.facetskeleton
end
Or is there some advantage of having both fields?
There was a problem hiding this comment.
But what if you have a hex mesh? Then you might want to have both an edge_skeleton and a face_skeleton. Then you will have a extra field called facet_skeleton which is identical to face_skeleton?
I realise now though that the dispatch approach might not work for shells in 3d... the dispatch for facetskeleton(..., grid{3}) would return the face_skeleton, even though you probably would want egde_skeleton...
There was a problem hiding this comment.
Ah okay! I think as of now, there is no construction of edge skeleton in 3D cases planned. Not sure if they are useful in some cases, but if we want to include them, then this makes sense
There was a problem hiding this comment.
Following my comments below, I think this should be facetskeleton and ridgeskeleton, that would allow having the edge skeleton in 3d, but for the current cases only facetskeleton would be used.
I would have thought the other way around, that the actual type is My main fear with those definitions (but they are what they are, so can't really change the name), is that I find |
|
@KnutAM Those are my concerns too. Regarding FacetValues... I think you are right, we might need to introduce a new type FacetValues, because EdgeValues typically does not store any normals, which is still something you might want for a 2d problem even though the boundary is technically an edge. |
|
I am fine with renaming the stuff, but we should be consistent to avoid further confusion. Hence, for consistency of the API we should also discuss the (re-)naming of:
I think the confusion can be somewhat complied with by providing additional docs and an FAQ. For me the bothering part when defining faces and edges by dimension (instead of by codimension) in the implementation is the implication on the data structures for the user-facing part. What exactly should |
revert some changes revert som changes
|
@termi-official and @KnutAM, can you make a review of this PR? I know you have thought about this problem, so I would like your input here. I think that I have convinced Fredrik and Kim that this is a good solution (although the breaking changes are of course not nice), but I imagine you have a different approach in mind . Note, there are still stuff to do in this PR, like deciding if we should rename FaceValues, docs, fixing tests.... But I think the point of the PR is clear. To reiterate the main idea of this PR. We keep the names vertices, edges and faces (e.g in
As an example, the sides of a Quadrilateral would no longer be called This change simplifies the code-base quite a lot at some places, especially in the dof-distribution, and solves most(all?) of the consistency issues we have The major drawback is now that users can no longer call e.g. Similar dispatches are used through out the code, e.g for |
|
Thanks for the work on solving the issue here!
Sure I will take a look. I do not think that I have a different solution per-se, but I share similar concerns as Kim about the entry barrier to Ferrite with these changes (and I mean with all possible solutions, not specific to this PR).
I already discussed this with Fredrik and want to bring this also up here: Now FaceIndex and EdgeIndex is a bit ambigous. Maybe we should introduce FacetIndex and RidgeIndex, because this is actually what these indices refer to on master.
I would rather reserve vertex for corners and nodes for points in general. I think this is a bit more consistent with our intent. We really need to clarify this in the docs, because I think this is already on master super confusing for users.
I think we all can agree on this definition.
I would drop the normal part here, because defining what the normal (field) is can be troublesome. E.g. what is the normal of a face in 4D or 2D?
For me personally this is also fine. But we should reflect this somewhere in the docs very carefully and visible for beginners. No idea how tho. |
|
Before I go for a full review, can you also update the remaining files in FEValues and interators.jl? How should we proceed with the current boundary indices Since we are at it
|
There was a problem hiding this comment.
Thanks for starting this huge effort @lijas!
Some initial comments/thoughts
My comments are based on my understanding that co-dimension is relative the spatial dimension (as opposed to the reference (element) dimension), so with my understanding, for RefLine embedded in 2d, the edge is the facet and the point is the ridge, but in a non-embedded case (1d), the edge is the cell and the point is the facet. Does this match your view?
Edit: Updated comments after I understood that we use the codim relative the topological/reference shape dimension.
I would have expected that for the grid we work with co-dimensional entities (basically as we were before, except we had the wrong names), i.e. we should just rename
FaceIndex -> FacetIndex
EdgeIndex -> RidgeIndex
VertexIndex -> PeakIndexinstead of changing depending on the sdim in the grid.
As @termi-official mentions, I think we should have FaceValues -> FacetValues and FaceQuadratureRule to FacetQuadratureRule. AFAIU, from the user perspective, one would for most cases only need to care about these names (cell, facet, ridge, peak).
Changing to using the dimensional entities (volume, face, edge, point) for interpolations, reference shapes, and dof distribution looks like it works great!
src/Grid/grid.jl
Outdated
| @inline facets(c::AbstractCell{<:AbstractRefShape{1}}) = vertices(c) | ||
| @inline facets(c::AbstractCell{<:AbstractRefShape{2}}) = edges(c) | ||
| @inline facets(c::AbstractCell{<:AbstractRefShape{3}}) = faces(c) |
There was a problem hiding this comment.
I think this is wrong, since facet refers to sdim - 1 and not rdim - 1?
There was a problem hiding this comment.
I unresolved this, because I think we should add this information to the documentation to have a single source of truth for the definitions. @lijas
| struct EdgeIterator{FC<:FaceCache} | ||
| fc::FC | ||
| set::Set{EdgeIndex} | ||
| end | ||
|
|
||
| function EdgeIterator(gridordh::Union{Grid,AbstractDofHandler}, | ||
| set::Set{EdgeIndex}, flags::UpdateFlags=UpdateFlags()) | ||
| if gridordh isa DofHandler | ||
| # Keep here to maintain same settings as for CellIterator | ||
| _check_same_celltype(get_grid(gridordh), set) | ||
| end | ||
| return EdgeIterator(FaceCache(gridordh, flags), set) | ||
| end | ||
|
|
||
| @inline _getcache(fi::EdgeIterator) = fi.fc | ||
| @inline _getset(fi::EdgeIterator) = fi.set | ||
|
|
||
| struct VertexIterator{FC<:FaceCache} | ||
| fc::FC | ||
| set::Set{VertexIndex} | ||
| end | ||
|
|
||
| function VertexIterator(gridordh::Union{Grid,AbstractDofHandler}, | ||
| set::Set{VertexIndex}, flags::UpdateFlags=UpdateFlags()) | ||
| if gridordh isa DofHandler | ||
| # Keep here to maintain same settings as for CellIterator | ||
| _check_same_celltype(get_grid(gridordh), set) | ||
| end | ||
| return VertexIterator(FaceCache(gridordh, flags), set) | ||
| end | ||
|
|
||
| @inline _getcache(fi::VertexIterator) = fi.fc | ||
| @inline _getset(fi::VertexIterator) = fi.set | ||
|
|
||
| FaceIterator(gridordh::Union{Grid,AbstractDofHandler}, set::Set{EdgeIndex}, flags::UpdateFlags=UpdateFlags()) = | ||
| EdgeIterator(gridordh, set, flags) | ||
|
|
||
| FaceIterator(gridordh::Union{Grid,AbstractDofHandler}, set::Set{VertexIndex}, flags::UpdateFlags=UpdateFlags()) = | ||
| VertexIterator(gridordh, set, flags) | ||
|
|
There was a problem hiding this comment.
Cannot a FacetIterator replace the current FaceIterator?
Elias and I discussed this based on https://defelement.com/ciarlet.html#Cell+subentities where they explicitly define the names based on codimension of the reference element (topological dimension) and not the spatial dimension (which they refer to as geometric dimension). So if we have 3 spatial dimension with hex and (embedded) quads we have 8 facets of the hex which are faces and we have 4 facets of the quad which are edges. |
Thanks for correcting me (I have to rethink a bit:)) I guess the implication for the example in #899 then would be that the "face" of the beam would not be accessible for adding boundary conditions to (using facets), but the endpoints would be available. And when thinking of this now, that makes sense because in the non-embedded case a load on the "face" of the beam would actually be a body load! I'll go through again and update my comments accordingly |
src/Grid/grid.jl
Outdated
| getfacetsets(grid::AbstractGrid{1}) = grid.vertexsets | ||
| getfacetsets(grid::AbstractGrid{2}) = grid.edgesets | ||
| getfacetsets(grid::AbstractGrid{3}) = grid.facesets |
There was a problem hiding this comment.
Since facet refers to rdim-1 this isn't always true then, right?
So it would be better to actually just store the peaksets, ridgesets, and facetsets directly?
Shouldn't this be renamed to "point" following the table, and "vertex" refers to the zero-dimensional reference shape? |
If we do |
I don't think AFAIU |
|
I just meant that we should keep the distinction between corners which are vertices currently and all 0 dimension entities of an element, which should be currently
https://defelement.com/elements/examples/quadrilateral-lagrange-gll-1.html Here for instance, they refer still to the corners as vertices, which is fine for me. I was just unsure if the table also gathers the inner 0 dimension entities as vertices. |
|
Interesting, so either I have further misunderstandings, or they are inconsistent 🤔 (But language-wise I'd switch vertex and point in the table) |
Yea, and thats why it make sense to use |
|
I mean that if it was relative the spatial dimension, then when integrating for example an embedded quad in 3d you would use a FacetValue and when integrating an embedded beam in 3d you would use RidgeValues. |
|
Aha, yes - I agree that codim relative rdim/topological dimension makes more sense! Just to put up the a summary of fenics' definitions based on their manual, I think the only difference from defelement is the Co-dimensionally defined quantities
Topologically dimensionally defined quantities
Other relevant names
|
In this definition of a vertex you set vertex == nodes or am I missing something? This is only valid for linear stuff. If we want to distinguish corner points from interior points |
That's not how I understand it (but these definitions are still confusing me a bit), when reading up on this I found this resource on the Ciarlet definition insightful, where the nodal values are defined as the support variables for the interpolation (https://finite-element.github.io/L2_fespaces.html). So as I understand, in general a node is not understood as a geometrical point at all (for example for the vector interpolations such as Nedelec where the nodal values are the integrated along an edge). However, in the case of e.g. Lagrange interpolations, all nodes coincide with the interpolated value at specific geometrical points. (This goes against my previous understanding of a node though, so I do see the risk confusion). When it comes to the geometrical nodes, since our interpolations are scalar we have the correspondance between a geometrical node and point (I guess in IGA this doesn't hold @lijas?). However, vertices are specific to the topology and refers specifically to the endpoints of intervals (i.e. the codim 1 of an interval/line). (But I agree that "MeshEntity of topological dimension 0" doesn't fully cover that) |
|
I'm just saying if we define something that includes all 0 dimension entities, there should be something for the corners, since they are very relevant for topological (geometry) constructions which was previously in the code base |
This is my (current) understanding too. |
|
I'm confused by the book. For triangles it seems that corners == vertices, but for quads all points == vertices (they are defined with a single index)
vs. |
|
Yes, I don't get their numbering scheme either, e.g. why 3 indices are used for 3rd order cases (I would have assumed
But I think we can agree that
? |
|
Maybe it is easier for everyone if we just rename the thing to "corner"? I think this is the most direct and unambiguous term. |
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #789 +/- ##
==========================================
- Coverage 93.30% 88.96% -4.35%
==========================================
Files 36 36
Lines 5257 5382 +125
==========================================
- Hits 4905 4788 -117
- Misses 352 594 +242 ☔ View full report in Codecov by Sentry. |
Overall, from what I can read, vertices and corners are synonyms, but yet I wouldn't call the end-points of a line for corners, but rather vertices? |
|
I suggest we use vertex then, and we can include a glossar in the docs with clarifications if needed. |
|
Knut pointed out that the dispatches are technically wrong, since facets is defined based on codimension of I also think we can get away with not renaming
In this way, we would not need to introduce the concept of "facet" to the user, it would only be an internal thing. |
In the following case, I would expect # Ferrite grid: Embedded line element is nr 2
# 4 ____ 3 3 3 ____ 6
# | | | | |
# | (1) | | | (3) |
# 1 ____ 2 2 2 ____ 5I would therefore argue that we should change the fields in
I think that while doing this change now, we should be fully consistent with the notation and avoid "exceptions", such as saying that
In this case, I would expect that I see that in this setting, it would also be relevant to know the edge tangent direction though, but this would be similar wrt. an embedded I took a quick look at the code, and a change that could make sense is to include a direction-field in
|
But surely you couldn't treat these three entities the same in the integration? You would in this case do something for the edges and something else for the vertex, right? I think we should just not allow |
True, I would have different loops for these. But if I ask for the set on the boundary, I expect both to be included, and then I would split the set via intersection with the subdofhandlers. I think this is how it would currently be on the master, except there we would ask for the faceset instead of the facetset.
I don't see the problem for keeping this functionality when we rename FaceIndex to FacetIndex, this would all just be a |
The main reasons for doing this is :
if dim== 2etc. We can also remove the "hacks" we have for embedded elements (shells, beams). In other words we will have better integration of embedded elements.Another indication that this might be the right change is that many other FE-related packages uses this definition, e.g Finics, Gmsh, defelements (also deal.ii i think)
The drawback is:
getedgeset()instead ofgetfaceset()(actually, this is not strictly necessary because the we can keep backwards compatibility by defining:getfacset(grid::Grid{2}, name::String) = grid.edgesets[name].... however, it might be a bit strange thatgetfacesetreturns an edgeset (for 2d-problems)).One possible way of keeping a "dimension independent" way of accessing boundary-sets (facesets/edgeset) is to introduce "facets", for example
getfacetset, addfacetset, with the following meaning:To be clear, this is only a suggestion. Please come with feedback :)