add dispatch for 2D edgeidx topology

[koehlerson]

Fixes #941. I don't know how to test if the correct codepath is tracked, since the outcome is the same for both code paths.

[codecov]

Codecov Report

Attention: Patch coverage is 80.00000% with 1 line in your changes missing coverage. Please review.

Project coverage is 93.73%. Comparing base (6dfa0e4) to head (2fafe56).

Files	Patch %	Lines
src/Grid/topology.jl	80.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #942      +/-   ##
==========================================
- Coverage   93.75%   93.73%   -0.02%     
==========================================
  Files          36       36              
  Lines        5441     5445       +4     
==========================================
+ Hits         5101     5104       +3     
- Misses        340      341       +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[KnutAM]

Is this more than an optimization for grids with sdim=2 which allows assuming that each edge is only shared between 2 cells?

I.e. how does this play with embedded shells, when we have AbstractGrid{3} but rdim=2? Would that still lead to the same error that you had for the AMR?

[koehlerson]

In future aspects, this is also a bugfix. However, we don't support currently mixed dimensional types in the topology. This was excluded for performance reasons, but can be brought back by adding mixed dimension tables #843

[termi-official]

In future aspects, this is also a bugfix. However, we don't support currently mixed dimensional types in the topology. This was excluded for performance reasons, but can be brought back by adding mixed dimension tables #843

Technically speaking, the latest changes should allow type stable queries on mixed dimensional grids.

[KnutAM]

In future aspects, this is also a bugfix.

I'm not following this though, because with the changes in #914, the topology is now completely independent of the spatial dimension of the grid. So AFAIU, the only point of this change is to speed up the computation when we can assume that each edge is only shared between two cells.
This assumption would be broken in the future when supporting embedded elements, since we could then have a membrane (Line) element B(2-3) like this

4 -- 3 3 3 -- 6
| A  | | | C  |
1 -- 2 2 2 -- 5

(The same would of course apply to faces, and that needs to be generalized to work on mixed grids, AFAIU the opposite of what this PR does to faces)

Actually, the above example works on master (shown below) but gives the wrong neighborhood with this branch

julia> nodes = [Node(Float64.(x)) for x in ((0,0), (1,0), (1,1), (0,1), (2,0), (2,1))];

julia> cells = [Quadrilateral((1, 2, 3, 4)), Line((2, 3)), Quadrilateral((2, 5, 6, 3))];

julia> grid = Grid(cells, nodes);

julia> top = ExclusiveTopology(grid);

julia> getneighborhood(top, grid, EdgeIndex(2,1), true)
3-element Vector{EdgeIndex}:
 EdgeIndex((3, 4))
 EdgeIndex((1, 2))
 EdgeIndex((2, 1))

julia> getneighborhood(top, grid, EdgeIndex(1,2), true)
3-element Vector{EdgeIndex}:
 EdgeIndex((3, 4))
 EdgeIndex((2, 1))
 EdgeIndex((1, 2))

# By purpose not allowed for grids with mixed reference dimensions:
julia> getneighborhood(top, grid, FacetIndex(1,2), true)
ERROR: ArgumentError: getneighborhood with FacetIndex is is only supported for grids containing cells with a common reference dimension.
...

[KnutAM]

What could be possible though, is instead of dispatching on the sdim would be to querry the reference dimension of the grid, something like this perhaps

function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, edgeidx::EdgeIndex, include_self=false)
    if get_reference_dimension(grid) == 2
        if include_self
            return [top.edge_edge_neighbor[edgeidx[1],edgeidx[2]].neighbor_info; edgeidx]
        else
            return top.edge_edge_neighbor[edgeidx[1],edgeidx[2]].neighbor_info
        end
    else
        # generic code as before
    end
end

get_reference_dimension will return :mixed if there are different reference dimensions. However, it is only "free of cost" to call for concrete cell types (or union cell types with the same reference dimension). So in that case, it might make more sense to call this during ExclusiveTopology construction, and save it.
(For type stability could perhaps save it as Symbol(get_reference_dimension)) Edit: Ugly, Union{Symbol, Int} is probably better.

[koehlerson]

In future aspects, this is also a bugfix.

I'm not following this though, because with the changes in #914, the topology is now completely independent of the spatial dimension of the grid. So AFAIU, the only point of this change is to speed up the computation when we can assume that each edge is only shared between two cells. This assumption would be broken in the future when supporting embedded elements, since we could then have a membrane (Line) element B(2-3) like this

