From 271091e36ef707394f6d0d8ad05215ac049d118a Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 3 Apr 2023 16:52:38 +0200 Subject: [PATCH 1/9] Update topology construction to allows Wedges. --- src/Grid/grid.jl | 91 +++++++++++++++++++------------- test/test_grid_dofhandler_vtk.jl | 36 +++++++------ 2 files changed, 74 insertions(+), 53 deletions(-) diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index ee862abab1..bac1c94728 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -135,19 +135,15 @@ function Base.show(io::IO, ::MIME"text/plain", n::EntityNeighborhood) end """ - face_npoints(::AbstractCell{dim,N,M) -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(::Cell{2,N,M}) where {N,M} = 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{dim,N,M) -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(::Cell{3,N,6}) where N = 4 -face_npoints(::Cell{3,N,4}) where N = 3 -edge_npoints(::Cell{3,N,M}) where {N,M} = 2 +nvertices_on_edge(cell::AbstractCell, local_edge_index::Int) = length(edges(cell)[local_edge_index]) getdim(::Cell{dim}) where dim = dim @@ -156,8 +152,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. @@ -172,7 +167,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) @@ -188,15 +183,15 @@ struct ExclusiveTopology <: AbstractTopology 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) + cell_vertex_table = vertices.(cells) #needs generic interface for <: AbstractCell + vertex_cell_table = Dict{Int,Set{Int}}() + + for (cellid, cell_vertices) in enumerate(cell_vertex_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 @@ -206,29 +201,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 = 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))) + + 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 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]) + for neighbor_cellid in cell_neighbors + cell_local_ids = findall(x->x in cell_vertex_table[neighbor_cellid], cell_vertex_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_vertex_table[cellid], cell_vertex_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_vertex_table[cellid], cell_vertex_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_vertex_table[cellid], cell_vertex_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]) diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 67c7e12d1c..5af7789044 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/test/test_grid_dofhandler_vtk.jl @@ -68,7 +68,7 @@ end add!(ch, dbc) end close!(ch) - update!(ch, 0.0) + Ferrite.update!(ch, 0.0) rng = StableRNG(1234) u = rand(rng, ndofs(dofhandler)) apply!(u, ch) @@ -192,7 +192,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 @@ -271,9 +270,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))] @@ -284,6 +283,8 @@ end @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 @@ -327,6 +328,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 | @@ -350,23 +353,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)) @@ -377,6 +380,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 | From 3184424b92a48dd8d72567b47de1b4555fca9555 Mon Sep 17 00:00:00 2001 From: termi-official Date: Tue, 4 Apr 2023 22:39:31 +0200 Subject: [PATCH 2/9] Revert accidental test hotfix. --- test/test_grid_dofhandler_vtk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 5af7789044..6c224b9b37 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/test/test_grid_dofhandler_vtk.jl @@ -68,7 +68,7 @@ end add!(ch, dbc) end close!(ch) - Ferrite.update!(ch, 0.0) + update!(ch, 0.0) rng = StableRNG(1234) u = rand(rng, ndofs(dofhandler)) apply!(u, ch) From dd05d2c7ff65bdf766e31c6dc20ab476996d4ab3 Mon Sep 17 00:00:00 2001 From: termi-official Date: Thu, 13 Apr 2023 14:41:53 +0200 Subject: [PATCH 3/9] Remove splatting. --- src/Grid/grid.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index bac1c94728..778adf1350 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -202,7 +202,7 @@ function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell cell_neighbor_table = Vector{EntityNeighborhood{CellIndex}}(undef, length(cells)) for (cellid, cell) in enumerate(cells) - cell_neighbors = union([Set{Int}(vertex_cell_table[vertex]) for vertex ∈ vertices(cell) if vertex_cell_table[vertex] != cellid]...) + 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))) face_neighbors = Set{Int}() From a66023995f6871b29d0f440389e735ede7b4e4b0 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Fri, 5 May 2023 12:35:24 +0200 Subject: [PATCH 4/9] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Maximilian Köhler --- src/Grid/grid.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index 778adf1350..65a40e3ae0 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -207,7 +207,7 @@ function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell 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]) + 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 @@ -217,7 +217,7 @@ function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell 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]) + 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 From 0a92bb95714becc46f322a617eb9e7f36160d27e Mon Sep 17 00:00:00 2001 From: termi-official Date: Fri, 5 May 2023 12:40:04 +0200 Subject: [PATCH 5/9] Revert variable rename. --- src/Grid/grid.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index 65a40e3ae0..c5bfc833ee 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -183,10 +183,10 @@ struct ExclusiveTopology <: AbstractTopology end function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell - cell_vertex_table = vertices.(cells) #needs generic interface for <: AbstractCell + cell_vertices_table = vertices.(cells) #needs generic interface for <: AbstractCell vertex_cell_table = Dict{Int,Set{Int}}() - for (cellid, cell_vertices) in enumerate(cell_vertex_table) + 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) @@ -229,18 +229,18 @@ function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell end for neighbor_cellid in cell_neighbors - cell_local_ids = findall(x->x in cell_vertex_table[neighbor_cellid], cell_vertex_table[cellid]) + cell_local_ids = findall(x->x in cell_vertices_table[neighbor_cellid], cell_vertices_table[cellid]) # vertex neighbor if length(cell_local_ids) == 1 - neighbor_local_ids = findall(x->x in cell_vertex_table[cellid], cell_vertex_table[neighbor_cellid]) + 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 neighbor_cellid ∈ face_neighbors - neighbor_local_ids = findall(x->x in cell_vertex_table[cellid], cell_vertex_table[neighbor_cellid]) + 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 && neighbor_cellid ∈ edge_neighbors - neighbor_local_ids = findall(x->x in cell_vertex_table[cellid], cell_vertex_table[neighbor_cellid]) + 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 From 9551a86880bd3d557ec3f9fc9462b451a6519c5a Mon Sep 17 00:00:00 2001 From: termi-official Date: Fri, 5 May 2023 13:53:07 +0200 Subject: [PATCH 6/9] Fix 1D topology construction. --- src/Grid/grid.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index 0e739a75ca..2c90983ce9 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 @@ -419,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 ? [edges(cells[cell])...] : [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) @@ -435,6 +437,7 @@ function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell face_skeleton_local = Vector{FaceIndex}() fs_length = length(face_skeleton_global) for (cellid,cell) in enumerate(cells) + getdim(cell) == 1 && continue for (local_face_id,face) in enumerate(faces(cell)) push!(face_skeleton_global, first(sortface(face))) fs_length_new = length(face_skeleton_global) From 83a6449f4eb165cc9e76ef54cffaee319b1ba7f1 Mon Sep 17 00:00:00 2001 From: termi-official Date: Fri, 5 May 2023 14:10:07 +0200 Subject: [PATCH 7/9] Fix 1D topology face skeleton construction. --- src/Dofs/DofHandler.jl | 4 +++- src/Grid/grid.jl | 1 - 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index abc3d25806..844308a0bc 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 2c90983ce9..5900fb2970 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -437,7 +437,6 @@ function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell face_skeleton_local = Vector{FaceIndex}() fs_length = length(face_skeleton_global) for (cellid,cell) in enumerate(cells) - getdim(cell) == 1 && continue for (local_face_id,face) in enumerate(faces(cell)) push!(face_skeleton_global, first(sortface(face))) fs_length_new = length(face_skeleton_global) From 4da91510d1b41b72d487f37b7dc78a36dd1b6222 Mon Sep 17 00:00:00 2001 From: termi-official Date: Fri, 5 May 2023 14:35:32 +0200 Subject: [PATCH 8/9] Remove splatting. --- src/Grid/grid.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index 5900fb2970..d9533f973b 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -419,7 +419,7 @@ 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 ? [edges(cells[cell])...] : [faces(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)] other_vertices = findfirst.(x->x!=global_vertexid,neighbor_connected_faces) any(other_vertices .=== nothing) && continue From 08cd31fc6858dcd774d6a4e2fcb743bf1b3ccfa9 Mon Sep 17 00:00:00 2001 From: termi-official Date: Fri, 5 May 2023 15:36:15 +0200 Subject: [PATCH 9/9] Add regression test for 1D topology --- test/test_grid_dofhandler_vtk.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 79abe65a89..68fa4299d0 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/test/test_grid_dofhandler_vtk.jl @@ -214,6 +214,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 |