ExclusiveTopology skeletons

[KnutAM]

While working with ExclusiveTopology in #974, I get the impression that we would only require a single facet_skeleton::Vector{FacetIndex}, instead of having separate vertex_skeleton, edge_skeleton, and face_skeleton fields saved. I don't want to do those changes there since unrelated, but unless there is anything I'm missing, I think simplifying that in a later PR would be good...

[termi-official]

This boils down to: Is there a cheap algorithm to deduplicate the edges and vertices, given the grid and facet skeleton. Can you sketch the algorithm which you have in mind to achieve this?

[KnutAM]

Currently, we only support ExclusiveTopology for grids where all cells have the same reference dimension, and then there is no need to de-duplicate as facet corresponds to vertex, edge, or face - same for all cells.

If we have a grid with cells of different reference dimension, what is actually the definition of any skeleton?
E.g. if we have

(7) --- (8) --- (9)
 |       |       |
 |  4    5    6  7
 |       |       |
(4) --- (5) --- (6)
 |       |       |
 |  1    2    3  |
 |       |       |
(1) --- (2) --- (3)

• Is edge (6,9) part of the edge skeleton? (shared between cell 6 and 7)
• Which vertices are part of the vertex skeleton? (All but 1, 3, 7 are shared with at least one other cell)
• For interfaces, the edge (5-8) is part of 3 cells, how should this be handled (currently the code assumes only a single neighbor)

So to me it sounds like the skeleton functionality is tricky for mixed reference dimensions and only makes sense for facets between cells of the same reference dimension. Therefore, we would then rather error for mixed refdim grids on creating the skeleton, even after supporting topology for mixed refdims.

If the need to combine these two features would come in the future, I guess there will be solutions, which just requires more information from the user side. E.g. the user could pass a cellset or similar to the topology constructor.

[lijas]

In adaptivity with hexahedrons, don't you need to find hanging nodes on faces and edges, so only keeping track of facets would not be enough?

[termi-official]

Currently, we only support ExclusiveTopology for grids where all cells have the same reference dimension, and then there is no need to de-duplicate as facet corresponds to vertex, edge, or face - same for all cells.

Taking your example grid above, if I iterate over the edge skeleton, how do I know that I have visited an edge already only with the information with in EdgeIndex?

If we have a grid with cells of different reference dimension, what is actually the definition of any skeleton?

A set of unique entities. E.g. for the face skeleton a set with all unique faces.

In adaptivity with hexahedrons, don't you need to find hanging nodes on faces and edges, so only keeping track of facets would not be enough?

The algorithms for adaptivity indeed need extra information, which is stored in the corresponding data structures. See e.g.

Ferrite.jl/src/Adaptivity/BWG.jl

