Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

topology fixes and docs clean up #455

Merged
merged 9 commits into from
Feb 20, 2023
Merged

topology fixes and docs clean up #455

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
b4c362d
fix include self for vertices and docs clean up
Jun 6, 2022
bc0ccf6
Update src/Grid/grid.jl
koehlerson Jun 8, 2022
9e4cdf9
sizehinting
Jun 8, 2022
6a7528b
add edge self neighborhood
Jun 8, 2022
b4c1778
fix grid in test
Jun 8, 2022
5bf0b72
fix for higher order
koehlerson Oct 12, 2022
0f9bcd0
add access hotfix and regression tests
koehlerson Oct 13, 2022
896d1a4
replace high order fix with @zkk960317 suggestion and add more regres…
koehlerson Oct 18, 2022
7b42664
add getneighborhood call with FaceIndex to tests
koehlerson Feb 19, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/src/reference/grid.md
View file
Open in desktop
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
koehlerson marked this conversation as resolved.
Show resolved Hide resolved
```

### Grid Sets Utility
Expand Down
86 changes: 73 additions & 13 deletions src/Grid/grid.jl
View file
Open in desktop
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`).
termi-official marked this conversation as resolved.
Show resolved Hide resolved
"""
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
Open in desktop
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)]
koehlerson marked this conversation as resolved.
Show resolved Hide resolved
@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