diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index f734ea2416..71e57d643e 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -568,7 +568,8 @@ function sortedge(edge::Tuple{Int,Int}) end """ - sortface(face::Tuple{Int,Int}) + sortface(face::Tuple{Int}) + sortface(face::Tuple{Int,Int}) sortface(face::Tuple{Int,Int,Int}) sortface(face::Tuple{Int,Int,Int,Int}) @@ -626,6 +627,7 @@ function sortface(face::Tuple{Int,Int,Int,Int}) return (a, b, c), SurfaceOrientationInfo() # TODO fill struct end +sortface(face::Tuple{Int}) = face, nothing """ find_field(dh::DofHandler, field_name::Symbol)::NTuple{2,Int} diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index b95f79580b..d9533f973b 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -90,7 +90,7 @@ function vertices(c::AbstractCell{1, RefLine}) end function faces(c::AbstractCell{1, RefLine}) ns = get_node_ids(c) - return (ns[1], ns[2]) # f1, f2 + return ((ns[1],), (ns[2],)) # f1, f2 end # RefTriangle (refdim = 2): vertices for vertexdofs, faces for facedofs (edgedofs) and BC @@ -259,19 +259,15 @@ function Base.show(io::IO, ::MIME"text/plain", n::EntityNeighborhood) end """ - face_npoints(::AbstractCell) -Specifies for each subtype of AbstractCell how many nodes form a face + nvertices_on_face(cell::AbstractCell, local_face_index::Int) +Specifies for each subtype of AbstractCell how many nodes form a face. """ -face_npoints(::AbstractCell{2}) = 2 -# face_npoints(::Cell{3,4,1}) = 4 #not sure how to handle embedded cells e.g. Quadrilateral3D +nvertices_on_face(cell::AbstractCell, local_face_index::Int) = length(faces(cell)[local_face_index]) """ - edge_npoints(::AbstractCell) -Specifies for each subtype of AbstractCell how many nodes form an edge + nvertices_on_edge(::AbstractCell, local_edge_index::Int) +Specifies for each subtype of AbstractCell how many nodes form an edge. """ -# edge_npoints(::Cell{3,4,1}) = 2 #not sure how to handle embedded cells e.g. Quadrilateral3D -face_npoints(::AbstractCell{3, RefHexahedron}) = 4 -face_npoints(::AbstractCell{3, RefTetrahedron}) = 3 -edge_npoints(::AbstractCell{3}) = 2 +nvertices_on_edge(cell::AbstractCell, local_edge_index::Int) = length(edges(cell)[local_edge_index]) getdim(::Union{AbstractCell{refdim},Type{<:AbstractCell{refdim}}}) where {refdim} = refdim @@ -280,8 +276,7 @@ abstract type AbstractTopology end """ ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell `ExclusiveTopology` saves topological (connectivity) data of the grid. The constructor works with an `AbstractCell` -vector for all cells that dispatch `vertices`, `faces` and in 3D `edges` as well as the utility functions -`face_npoints` and `edge_npoints`. +vector for all cells that dispatch `vertices`, `faces` and in 3D `edges`. The struct saves the highest dimensional neighborhood, i.e. if something is connected by a face and an edge only the face neighborhood is saved. The lower dimensional neighborhood is recomputed, if needed. @@ -296,7 +291,7 @@ The struct saves the highest dimensional neighborhood, i.e. if something is conn """ struct ExclusiveTopology <: AbstractTopology # maps a global vertex id to all cells containing the vertex - vertex_to_cell::Dict{Int,Vector{Int}} + vertex_to_cell::Dict{Int,Set{Int}} # index of the vector = cell id -> all other connected cells cell_neighbor::Vector{EntityNeighborhood{CellIndex}} # face_neighbor[cellid,local_face_id] -> exclusive connected entities (not restricted to one entity) @@ -313,14 +308,14 @@ end function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell cell_vertices_table = vertices.(cells) #needs generic interface for <: AbstractCell - vertex_cell_table = Dict{Int,Vector{Int}}() - - for (cellid, cell_nodes) in enumerate(cell_vertices_table) - for node in cell_nodes - if haskey(vertex_cell_table, node) - push!(vertex_cell_table[node], cellid) + vertex_cell_table = Dict{Int,Set{Int}}() + + for (cellid, cell_vertices) in enumerate(cell_vertices_table) + for vertex in cell_vertices + if haskey(vertex_cell_table, vertex) + push!(vertex_cell_table[vertex], cellid) else - vertex_cell_table[node] = [cellid] + vertex_cell_table[vertex] = Set([cellid]) end end end @@ -330,29 +325,51 @@ function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell I_vertex = Int[]; J_vertex = Int[]; V_vertex = EntityNeighborhood[] cell_neighbor_table = Vector{EntityNeighborhood{CellIndex}}(undef, length(cells)) - for (cellid, cell) in enumerate(cells) - #cell neighborhood - cell_neighbors = getindex.((vertex_cell_table,), cell_vertices_table[cellid]) # cell -> vertex -> cell - cell_neighbors = unique(reduce(vcat,cell_neighbors)) # non unique list initially - filter!(x->x!=cellid, cell_neighbors) # get rid of self neighborhood - cell_neighbor_table[cellid] = EntityNeighborhood(CellIndex.(cell_neighbors)) + for (cellid, cell) in enumerate(cells) + cell_neighbors = reduce(union!, [Set{Int}(vertex_cell_table[vertex]) for vertex ∈ vertices(cell) if vertex_cell_table[vertex] != cellid]) + cell_neighbor_table[cellid] = EntityNeighborhood(CellIndex.(collect(cell_neighbors))) - for neighbor in cell_neighbors - 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]) + face_neighbors = Set{Int}() + for (face_idx,face) ∈ enumerate(faces(cell)) + neighbor_candidates = Set{Int}(c for c ∈ vertex_cell_table[face[1]] if c != cellid) + for face_vertex ∈ face[2:end] + intersect!(neighbor_candidates, vertex_cell_table[face_vertex]) + end + union!(face_neighbors, neighbor_candidates) + end + + if getdim(cell) > 2 + edge_neighbors = Set{Int}() + for (edge_idx,edge) ∈ enumerate(edges(cell)) + neighbor_candidates = Set{Int}(c for c ∈ vertex_cell_table[edge[1]] if c != cellid) + for edge_vertex ∈ edge[2:end] + edge_neighbor = vertex_cell_table[edge_vertex] + if edge_neighbor != cellid && edge_neighbor ∉ face_neighbors + intersect!(neighbor_candidates, edge_neighbor) + end + end + union!(edge_neighbors, neighbor_candidates) + end + end + + for neighbor_cellid in cell_neighbors + cell_local_ids = findall(x->x in cell_vertices_table[neighbor_cellid], cell_vertices_table[cellid]) # vertex neighbor if length(cell_local_ids) == 1 - _vertex_neighbor!(V_vertex, I_vertex, J_vertex, cellid, cell, neighbor_local_ids, neighbor, cells[neighbor]) + neighbor_local_ids = findall(x->x in cell_vertices_table[cellid], cell_vertices_table[neighbor_cellid]) + _vertex_neighbor!(V_vertex, I_vertex, J_vertex, cellid, cell, neighbor_local_ids, neighbor_cellid, cells[neighbor_cellid]) # face neighbor - elseif length(cell_local_ids) == face_npoints(cell) - _face_neighbor!(V_face, I_face, J_face, cellid, cell, neighbor_local_ids, neighbor, cells[neighbor]) + elseif neighbor_cellid ∈ face_neighbors + neighbor_local_ids = findall(x->x in cell_vertices_table[cellid], cell_vertices_table[neighbor_cellid]) + _face_neighbor!(V_face, I_face, J_face, cellid, cell, neighbor_local_ids, neighbor_cellid, cells[neighbor_cellid]) # edge neighbor - elseif getdim(cell) > 2 && length(cell_local_ids) == edge_npoints(cell) - _edge_neighbor!(V_edge, I_edge, J_edge, cellid, cell, neighbor_local_ids, neighbor, cells[neighbor]) + elseif getdim(cell) > 2 && neighbor_cellid ∈ edge_neighbors + neighbor_local_ids = findall(x->x in cell_vertices_table[cellid], cell_vertices_table[neighbor_cellid]) + _edge_neighbor!(V_edge, I_edge, J_edge, cellid, cell, neighbor_local_ids, neighbor_cellid, cells[neighbor_cellid]) end - end + end end - + celltype = eltype(cells) if isconcretetype(celltype) dim = getdim(cells[1]) @@ -402,9 +419,11 @@ function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell vertex_neighbors_local = VertexIndex[] vertex_neighbors_global = Int[] for cell in cellset - neighbor_boundary = getdim(cells[cell]) == 2 ? [faces(cells[cell])...] : [edges(cells[cell])...] #get lowest dimension boundary + neighbor_boundary = getdim(cells[cell]) > 2 ? collect(edges(cells[cell])) : collect(faces(cells[cell])) #get lowest dimension boundary neighbor_connected_faces = neighbor_boundary[findall(x->global_vertexid in x, neighbor_boundary)] - neighbor_vertices_global = getindex.(neighbor_connected_faces, findfirst.(x->x!=global_vertexid,neighbor_connected_faces)) + other_vertices = findfirst.(x->x!=global_vertexid,neighbor_connected_faces) + any(other_vertices .=== nothing) && continue + neighbor_vertices_global = getindex.(neighbor_connected_faces, other_vertices) neighbor_vertices_local= [VertexIndex(cell,local_vertex) for local_vertex in findall(x->x in neighbor_vertices_global, vertices(cells[cell]))] append!(vertex_neighbors_local, neighbor_vertices_local) append!(vertex_neighbors_global, neighbor_vertices_global) diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 439306a716..000f0a8129 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/test/test_grid_dofhandler_vtk.jl @@ -193,7 +193,6 @@ end end @testset "Grid sets" begin - grid = Ferrite.generate_grid(Hexahedron, (1, 1, 1), Vec((0.,0., 0.)), Vec((1.,1.,1.))) #Test manual add @@ -216,6 +215,18 @@ end end @testset "Grid topology" begin +# +# (1) (2) (3) (4) +# +---+---+---+ +# + linegrid = generate_grid(Line,(3,)) + linetopo = ExclusiveTopology(linegrid) + @test linetopo.vertex_neighbor[1,2] == Ferrite.EntityNeighborhood(VertexIndex(2,1)) + @test linetopo.vertex_neighbor[2,1] == Ferrite.EntityNeighborhood(VertexIndex(1,2)) + @test linetopo.vertex_neighbor[2,2] == Ferrite.EntityNeighborhood(VertexIndex(3,1)) + @test linetopo.vertex_neighbor[3,1] == Ferrite.EntityNeighborhood(VertexIndex(2,2)) + @test length(linetopo.face_skeleton) == 4 + # (11) # (10)+-----+-----+(12) # | 5 | 6 | @@ -272,9 +283,9 @@ end hexgrid = generate_grid(Hexahedron,(2,2,1)) topology = ExclusiveTopology(hexgrid) @test topology.edge_neighbor[1,11] == Ferrite.EntityNeighborhood(EdgeIndex(4,9)) - @test getneighborhood(topology,hexgrid,EdgeIndex(1,11),true) == [EdgeIndex(4,9),EdgeIndex(2,12),EdgeIndex(3,10),EdgeIndex(1,11)] + @test Set(getneighborhood(topology,hexgrid,EdgeIndex(1,11),true)) == Set([EdgeIndex(4,9),EdgeIndex(2,12),EdgeIndex(3,10),EdgeIndex(1,11)]) @test topology.edge_neighbor[2,12] == Ferrite.EntityNeighborhood(EdgeIndex(3,10)) - @test getneighborhood(topology,hexgrid,EdgeIndex(2,12),true) == [EdgeIndex(3,10),EdgeIndex(1,11),EdgeIndex(4,9),EdgeIndex(2,12)] + @test Set(getneighborhood(topology,hexgrid,EdgeIndex(2,12),true)) == Set([EdgeIndex(3,10),EdgeIndex(1,11),EdgeIndex(4,9),EdgeIndex(2,12)]) @test topology.edge_neighbor[3,10] == Ferrite.EntityNeighborhood(EdgeIndex(2,12)) @test topology.edge_neighbor[4,9] == Ferrite.EntityNeighborhood(EdgeIndex(1,11)) @test getneighborhood(topology,hexgrid,FaceIndex((1,3))) == [FaceIndex((2,5))] @@ -328,6 +339,8 @@ end @test topology.edge_neighbor[1,1] == topology.edge_neighbor[1,3] == zero(Ferrite.EntityNeighborhood{FaceIndex}) @test topology.face_neighbor[3,1] == topology.face_neighbor[3,3] == zero(Ferrite.EntityNeighborhood{FaceIndex}) @test topology.face_neighbor[4,1] == topology.face_neighbor[4,3] == zero(Ferrite.EntityNeighborhood{FaceIndex}) + + # # +-----+-----+-----+ # | 7 | 8 | 9 | @@ -351,23 +364,23 @@ end @test issubset([7,4,5,6,9], patches[8]) @test issubset([8,5,6], patches[9]) - @test Ferrite.getneighborhood(topology, quadgrid, VertexIndex(1,1)) == [VertexIndex(1,2), VertexIndex(1,4)] - @test Ferrite.getneighborhood(topology, quadgrid, VertexIndex(2,1)) == [VertexIndex(1,1), VertexIndex(1,3), VertexIndex(2,2), VertexIndex(2,4)] - @test Ferrite.getneighborhood(topology, quadgrid, VertexIndex(5,4)) == [VertexIndex(4,2), VertexIndex(4,4), VertexIndex(5,1), VertexIndex(5,3), VertexIndex(7,1), VertexIndex(7,3), VertexIndex(8,2), VertexIndex(8,4)] - @test Ferrite.getneighborhood(topology, quadgrid, VertexIndex(1,1),true) == [VertexIndex(1,2), VertexIndex(1,4), VertexIndex(1,1)] - @test Ferrite.getneighborhood(topology, quadgrid, VertexIndex(2,1),true) == [VertexIndex(1,1), VertexIndex(1,3), VertexIndex(2,2), VertexIndex(2,4), VertexIndex(1,2), VertexIndex(2,1)] - @test Ferrite.getneighborhood(topology, quadgrid, VertexIndex(5,4),true) == [VertexIndex(4,2), VertexIndex(4,4), VertexIndex(5,1), VertexIndex(5,3), VertexIndex(7,1), VertexIndex(7,3), VertexIndex(8,2), VertexIndex(8,4), VertexIndex(4,3), VertexIndex(5,4), VertexIndex(7,2), VertexIndex(8,1)] - @test Ferrite.toglobal(quadgrid, Ferrite.getneighborhood(topology, quadgrid, VertexIndex(1,1))) == [2,5] - @test Ferrite.toglobal(quadgrid, Ferrite.getneighborhood(topology, quadgrid, VertexIndex(2,1))) == [1,6,3] - @test Ferrite.toglobal(quadgrid, Ferrite.getneighborhood(topology, quadgrid, VertexIndex(5,4))) == [6,9,11,14] - - @test topology.face_skeleton == [FaceIndex(1,1),FaceIndex(1,2),FaceIndex(1,3),FaceIndex(1,4), + @test Set(Ferrite.getneighborhood(topology, quadgrid, VertexIndex(1,1))) == Set([VertexIndex(1,2), VertexIndex(1,4)]) + @test Set(Ferrite.getneighborhood(topology, quadgrid, VertexIndex(2,1))) == Set([VertexIndex(1,1), VertexIndex(1,3), VertexIndex(2,2), VertexIndex(2,4)]) + @test Set(Ferrite.getneighborhood(topology, quadgrid, VertexIndex(5,4))) == Set([VertexIndex(4,2), VertexIndex(4,4), VertexIndex(5,1), VertexIndex(5,3), VertexIndex(7,1), VertexIndex(7,3), VertexIndex(8,2), VertexIndex(8,4)]) + @test Set(Ferrite.getneighborhood(topology, quadgrid, VertexIndex(1,1),true)) == Set([VertexIndex(1,2), VertexIndex(1,4), VertexIndex(1,1)]) + @test Set(Ferrite.getneighborhood(topology, quadgrid, VertexIndex(2,1),true)) == Set([VertexIndex(1,1), VertexIndex(1,3), VertexIndex(2,2), VertexIndex(2,4), VertexIndex(1,2), VertexIndex(2,1)]) + @test Set(Ferrite.getneighborhood(topology, quadgrid, VertexIndex(5,4),true)) == Set([VertexIndex(4,2), VertexIndex(4,4), VertexIndex(5,1), VertexIndex(5,3), VertexIndex(7,1), VertexIndex(7,3), VertexIndex(8,2), VertexIndex(8,4), VertexIndex(4,3), VertexIndex(5,4), VertexIndex(7,2), VertexIndex(8,1)]) + @test Set(Ferrite.toglobal(quadgrid, Ferrite.getneighborhood(topology, quadgrid, VertexIndex(1,1)))) == Set([2,5]) + @test Set(Ferrite.toglobal(quadgrid, Ferrite.getneighborhood(topology, quadgrid, VertexIndex(2,1)))) == Set([1,6,3]) + @test Set(Ferrite.toglobal(quadgrid, Ferrite.getneighborhood(topology, quadgrid, VertexIndex(5,4)))) == Set([6,9,11,14]) + + @test Set(topology.face_skeleton) == Set([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), FaceIndex(4,2),FaceIndex(4,3),FaceIndex(4,4), FaceIndex(5,2),FaceIndex(5,3),FaceIndex(6,2),FaceIndex(6,3), FaceIndex(7,2),FaceIndex(7,3),FaceIndex(7,4), - FaceIndex(8,2),FaceIndex(8,3),FaceIndex(9,2),FaceIndex(9,3)] + FaceIndex(8,2),FaceIndex(8,3),FaceIndex(9,2),FaceIndex(9,3)]) @test length(topology.face_skeleton) == 4*3 + 3*4 quadratic_quadgrid = generate_grid(QuadraticQuadrilateral,(3,3)) @@ -378,6 +391,7 @@ end @test all(quadgrid_topology.vertex_neighbor .== topology.vertex_neighbor) quadratic_patches = Vector{Int}[Ferrite.getneighborhood(quadgrid_topology, quadratic_quadgrid, CellIndex(i)) for i in 1:getncells(quadratic_quadgrid)] @test all(patches .== quadratic_patches) + # # +-----+-----+-----+ # | 7 | 8 | 9 |