Lines 795 to 919 in b913af4

	"""
	hangingnodes(forest,nodeids,nodeowners)
	Constructs a map from constrained nodeids to the ones that constraint
	"""
	function hangingnodes(forest::ForestBWG{dim}, nodeids, nodeowners) where dim
	_perm = dim == 2 ? 𝒱₂_perm : 𝒱₃_perm
	_perminv = dim == 2 ? 𝒱₂_perm_inv : 𝒱₃_perm_inv
	facetable = dim == 2 ? 𝒱₂ : 𝒱₃
	opposite_face = dim == 2 ? opposite_face_2 : opposite_face_3
	hnodes = Dict{Int,Vector{Int}}()
	facet_neighborhood = Ferrite.Ferrite.get_facet_facet_neighborhood(forest)
	for (k,tree) in enumerate(forest.cells)
	rootfaces = faces(root(dim),tree.b)
	for (l,leaf) in enumerate(tree.leaves)
	c̃ = child_id(leaf,tree.b)
	if leaf == root(dim)
	continue
	end
	for (ci,c) in enumerate(vertices(leaf,tree.b))
	parent_ = parent(leaf,tree.b)
	parentfaces = faces(parent_,tree.b)
	for (pface_i, pface) in enumerate(parentfaces)
	if iscenter(c,pface) #hanging node candidate
	neighbor_candidate = facet_neighbor(parent_, pface_i, tree.b)
	if inside(tree,neighbor_candidate) #intraoctree branch
	neighbor_candidate_idx = findfirst(x->x==neighbor_candidate,tree.leaves)
	if neighbor_candidate_idx !== nothing
	neighbor_candidate_faces = faces(neighbor_candidate,tree.b)
	nf = findfirst(x->x==pface,neighbor_candidate_faces)
	hnodes[nodeids[nodeowners[(k,c)]]] = [nodeids[nodeowners[(k,nc)]] for nc in neighbor_candidate_faces[nf]]
	if dim > 2
	vs = vertices(leaf,tree.b)
	for ξ ∈ 1:ncorners_face3D
	c′ = facetable[pface_i, ξ]
	if c′ ∉ (c̃,ci)
	neighbor_candidate_edges = edges(neighbor_candidate,tree.b)
	ne = findfirst(x->iscenter(vs[c′],x),neighbor_candidate_edges)
	if ne !== nothing
	hnodes[nodeids[nodeowners[(k,vs[c′])]]] = [nodeids[nodeowners[(k,ne)]] for ne in neighbor_candidate_edges[ne]]
	end
	end
	end
	end
	break
	end
	else #interoctree branch
	for (ri,rf) in enumerate(rootfaces)
	facet_neighbor_ = facet_neighborhood[k,_perm[ri]]
	if length(facet_neighbor_) == 0
	continue
	end
	if contains_facet(rf, pface)
	k′ = facet_neighbor_[1][1]
	ri′ = _perminv[facet_neighbor_[1][2]]
	interoctree_neighbor = transform_facet(forest, k′, ri′, neighbor_candidate)
	interoctree_neighbor_candidate_idx = findfirst(x->x==interoctree_neighbor,forest.cells[k′].leaves)
	if interoctree_neighbor_candidate_idx !== nothing
	r = compute_face_orientation(forest,k,pface_i)
	neighbor_candidate_faces = faces(neighbor_candidate,forest.cells[k′].b)
	transformed_neighbor_faces = faces(interoctree_neighbor,forest.cells[k′].b)
	fnodes = transformed_neighbor_faces[ri′]
	vs = vertices(leaf,tree.b)
	if dim > 2
	rotated_ξ = ntuple(i->rotation_permutation(ri′,ri,r,i),ncorners_face3D)
	else
	rotated_ξ = ntuple(i->rotation_permutation(r,i),ncorners_face2D)
	end
	hnodes[nodeids[nodeowners[(k,c)]]] = [nodeids[nodeowners[(k′,fnodes[ξ])]] for ξ in rotated_ξ]
	if dim > 2
	for ξ in rotated_ξ
	c′ = facetable[pface_i, ξ]
	if c′ ∉ (c̃,c)
	neighbor_candidate_edges = edges(interoctree_neighbor,tree.b)
	ne = findfirst(x->iscenter(vs[c′],x),neighbor_candidate_edges)
	if ne !== nothing
	hnodes[nodeids[nodeowners[(k,vs[c′])]]] = [nodeids[nodeowners[(k,ne)]] for ne in neighbor_candidate_edges[ne]]
	end
	end
	end
	end
	break
	end
	end
	end
	end
	end
	end
	if dim > 2
	parentedges = edges(parent_,tree.b)
	for (pedge_i, pedge) in enumerate(parentedges)
	if iscenter(c,pedge) #hanging node candidate
	neighbor_candidate = edge_neighbor(parent_, pedge_i, tree.b)
	for (ri,re) in enumerate(edges(root(dim),tree.b))
	edge_neighbor_ = forest.topology.edge_edge_neighbor[k,edge_perm[ri]]
	if length(edge_neighbor_) == 0
	continue
	end
	if contains_edge(re, pedge)
	k′ = edge_neighbor_[1][1]
	ri′ = edge_perm_inv[edge_neighbor_[1][2]]
	interoctree_neighbor = transform_edge(forest, k′, ri′, neighbor_candidate, true)
	interoctree_neighbor_candidate_idx = findfirst(x->x==interoctree_neighbor,forest.cells[k′].leaves)
	if interoctree_neighbor_candidate_idx !== nothing
	neighbor_candidate_edges = edges(neighbor_candidate,forest.cells[k′].b)
	transformed_neighbor_edges = edges(interoctree_neighbor,forest.cells[k′].b)
	ne = findfirst(x->iscenter(c,x),neighbor_candidate_edges)
	if ne !== nothing
	hnodes[nodeids[nodeowners[(k,c)]]] = [nodeids[nodeowners[(k′,nc)]] for nc in transformed_neighbor_edges[ne]]
	else
	@error "things are messed up for hanging edge constraint interoctree at octree $k edge $(_perm[ri])"
	end
	break
	end
	end
	end
	end
	end
	end
	end
	end
	end
	return hnodes
	end

