diff --git a/docs/src/reference/grid.md b/docs/src/reference/grid.md index 4e834d9770..468e5856be 100644 --- a/docs/src/reference/grid.md +++ b/docs/src/reference/grid.md @@ -42,6 +42,7 @@ getcoordinates! Ferrite.ExclusiveTopology Ferrite.getneighborhood Ferrite.faceskeleton +Ferrite.toglobal ``` ### Grid Sets Utility diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index a5f94965e7..9edf6d0a91 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -32,6 +32,8 @@ struct Cell{dim,N,M} <: AbstractCell{dim,N,M} end nfaces(c::C) where {C<:AbstractCell} = nfaces(typeof(c)) nfaces(::Type{<:AbstractCell{dim,N,M}}) where {dim,N,M} = M +nedges(c::C) where {C<:AbstractCell} = length(edges(c)) +nvertices(c::C) where {C<:AbstractCell} = length(vertices(c)) nnodes(c::C) where {C<:AbstractCell} = nnodes(typeof(c)) nnodes(::Type{<:AbstractCell{dim,N,M}}) where {dim,N,M} = N @@ -193,8 +195,8 @@ function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell cell_neighbor_table[cellid] = EntityNeighborhood(CellIndex.(cell_neighbors)) for neighbor in cell_neighbors - neighbor_local_ids = findall(x->x in cell.nodes, cells[neighbor].nodes) - cell_local_ids = findall(x->x in cells[neighbor].nodes, cell.nodes) + 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]) # vertex neighbor if length(cell_local_ids) == 1 _vertex_neighbor!(V_vertex, I_vertex, J_vertex, cellid, cell, neighbor_local_ids, neighbor, cells[neighbor]) @@ -207,7 +209,43 @@ function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell end end end - + + celltype = eltype(cells) + if isconcretetype(celltype) + dim = getdim(cells[1]) + _nvertices = nvertices(cells[1]) + push!(V_vertex,zero(EntityNeighborhood{VertexIndex})) + push!(I_vertex,1); push!(J_vertex,_nvertices) + if dim > 1 + _nfaces = nfaces(cells[1]) + push!(V_face,zero(EntityNeighborhood{FaceIndex})) + push!(I_face,1); push!(J_face,_nfaces) + end + if dim > 2 + _nedges = nedges(cells[1]) + push!(V_edge,zero(EntityNeighborhood{EdgeIndex})) + push!(I_edge,1); push!(J_edge,_nedges) + end + else + celltypes = typeof.(cells) + for celltype in celltypes + celltypeidx = findfirst(x->typeof(x)==celltype,cells) + dim = getdim(cells[celltypeidx]) + _nvertices = nvertices(cells[celltypeidx]) + push!(V_vertex,zero(EntityNeighborhood{VertexIndex})) + push!(I_vertex,celltypeidx); push!(J_vertex,_nvertices) + if dim > 1 + _nfaces = nfaces(cells[celltypeidx]) + push!(V_face,zero(EntityNeighborhood{FaceIndex})) + push!(I_face,celltypeidx); push!(J_face,_nfaces) + end + if dim > 2 + _nedges = nedges(cells[celltypeidx]) + push!(V_edge,zero(EntityNeighborhood{EdgeIndex})) + push!(I_edge,celltypeidx); push!(J_edge,_nedges) + end + end + end face_neighbor = sparse(I_face,J_face,V_face) vertex_neighbor = sparse(I_vertex,J_vertex,V_vertex) edge_neighbor = sparse(I_edge,J_edge,V_edge) @@ -341,10 +379,10 @@ end # Grid utility functions # ########################## """ - getneighborhood(top::ExclusiveTopology, grid::Grid{dim,C,T}, cellidx::CellIndex, include_self=false) - getneighborhood(top::ExclusiveTopology, grid::Grid{dim,C,T}, faceidx::FaceIndex, include_self=false) - getneighborhood(top::ExclusiveTopology, grid::Grid{dim,C,T}, vertexidx::VertexIndex, include_self=false) - getneighborhood(top::ExclusiveTopology, grid::Grid{dim,C,T}, edgeidx::EdgeIndex, include_self=false) + getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, cellidx::CellIndex, include_self=false) + getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, faceidx::FaceIndex, include_self=false) + getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, vertexidx::VertexIndex, include_self=false) + getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, edgeidx::EdgeIndex, include_self=false) Returns all directly connected entities of the same type, i.e. calling the function with a `VertexIndex` will return a list of directly connected vertices (connected via face/edge). If `include_self` is true, the given `*Index` is included @@ -375,17 +413,34 @@ function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, vertexidx:: cell_vertices = vertices(getcells(grid,cellid)) global_vertexid = cell_vertices[local_vertexid] if include_self - return [top.vertex_vertex_neighbor[global_vertexid].neighbor_info; vertexidx] + vertex_to_cell = top.vertex_to_cell[global_vertexid] + self_reference_local = Vector{VertexIndex}(undef,length(vertex_to_cell)) + for (i,cellid) in enumerate(vertex_to_cell) + local_vertex = VertexIndex(cellid,findfirst(x->x==global_vertexid,vertices(getcells(grid,cellid)))) + self_reference_local[i] = local_vertex + end + return [top.vertex_vertex_neighbor[global_vertexid].neighbor_info; self_reference_local] else return top.vertex_vertex_neighbor[global_vertexid].neighbor_info end end function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid{3}, edgeidx::EdgeIndex, include_self=false) - if include_self - return [top.edge_neighbor[edgeidx[1],edgeidx[2]].neighbor_info; edgeidx] + cellid, local_edgeidx = edgeidx[1], edgeidx[2] + cell_edges = edges(getcells(grid,cellid)) + nonlocal_edgeid = cell_edges[local_edgeidx] + cell_neighbors = getneighborhood(top,grid,CellIndex(cellid)) + self_reference_local = EdgeIndex[] + for cellid in cell_neighbors + local_neighbor_edgeid = findfirst(x->issubset(x,nonlocal_edgeid),edges(getcells(grid,cellid))) + local_neighbor_edgeid === nothing && continue + local_edge = EdgeIndex(cellid,local_neighbor_edgeid) + push!(self_reference_local, local_edge) + end + if include_self + return unique([top.edge_neighbor[cellid, local_edgeidx].neighbor_info; self_reference_local; edgeidx]) else - return top.edge_neighbor[edgeidx[1],edgeidx[2]].neighbor_info + return unique([top.edge_neighbor[cellid, local_edgeidx].neighbor_info; self_reference_local]) end end @@ -396,8 +451,13 @@ Returns an iterateable face skeleton. The skeleton consists of `FaceIndex` that """ faceskeleton(top::ExclusiveTopology, grid::AbstractGrid) = top.face_skeleton -toglobal(grid::Grid,vertexidx::VertexIndex) = vertices(getcells(grid,vertexidx[1]))[vertexidx[2]] -toglobal(grid::Grid,vertexidx::Vector{VertexIndex}) = unique(toglobal.((grid,),vertexidx)) +""" + toglobal(grid::AbstractGrid, vertexidx::VertexIndex) -> Int + toglobal(grid::AbstractGrid, vertexidx::Vector{VertexIndex}) -> Vector{Int} +This function takes the local vertex representation (a `VertexIndex`) and looks up the unique global id (an `Int`). +""" +toglobal(grid::AbstractGrid,vertexidx::VertexIndex) = vertices(getcells(grid,vertexidx[1]))[vertexidx[2]] +toglobal(grid::AbstractGrid,vertexidx::Vector{VertexIndex}) = unique(toglobal.((grid,),vertexidx)) @inline getdim(::AbstractGrid{dim}) where {dim} = dim """ diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 97d42f6645..46851bb43f 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/test/test_grid_dofhandler_vtk.jl @@ -235,7 +235,7 @@ end @test topology.vertex_neighbor[4,4] == Ferrite.EntityNeighborhood(VertexIndex(5,2)) @test topology.vertex_neighbor[5,2] == Ferrite.EntityNeighborhood(VertexIndex(4,4)) @test topology.vertex_neighbor[6,1] == Ferrite.EntityNeighborhood(VertexIndex(3,3)) - #test face neighbor maps cellid and local face id to neighbor id and neighbor local face id + #test face neighbor maps cellid and local face id to neighbor id and neighbor local face id @test topology.face_neighbor[1,2] == Ferrite.EntityNeighborhood(FaceIndex(2,4)) @test topology.face_neighbor[1,3] == Ferrite.EntityNeighborhood(FaceIndex(3,1)) @test topology.face_neighbor[2,3] == Ferrite.EntityNeighborhood(FaceIndex(4,1)) @@ -268,21 +268,27 @@ end # | 1 | 2 | # (10) +-----+-----+(12) # (11) - hexgrid = generate_grid(Hexahedron,(2,2,1)) + 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 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 topology.edge_neighbor[3,10] == Ferrite.EntityNeighborhood(EdgeIndex(2,12)) @test topology.edge_neighbor[4,9] == Ferrite.EntityNeighborhood(EdgeIndex(1,11)) - @test isempty(topology.vertex_neighbor) - @test topology.face_neighbor[1,3] == Ferrite.EntityNeighborhood(FaceIndex(2,5)) - @test topology.face_neighbor[1,4] == Ferrite.EntityNeighborhood(FaceIndex(3,2)) - @test topology.face_neighbor[2,4] == Ferrite.EntityNeighborhood(FaceIndex(4,2)) - @test topology.face_neighbor[2,5] == Ferrite.EntityNeighborhood(FaceIndex(1,3)) - @test topology.face_neighbor[3,2] == Ferrite.EntityNeighborhood(FaceIndex(1,4)) - @test topology.face_neighbor[3,3] == Ferrite.EntityNeighborhood(FaceIndex(4,5)) - @test topology.face_neighbor[4,2] == Ferrite.EntityNeighborhood(FaceIndex(2,4)) - @test topology.face_neighbor[4,5] == Ferrite.EntityNeighborhood(FaceIndex(3,3)) + @test getneighborhood(topology,hexgrid,FaceIndex((1,3))) == [FaceIndex((2,5))] + @test getneighborhood(topology,hexgrid,FaceIndex((1,4))) == [FaceIndex((3,2))] + @test getneighborhood(topology,hexgrid,FaceIndex((2,4))) == [FaceIndex((4,2))] + @test getneighborhood(topology,hexgrid,FaceIndex((2,5))) == [FaceIndex((1,3))] + @test getneighborhood(topology,hexgrid,FaceIndex((3,2))) == [FaceIndex((1,4))] + @test getneighborhood(topology,hexgrid,FaceIndex((3,3))) == [FaceIndex((4,5))] + @test getneighborhood(topology,hexgrid,FaceIndex((4,2))) == [FaceIndex((2,4))] + @test getneighborhood(topology,hexgrid,FaceIndex((4,5))) == [FaceIndex((3,3))] + serendipitygrid = generate_grid(Cell{3,20,6},(2,2,1)) + stopology = ExclusiveTopology(serendipitygrid) + # regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/518 + @test all(stopology.face_neighbor .== topology.face_neighbor) + @test all(stopology.vertex_neighbor .== topology.vertex_neighbor) # +-----+-----+ # |\ 6 |\ 8 | @@ -297,19 +303,31 @@ end trigrid = generate_grid(Triangle,(2,2)) topology = ExclusiveTopology(trigrid) @test topology.vertex_neighbor[3,3] == Ferrite.EntityNeighborhood([VertexIndex(5,2),VertexIndex(6,1),VertexIndex(7,1)]) + quadtrigrid = generate_grid(QuadraticTriangle,(2,2)) + quadtopology = ExclusiveTopology(trigrid) + # add more regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/518 + @test all(quadtopology.face_neighbor .== topology.face_neighbor) + @test all(quadtopology.vertex_neighbor .== topology.vertex_neighbor) # test mixed grid cells = [ Hexahedron((1, 2, 3, 4, 5, 6, 7, 8)), - Quadrilateral((3, 2, 9, 10)), + 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{2,Float64}, 10)] + nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 18)] grid = Grid(cells, nodes) topology = ExclusiveTopology(grid) - @test isempty(topology.vertex_neighbor) - @test topology.face_neighbor[2,1] == Ferrite.EntityNeighborhood(EdgeIndex(1,2)) - @test topology.edge_neighbor[1,2] == Ferrite.EntityNeighborhood(FaceIndex(2,1)) -# + + @test topology.face_neighbor[3,4] == Ferrite.EntityNeighborhood(EdgeIndex(1,2)) + @test topology.edge_neighbor[1,2] == Ferrite.EntityNeighborhood(FaceIndex(3,4)) + # regression that it doesn't error for boundary faces, see https://github.com/Ferrite-FEM/Ferrite.jl/issues/518 + @test topology.face_neighbor[1,6] == topology.face_neighbor[1,1] == zero(Ferrite.EntityNeighborhood{FaceIndex}) + @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 | # +-----+-----+-----+ @@ -317,7 +335,7 @@ end # +-----+-----+-----+ # | 1 | 2 | 3 | # +-----+-----+-----+ -# test application: form level 1 neighborhood patches of elements +# test application: form level 1 neighborhood patches of elements quadgrid = generate_grid(Quadrilateral,(3,3)) topology = ExclusiveTopology(quadgrid) patches = Vector{Int}[Ferrite.getneighborhood(topology, quadgrid, CellIndex(i)) for i in 1:getncells(quadgrid)] @@ -331,14 +349,17 @@ end @test issubset([4,5,8], patches[7]) @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), FaceIndex(2,1),FaceIndex(2,2),FaceIndex(2,3), FaceIndex(3,1),FaceIndex(3,2),FaceIndex(3,3), @@ -349,9 +370,14 @@ end @test length(topology.face_skeleton) == 4*3 + 3*4 quadratic_quadgrid = generate_grid(QuadraticQuadrilateral,(3,3)) - quadgrid_topology = ExclusiveTopology(grid) - quadgrid_topology.face_skeleton == topology.face_skeleton -# + quadgrid_topology = ExclusiveTopology(quadratic_quadgrid) + @test quadgrid_topology.face_skeleton == topology.face_skeleton + # add more regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/518 + @test all(quadgrid_topology.face_neighbor .== topology.face_neighbor) + @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 | # +-----+-----+-----+ @@ -365,7 +391,7 @@ end Ferrite.reinit!(fv, coords, faceid) end reinit!(fv::FaceValues, faceid::FaceIndex, grid) = reinit!(fv,faceid[1],faceid[2],grid) # wrapper for reinit!(fv,cellid,faceid,grid) - face_neighbors_ele5 = nonzeros(topology.face_neighbor[5,:]) + face_neighbors_ele5 = nonzeros(topology.face_neighbor[5,:]) ip = Lagrange{2, RefCube, 1}() qr_face = QuadratureRule{1, RefCube}(2) fv_ele = FaceVectorValues(qr_face, ip) @@ -382,11 +408,11 @@ end u_5_n = function_value(fv_ele, q_point, u_ele5) ⋅ normal_5 for neighbor_entity in face_neighbors_ele5[ele_faceid].neighbor_info # only one entity can be changed to get rid of the for loop reinit!(fv_neighbor, neighbor_entity, quadgrid) - normal_neighbor = getnormal(fv_neighbor, q_point) + normal_neighbor = getnormal(fv_neighbor, q_point) u_neighbor = function_value(fv_neighbor, q_point, u_neighbors) ⋅ normal_neighbor jump_int += (u_5_n + u_neighbor) * dΩ jump_abs += abs(u_5_n + u_neighbor) * dΩ - end + end end end @test isapprox(jump_abs, 2/3*2*4,atol=1e-6) # 2*4*0.66666, jump is always 2, 4 sides, length =0.66