Topology Bug

[termi-official]

> grid = generate_grid(Hexahedron, (2, 1, 1))
> topo = ExclusiveTopology(grid)
> getneighborhood(topo, grid, FaceIndex(2,6), true)

ERROR: BoundsError
Stacktrace:
 [1] getindex(A::SparseArrays.SparseMatrixCSC{Ferrite.EntityNeighborhood, Int64}, i0::Int64, i1::Int64)
   @ SparseArrays ~/Tools/julia-1.7.2/share/julia/stdlib/v1.7/SparseArrays/src/sparsematrix.jl:2159
 [2] getneighborhood(top::ExclusiveTopology, grid::Grid{3, Hexahedron, Float64}, faceidx::FaceIndex, include_self::Bool)
   @ Ferrite ~/.julia/packages/Ferrite/X5BVa/src/Grid/grid.jl:0
 [3] top-level scope
   @ REPL[15]:1

and

julia> grid = generate_grid(QuadraticQuadrilateral, (3, 1));

julia> topo = ExclusiveTopology(grid)
ExclusiveTopology(Dict(5 => [2, 3], 15 => [1], 7 => [3], 21 => [3], 17 => [1, 2], 3 => [1, 2], 1 => [1], 19 => [2, 3]), Ferrite.EntityNeighborhood{CellIndex}[Ferrite.EntityNeighborhood{CellIndex}(CellIndex[CellIndex(2)]), Ferrite.EntityNeighborhood{CellIndex}(CellIndex[CellIndex(1), CellIndex(3)]), Ferrite.EntityNeighborhood{CellIndex}(CellIndex[CellIndex(2)])], sparse(Int64[], Int64[], Ferrite.EntityNeighborhood[], 0, 0), sparse(Int64[], Int64[], Ferrite.EntityNeighborhood[], 0, 0), sparse(Int64[], Int64[], Ferrite.EntityNeighborhood[], 0, 0), Dict{Int64, Ferrite.EntityNeighborhood{VertexIndex}}(5 => Ferrite.EntityNeighborhood{VertexIndex}(VertexIndex[VertexIndex((2, 1)), VertexIndex((2, 3)), VertexIndex((3, 2)), VertexIndex((3, 4))]), 15 => Ferrite.EntityNeighborhood{VertexIndex}(VertexIndex[VertexIndex((1, 1)), VertexIndex((1, 3))]), 7 => Ferrite.EntityNeighborhood{VertexIndex}(VertexIndex[VertexIndex((3, 1)), VertexIndex((3, 3))]), 21 => Ferrite.EntityNeighborhood{VertexIndex}(VertexIndex[VertexIndex((3, 2)), VertexIndex((3, 4))]), 17 => Ferrite.EntityNeighborhood{VertexIndex}(VertexIndex[VertexIndex((1, 2)), VertexIndex((1, 4)), VertexIndex((2, 1)), VertexIndex((2, 3))]), 3 => Ferrite.EntityNeighborhood{VertexIndex}(VertexIndex[VertexIndex((1, 1)), VertexIndex((1, 3)), VertexIndex((2, 2)), VertexIndex((2, 4))]), 1 => Ferrite.EntityNeighborhood{VertexIndex}(VertexIndex[VertexIndex((1, 2)), VertexIndex((1, 4))]), 19 => Ferrite.EntityNeighborhood{VertexIndex}(VertexIndex[VertexIndex((2, 2)), VertexIndex((2, 4)), VertexIndex((3, 1)), VertexIndex((3, 3))])), FaceIndex[FaceIndex((1, 1)), FaceIndex((1, 2)), FaceIndex((1, 3)), FaceIndex((1, 4)), FaceIndex((2, 1)), FaceIndex((2, 2)), FaceIndex((2, 3)), FaceIndex((3, 1)), FaceIndex((3, 2)), FaceIndex((3, 3))])

julia> topo.
cell_neighbor          edge_neighbor           face_neighbor           face_skeleton           vertex_neighbor         vertex_to_cell          vertex_vertex_neighbor
julia> topo.face_neighbor
0×0 SparseArrays.SparseMatrixCSC{Ferrite.EntityNeighborhood, Int64} with 0 stored entries

[zkk960317]

I find another bug with function ExclusiveTopology at Line 193 and 194 in Grid.jl

neighbor_local_ids = findall(x->x in cell.nodes, cells[neighbor].nodes)
cell_local_ids = findall(x->x in cells[neighbor].nodes, cell.nodes)

This is not right for non-first-order cells, because the number of nodes shared by neighboring cells is greater than the number of vertexes. This causes a problem for the following conditional judgments.　I think the correct codes are:

neighbor_local_ids = findall(x->x in cell_vertices_table[cellid], cell_vertices_table[neighbor])
cell_local_ids = findall(x->x in cell_vertices_table[neighbor], cell_vertices_table[cellid])

@termi-official @koehlerson

[koehlerson]

you are absolutely right, I fixed it slightly different here:
https://github.com/Ferrite-FEM/Ferrite.jl/pull/455/files#diff-13ae190cb008d749905405bddca2a9ce395b303e12c781186df0d908efd05a90R195-R196

@zkk960317 could you try out that PR for your application and let me know if it works as expected?

[zkk960317]

@koehlerson Yeah, it seems that the result is correct for QuadraticTriangle cells. I'm not sure for other celltypes. But what I want to know is whether my code and your code have the same effect? Are there computational efficiency considerations involved?

[koehlerson]

Ok nice, I already added a regression test for quadratic quadrilaterals here https://github.com/Ferrite-FEM/Ferrite.jl/pull/455/files#diff-796a5ea6f5cdd76d08a6b0a07f1bd1813fa58607624b572ef71f1cd8ffdddf8aR363-R369. Probably I should add some for three dimensional cells as well.

I think our codes do the very same, but I think I should stick to your code suggestion, since it doesn't call again vertices but rather accesses the already built cell to vertex table with the correct indices :) I'll try it today in the evening, but from my gut feeling I would say your provided code snippet should perform better