ExclusiveTopology with ArrayOfVectorViews backend

[KnutAM]

See discussion in #965, this PR introduces a cleaned-up version.

Most important changes in this PR compared to master are:

• Storage from Array{EntityNeighborhood{T}, N} -> ArrayOfVectorViews{T, N} <: AbstractArray{"View{Vector{T}}", N}
• getcells dispatches removed (since these would now act on Vector{<:BoundaryIndex}). In any case, getcells normally returns ::AbstractCells, but in topology they return ::Ints on mater, if any, cellid should be overloaded).
• EntityNeighborhood is gone

Compared to #965, the support for hashed indexing is removed and only ExclusiveTopology is implemented. Before finalizing, I would like some feedback from the topology side, and then I'll add

• Separate tests for ArrayOfVectorViews
• Documentation of ArrayOfVectorViews
• Add warning for non-conforming grids, ref fix topology getneighborhood call in 2D #941 (comment) (and for embedded grids)
• Update docstrings for ExclusiveTopology
• Changelog

Let me know if more details are required for the initial review.

Note that in relation to #942, this PR explicitly errors if trying to construct ExclusiveTopology with grids containing embedded cells, thus avoiding silent bugs.

Benchmarks

using Ferrite, BenchmarkTools
function sum_indices(nbh::AbstractMatrix)
    v = 0
    for nbs in nbh
        for n in nbs
            v += n[1] + n[2]
        end
    end
    return v
end

function sum_indices(top::ExclusiveTopology)
    v  = sum_indices(top.vertex_vertex_neighbor)
    v += sum_indices(top.edge_edge_neighbor)
    v += sum_indices(top.face_face_neighbor)
    return v
end
function sum_skeleton_facet_neighbors(top, grid)
    skeleton = Ferrite.facetskeleton(top, grid)
    v = 0
    for facetidx in skeleton
        v += sum(n[1] + n[2] for n in getneighborhood(top, grid, facetidx); init = 0)
    end
    return v
end

grid = generate_grid(Hexahedron, (10, 10, 10));
top = ExclusiveTopology(grid); # Compilation

@btime sum_indices($top); 
# Master: 55.700 μs (0 allocations: 0 bytes)
#     PR: 36.000 μs (0 allocations: 0 bytes)
Ferrite.facetskeleton(top, grid); # Compilation

grid = generate_grid(Hexahedron, (100, 100, 100));
@time top = ExclusiveTopology(grid);
# Master: 27.150434 seconds (64.58 M allocations: 7.195 GiB, 72.43% gc time)
#     PR:  3.799651 seconds (45 allocations: 1.895 GiB, 13.84% gc time)
@time Ferrite.facetskeleton(top, grid);
# Master: 0.143506 seconds (2 allocations: 46.234 MiB)
#     PR: 0.032492 seconds (2 allocations: 46.234 MiB)
@btime sum_skeleton_facet_neighbors($top, $grid);
# Master: 59.683 ms (0 allocations: 0 bytes)
#     PR: 21.887 ms (0 allocations: 0 bytes)

(PR refers to b7e2add)

[codecov]

Codecov Report

Attention: Patch coverage is 98.99497% with 2 lines in your changes missing coverage. Please review.

Project coverage is 93.58%. Comparing base (2c8b751) to head (2961031).
Report is 1 commits behind head on master.

Files	Patch %	Lines
src/Grid/topology.jl	98.57%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #974      +/-   ##
==========================================
+ Coverage   93.20%   93.58%   +0.37%     
==========================================
  Files          38       39       +1     
  Lines        5816     5888      +72     
==========================================
+ Hits         5421     5510      +89     
+ Misses        395      378      -17     

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

[KnutAM]

c276cb2 brings construction in benchmark further down from
5.527448 seconds (58 allocations: 2.107 GiB, 9.91% gc time)
to
4.217167 seconds (45 allocations: 1.895 GiB, 12.49% gc time)

[KnutAM]

When rechecking the timings, I saw some worse performance now.
So took @termi-official's suggestion and optimized the sizehints, benchmark results are updated with this optimization.

Script

function analyze_top_sizes(grid; show_avg=false)
    top = ExclusiveTopology(grid);
    for key in fieldnames(typeof(top))
        nbh = getfield(top, key)
        nbh isa Ferrite.ArrayOfVectorViews || continue
        max_size = maximum(length, nbh; init=0)
        avg_size = sum(length, nbh; init=0)/length(nbh)
        if show_avg
            println(key)
            println("  max: ", max_size, " entities")
            println("  avg: ", avg_size, " entities")
        else
            println(key, ": ", max_size)
        end
    end
end

function print_sizes(n=10, celltypes=[Line, Triangle, Quadrilateral, Tetrahedron, Hexahedron])
   for CT in celltypes
       grid = generate_grid(CT, ntuple(_->n, Ferrite.getrefdim(CT)))
       println(CT)
       analyze_top_sizes(grid)
       println()
   end
end

Max number of entities for different generated grids

field	Line	Triangle	Quadrilateral	Tetrahedron	Hexahedron
vertex_to_cell	2	6	4	24	8
cell_neighbor	2	12	8	70	26
face_face_neighbor	0	0	0	1	1
edge_edge_neighbor	0	1	1	3	1
vertex_vertex_neighbor	1	3	1	13	1

[KnutAM]

I will fix #988 after merging this PR to avoid annoying conflicts.

[termi-official]

Resolves #617 ?

[termi-official]

Great improvement on the performance and I think the data structure will also help in the long run. I left a few comments.

[KnutAM]

I left a few comments.

Thanks!

[termi-official]

LGTM. Can you dump the benchmark code in the benchmark folder (here or separate PR)?

[KnutAM]

LGTM. Can you dump the benchmark code in the benchmark folder (here or separate PR)?

Thanks! I haven't looked so much at the benchmark setup, will just open an issue for now to not forget to do it later unless someone solves it before I have time.