4 -- 3 3 3 -- 6
| A  | | | C  |
1 -- 2 2 2 -- 5

(The same would of course apply to faces, and that needs to be generalized to work on mixed grids, AFAIU the opposite of what this PR does to faces)

Actually, the above example works on master (shown below) but gives the wrong neighborhood with this branch

julia> nodes = [Node(Float64.(x)) for x in ((0,0), (1,0), (1,1), (0,1), (2,0), (2,1))];

julia> cells = [Quadrilateral((1, 2, 3, 4)), Line((2, 3)), Quadrilateral((2, 5, 6, 3))];

julia> grid = Grid(cells, nodes);

julia> top = ExclusiveTopology(grid);

julia> getneighborhood(top, grid, EdgeIndex(2,1), true)
3-element Vector{EdgeIndex}:
 EdgeIndex((3, 4))
 EdgeIndex((1, 2))
 EdgeIndex((2, 1))

julia> getneighborhood(top, grid, EdgeIndex(1,2), true)
3-element Vector{EdgeIndex}:
 EdgeIndex((3, 4))
 EdgeIndex((2, 1))
 EdgeIndex((1, 2))

# By purpose not allowed for grids with mixed reference dimensions:
julia> getneighborhood(top, grid, FacetIndex(1,2), true)
ERROR: ArgumentError: getneighborhood with FacetIndex is is only supported for grids containing cells with a common reference dimension.
...

That example works because of the happy coincidence that the code branch does the correct thing even though it was intended for something different. This doesn't mean that mixed dimensional grids work in general on master.

From the test suite (currently commented out due to unsupport of mixed dimensional grids)

  cells = [
        Hexahedron((1, 2, 3, 4, 5, 6, 7, 8)),
        Hexahedron((11, 13, 14, 12, 15, 16, 17, 18)),
        Quadrilateral((2, 9, 10, 3)),
        Quadrilateral((9, 11, 12, 10)),
        ]
    nodes = [Node(coord) for coord in zeros(Vec{3,Float64}, 18)]
    grid = Grid(cells, nodes)
    topology = ExclusiveTopology(grid)
    getneighborhood(top,grid,FaceIndex(3,4)) # returns nothing, should be EdgeIndex(1,2)

[koehlerson]

If you are sure that the code branch always does the correct thing for every possible mixed dimensional grid we want to support, then it's okay to leave it as is and close the PR. I'm just not sure that's the case since I wrote it with something else in mind

[KnutAM]

If you are sure that the code branch always does the correct thing for every possible mixed dimensional grid we want to support, then it's okay to leave it as is and close the PR. I'm just not sure that's the case since I wrote it with something different in mind

No, I'm not sure :) However, it seems to me (with the above example) that this PR doesn't make it more correct, or do I miss something? As I understand, the 3d code does not assume that there is only one edge, while the 2d code does.

Another option, for solving the issue with AMR, could even be to dispatch only the case with FacetIndex to the special branch, since when called with FacetIndex, it is guaranteed to not be embedded mixed.

[koehlerson]

It makes it more correct in the sense of restoring exactly the code branch that has been called before #914 I will go in the AMR branch probably with get_facet_facet_neighborhood since I just need the tables and nothing reconstructed there

[KnutAM]

It makes it more correct in the sense of restoring exactly the code branch that has been called before #914 I will go in the AMR branch probably with get_facet_facet_neighborhood since I just need the tables and nothing reconstructed there.

In that case we should error on ExclusiveTopology (as before 914) until #843 is fixed, which sounds like a good solution to me. With that additional change, I think this PR good! But of course, considering the discussion I wouldn't count on it to stay after 843, but I could be wrong here!

Concretely, I suggest changing

Ferrite.jl/src/Grid/topology.jl

Line 247 in cd82d86

ExclusiveTopology(grid::AbstractGrid) = ExclusiveTopology(getcells(grid))

to

function ExclusiveTopology(grid::AbstractGrid)
    get_reference_dimension(grid) == getspatialdim(grid) || throw(ArgumentError("Exclusive topology is currently not supported for grids containing cells with different reference dimensions"))
    return ExclusiveTopology(getcells(grid)) 
end

[KnutAM]

Sorry for misunderstanding this earlier, after working more with the topology, I think I understand the problem :)

In #974 I enforce rdim=sdim for all cells, so IMO we can skip that here and do it there instead, and this can from my perspective be merged as-is!