[KnutAM]

Actually, it seems like vertexskeleton really doesn't work as intended.

grid = generate_grid(Quadrilateral, (2,1));
Ferrite.vertexskeleton(ExclusiveTopology(grid), grid) # Gives 8-element Vector{VertexIndex}:

whereas there are only 6 unique vertices in this grid. This should be fixed!

A set of unique entities. E.g. for the face skeleton a set with all unique faces.

Ah, should have known, I got confused since entities without neighbors are skipped by the InterfaceIterator. I guess the problem is rather there then, and I simply assumed that this was the only use for the skeleton.

Sidetrack, but is it only me who finds skeleton confusing? (It makes sense for facets, but for vertices these are not connected and its sounds strange.) Are unique_faces, unique_edges, and unique_vertices perhaps better?

[termi-official]

Actually, it seems like vertexskeleton really doesn't work as intended.

grid = generate_grid(Quadrilateral, (2,1));
Ferrite.vertexskeleton(ExclusiveTopology(grid), grid) # Gives 8-element Vector{VertexIndex}:

whereas there are only 6 unique vertices in this grid. This should be fixed!

Where does this function even come from? I can only remember providing the face skeleton for DG.

A set of unique entities. E.g. for the face skeleton a set with all unique faces.

Ah, should have known, I got confused since entities without neighbors are skipped by the InterfaceIterator. I guess the problem is rather there then, and I simply assumed that this was the only use for the skeleton.

I do not understand the problem here. Can you elaborate?

Sidetrack, but is it only me who finds skeleton confusing? (It makes sense for facets, but for vertices these are not connected and its sounds strange.) Are unique_faces, unique_edges, and unique_vertices perhaps better?

Skeleton is the name used in literature, so we sticked to it while implementing it.

[KnutAM]

Where does this function even come from? I can only remember providing the face skeleton for DG.

This was included in the massive face -> facet pr, because when DG was implemented, facets were called faces, and yes, then it was probably only implemented for faces (now facets).
This is part of the background as to why I was asking if only providing it for facets would suffice, since that seemed to be the original intention...

Skeleton is the name used in literature, so we sticked to it while implementing it.

Makes sense!

I do not understand the problem here. Can you elaborate?

The InterfaceIterator only supports the case of single refdim, hence to me it made sense that the skeleton should also only support single refdim (and then facetskeleton should be enough). But since that errors for mixed refdim, I guess this is ok for now.

[termi-official]

Where does this function even come from? I can only remember providing the face skeleton for DG.

This was included in the massive face -> facet pr, because when DG was implemented, facets were called faces, and yes, then it was probably only implemented for faces (now facets). This is part of the background as to why I was asking if only providing it for facets would suffice, since that seemed to be the original intention...

Ah, I see. Yes, for now we mostly use the facet skeleton for operations on the trace space (e.g. interface integrals in DG). I know also use-cases where the edge skeleton or the vertex skeleton is needed, but these are really niche (e.g. interaction with generalized finite difference stencils or anisotropic preconditioning) so I think it is okay to just remove the edge and vertex skeleton for now.

I do not understand the problem here. Can you elaborate?

The InterfaceIterator only supports the case of single refdim, hence to me it made sense that the skeleton should also only support single refdim (and then facetskeleton should be enough). But since that errors for mixed refdim, I guess this is ok for now.

Got it, thanks.