Skip to content

Commit

Permalink
Topology fixes and docs clean up, fixes #518, fixes #453. (#455)
Browse files Browse the repository at this point in the history
  • Loading branch information
@koehlerson
koehlerson authored Feb 20, 2023
1 parent 9ce3367 commit 88c619c
Show file tree
Hide file tree
Showing 3 changed files with 126 additions and 39 deletions.
1 change: 1 addition & 0 deletions docs/src/reference/grid.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ getcoordinates!
Ferrite.ExclusiveTopology
Ferrite.getneighborhood
Ferrite.faceskeleton
Ferrite.toglobal
```

### Grid Sets Utility
Expand Down
86 changes: 73 additions & 13 deletions src/Grid/grid.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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])
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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
"""
Expand Down
78 changes: 52 additions & 26 deletions test/test_grid_dofhandler_vtk.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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 |
Expand All @@ -297,27 +303,39 @@ 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 |
# +-----+-----+-----+
# | 4 | 5 | 6 |
# +-----+-----+-----+
# | 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)]
Expand All @@ -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),
Expand All @@ -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 |
# +-----+-----+-----+
Expand All @@ -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)
Expand All @@ -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
Expand Down

0 comments on commit 88c619c

Please sign in to comment.