From 64eab2d13c930bf5ec57cbb6e6d401c5a5f7b959 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Thu, 7 Nov 2024 21:04:54 +0100 Subject: [PATCH 01/28] homeoffice push --- src/Ferrite.jl | 4 +++ src/iterators.jl | 63 +++++++++++++++++++++++++++++------------------- 2 files changed, 42 insertions(+), 25 deletions(-) diff --git a/src/Ferrite.jl b/src/Ferrite.jl index 7b358463be..0b72310434 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -107,6 +107,10 @@ struct FacetIndex <: BoundaryIndex idx::Tuple{Int, Int} # cell and side end +struct InterfaceIndex <: BoundaryIndex + idx::Tuple{Int, Int, Int, Int} +end + const AbstractVecOrSet{T} = Union{AbstractSet{T}, AbstractVector{T}} const IntegerCollection = AbstractVecOrSet{<:Integer} diff --git a/src/iterators.jl b/src/iterators.jl index 0abf2972b1..8dc54eafef 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -176,9 +176,9 @@ function InterfaceCache(gridordh::Union{AbstractGrid, AbstractDofHandler}) return InterfaceCache(fc_a, fc_b, Int[]) end -function reinit!(cache::InterfaceCache, facet_a::BoundaryIndex, facet_b::BoundaryIndex) - reinit!(cache.a, facet_a) - reinit!(cache.b, facet_b) +function reinit!(cache::InterfaceCache, interface::InterfaceIndex) + reinit!(cache.a, FacetIndex(interface.idx[1], interface.idx[2])) + reinit!(cache.b, FacetIndex(interface.idx[3], interface.idx[4])) resize!(cache.dofs, length(celldofs(cache.a)) + length(celldofs(cache.b))) for (i, d) in pairs(cache.a.dofs) cache.dofs[i] = d @@ -341,39 +341,52 @@ end `InterfaceIterator` is stateful and should not be used for things other than `for`-looping (e.g. broadcasting over, or collecting the iterator may yield unexpected results). """ -struct InterfaceIterator{IC <: InterfaceCache, G <: Grid} +struct InterfaceIterator{IC <: InterfaceCache} cache::IC - grid::G - topology::ExclusiveTopology + set::OrderedSet{InterfaceIndex} end +@inline _getcache(ii::InterfaceIterator) = ii.cache +@inline _getset(ii::InterfaceIterator) = ii.set + function InterfaceIterator( gridordh::Union{Grid, AbstractDofHandler}, + set_here::AbstractVecOrSet{FacetIndex}, topology::ExclusiveTopology = ExclusiveTopology(gridordh isa Grid ? gridordh : get_grid(gridordh)) ) grid = gridordh isa Grid ? gridordh : get_grid(gridordh) - return InterfaceIterator(InterfaceCache(gridordh), grid, topology) -end - -# Iterator interface -function Base.iterate(ii::InterfaceIterator{<:Any, <:Grid{sdim}}, state...) where {sdim} - neighborhood = get_facet_facet_neighborhood(ii.topology, ii.grid) # TODO: This could be moved to InterfaceIterator constructor (potentially type-instable for non-union or mixed grids) - while true - it = iterate(facetskeleton(ii.topology, ii.grid), state...) - it === nothing && return nothing - facet_a, state = it - if isempty(neighborhood[facet_a[1], facet_a[2]]) - continue - end - neighbors = neighborhood[facet_a[1], facet_a[2]] - length(neighbors) > 1 && error("multiple neighboring faces not supported yet") - facet_b = neighbors[1] - reinit!(ii.cache, facet_a, facet_b) - return (ii.cache, state) + if gridordh isa DofHandler + # Keep here to maintain same settings as for CellIterator + _check_same_celltype(grid, set_here) end - return + neighborhood = get_facet_facet_neighborhood(topology, grid) + @assert all(facet -> !isempty(neighborhood[facet[1], facet[2]]), set) "Facets in the set must have neighbors" + set = OrderedSet{InterfaceIndex}() + sizehint!(set, ninterfaces) + for facet in set_here + neighbor = neighborhood[facet[1], facet[2]][] + push!(set, InterfaceIndex(facet[1], facet[2], neighbor[1], neighbor[2])) + end + return InterfaceIterator(InterfaceCache(gridordh), set) end +function InterfaceIterator( + gridordh::Union{Grid, AbstractDofHandler}, + topology::ExclusiveTopology = ExclusiveTopology(gridordh isa Grid ? gridordh : get_grid(gridordh)) + ) + grid = gridordh isa Grid ? gridordh : get_grid(gridordh) + if gridordh isa DofHandler + # Keep here to maintain same settings as for CellIterator + _check_same_celltype(grid, 1:getncells(grid)) + end + neighborhood = get_facet_facet_neighborhood(topology, grid) + fs = facetskeleton(topology, grid) + ninterfaces = count(facet -> !isempty(neighborhood[facet[1], facet[2]]), fs) + set = OrderedSet{InterfaceIndex}() + sizehint!(set, ninterfaces) + map!(facet -> InterfaceIndex((facet[1], facet[2], neighborhood[facet[1], facet[2]][1][1], neighborhood[facet[1], facet[2]][1][2])), set, filter(facet -> !isempty(neighborhood[facet[1], facet[2]]), fs)) + return InterfaceIterator(InterfaceCache(gridordh), set) +end # Iterator interface for CellIterator/FacetIterator const GridIterators{C} = Union{CellIterator{C}, FacetIterator{C}, InterfaceIterator{C}} From 050f22622516488a73f6040633dc4fe287244d05 Mon Sep 17 00:00:00 2001 From: Abdulaziz Hamid Mohamed Date: Thu, 7 Nov 2024 22:44:20 +0100 Subject: [PATCH 02/28] oopsie --- src/iterators.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/iterators.jl b/src/iterators.jl index 8dc54eafef..d722ab6357 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -384,7 +384,10 @@ function InterfaceIterator( ninterfaces = count(facet -> !isempty(neighborhood[facet[1], facet[2]]), fs) set = OrderedSet{InterfaceIndex}() sizehint!(set, ninterfaces) - map!(facet -> InterfaceIndex((facet[1], facet[2], neighborhood[facet[1], facet[2]][1][1], neighborhood[facet[1], facet[2]][1][2])), set, filter(facet -> !isempty(neighborhood[facet[1], facet[2]]), fs)) + for facet in fs + isempty(neighborhood[facet[1], facet[2]]) && continue + push!(set, InterfaceIndex((facet[1], facet[2], neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2]))) + end return InterfaceIterator(InterfaceCache(gridordh), set) end From ec96036f03441f5086754503e50d98ee5ec9a89c Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 8 Nov 2024 17:05:45 +0100 Subject: [PATCH 03/28] use Set --- src/iterators.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index d722ab6357..b2064190a2 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -341,9 +341,9 @@ end `InterfaceIterator` is stateful and should not be used for things other than `for`-looping (e.g. broadcasting over, or collecting the iterator may yield unexpected results). """ -struct InterfaceIterator{IC <: InterfaceCache} +struct InterfaceIterator{IC <: InterfaceCache, SetType <: AbstractSet{InterfaceIndex}} cache::IC - set::OrderedSet{InterfaceIndex} + set::SetType end @inline _getcache(ii::InterfaceIterator) = ii.cache @@ -360,12 +360,12 @@ function InterfaceIterator( _check_same_celltype(grid, set_here) end neighborhood = get_facet_facet_neighborhood(topology, grid) - @assert all(facet -> !isempty(neighborhood[facet[1], facet[2]]), set) "Facets in the set must have neighbors" - set = OrderedSet{InterfaceIndex}() - sizehint!(set, ninterfaces) + @assert all(facet -> !isempty(neighborhood[facet[1], facet[2]]), set_here) "Facets in the set must have neighbors" + set = Set{InterfaceIndex}() + sizehint!(set, length(set_here)) for facet in set_here neighbor = neighborhood[facet[1], facet[2]][] - push!(set, InterfaceIndex(facet[1], facet[2], neighbor[1], neighbor[2])) + push!(set, InterfaceIndex((facet[1], facet[2], neighbor[1], neighbor[2]))) end return InterfaceIterator(InterfaceCache(gridordh), set) end @@ -382,7 +382,7 @@ function InterfaceIterator( neighborhood = get_facet_facet_neighborhood(topology, grid) fs = facetskeleton(topology, grid) ninterfaces = count(facet -> !isempty(neighborhood[facet[1], facet[2]]), fs) - set = OrderedSet{InterfaceIndex}() + set = Set{InterfaceIndex}() sizehint!(set, ninterfaces) for facet in fs isempty(neighborhood[facet[1], facet[2]]) && continue From 6f4f04027bf341851a145c8da31e63e4955bceee Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 8 Nov 2024 17:45:13 +0100 Subject: [PATCH 04/28] now tests shall pass --- src/Ferrite.jl | 5 ++++- src/Grid/grid.jl | 8 ++++---- src/exports.jl | 1 + src/iterators.jl | 4 ++-- test/test_grid_dofhandler_vtk.jl | 6 +++--- 5 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/Ferrite.jl b/src/Ferrite.jl index 0b72310434..3b220ac5c1 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -107,8 +107,11 @@ struct FacetIndex <: BoundaryIndex idx::Tuple{Int, Int} # cell and side end +""" +An `InterfaceIndex` wraps an (Int, Int, Int, Int) and defines an interface by pointing to a (cell_here, facet_here, cell_there, facet_there). +""" struct InterfaceIndex <: BoundaryIndex - idx::Tuple{Int, Int, Int, Int} + idx::Tuple{Int, Int, Int, Int} # cell - side - cell - side end const AbstractVecOrSet{T} = Union{AbstractSet{T}, AbstractVector{T}} diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index 3d27431bd1..2c67bf3101 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -641,10 +641,10 @@ boundaryfunction(::Type{EdgeIndex}) = edges boundaryfunction(::Type{VertexIndex}) = vertices boundaryfunction(::Type{FacetIndex}) = facets -for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex, :FacetIndex) +for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex, :FacetIndex, :InterfaceIndex) @eval begin #Constructor - ($INDEX)(a::Int, b::Int) = ($INDEX)((a, b)) + ($INDEX)(args...) = ($INDEX)((args...,)) Base.getindex(I::($INDEX), i::Int) = I.idx[i] @@ -653,8 +653,8 @@ for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex, :FacetIndex) # Necessary to check if, e.g. `(cellid, faceidx) in faceset` Base.isequal(x::$INDEX, y::$INDEX) = x.idx == y.idx - Base.isequal(x::Tuple{Int, Int}, y::$INDEX) = x[1] == y.idx[1] && x[2] == y.idx[2] - Base.isequal(y::$INDEX, x::Tuple{Int, Int}) = x[1] == y.idx[1] && x[2] == y.idx[2] + Base.isequal(x::NTuple{N, Int}, y::$INDEX) where N = all(i -> x[i] == y.idx[i], 1:N) + Base.isequal(y::$INDEX, x::NTuple{N, Int}) where N = all(i -> x[i] == y.idx[i], 1:N) Base.hash(x::$INDEX, h::UInt) = hash(x.idx, h) end end diff --git a/src/exports.jl b/src/exports.jl index c710861c15..04afb094aa 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -77,6 +77,7 @@ export EdgeIndex, VertexIndex, FacetIndex, + InterfaceIndex, geometric_interpolation, ExclusiveTopology, getneighborhood, diff --git a/src/iterators.jl b/src/iterators.jl index b2064190a2..5c9488e04d 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -365,7 +365,7 @@ function InterfaceIterator( sizehint!(set, length(set_here)) for facet in set_here neighbor = neighborhood[facet[1], facet[2]][] - push!(set, InterfaceIndex((facet[1], facet[2], neighbor[1], neighbor[2]))) + push!(set, InterfaceIndex(facet[1], facet[2], neighbor[1], neighbor[2])) end return InterfaceIterator(InterfaceCache(gridordh), set) end @@ -386,7 +386,7 @@ function InterfaceIterator( sizehint!(set, ninterfaces) for facet in fs isempty(neighborhood[facet[1], facet[2]]) && continue - push!(set, InterfaceIndex((facet[1], facet[2], neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2]))) + push!(set, InterfaceIndex(facet[1], facet[2], neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2])) end return InterfaceIterator(InterfaceCache(gridordh), set) end diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 9f1ad073cd..0ac5da3faf 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/test/test_grid_dofhandler_vtk.jl @@ -232,12 +232,12 @@ end # InterfaceCache grid = generate_grid(Quadrilateral, (2, 1)) ic = InterfaceCache(grid) - reinit!(ic, FaceIndex(1, 2), FaceIndex(2, 4)) + reinit!(ic, InterfaceIndex(1, 2, 2, 4)) @test interfacedofs(ic) == Int[] # Empty because no DofHandler given ip = DiscontinuousLagrange{RefQuadrilateral, 1}() dh = DofHandler(grid); add!(dh, :u, ip); close!(dh) ic = InterfaceCache(dh) - reinit!(ic, FaceIndex(1, 2), FaceIndex(2, 4)) + reinit!(ic, InterfaceIndex(1, 2, 2, 4)) @test interfacedofs(ic) == collect(1:8) # Mixed Elements dim = 2 @@ -254,7 +254,7 @@ end sdh2 = SubDofHandler(dh, Set([2])); add!(sdh2, :u, ip2) close!(dh) ic = InterfaceCache(dh) - reinit!(ic, FaceIndex(1, 2), FaceIndex(2, 3)) + reinit!(ic, InterfaceIndex(1, 2, 2, 3)) @test interfacedofs(ic) == collect(1:7) # Unit test of some utilities mixed_grid = Grid( From 932b3664974dd85bf74f592a6d4b43e9c7a975c1 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 8 Nov 2024 18:51:00 +0100 Subject: [PATCH 05/28] runic --- 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 2c67bf3101..29a265928a 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -653,8 +653,8 @@ for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex, :FacetIndex, :InterfaceIndex # Necessary to check if, e.g. `(cellid, faceidx) in faceset` Base.isequal(x::$INDEX, y::$INDEX) = x.idx == y.idx - Base.isequal(x::NTuple{N, Int}, y::$INDEX) where N = all(i -> x[i] == y.idx[i], 1:N) - Base.isequal(y::$INDEX, x::NTuple{N, Int}) where N = all(i -> x[i] == y.idx[i], 1:N) + Base.isequal(x::NTuple{N, Int}, y::$INDEX) where {N} = all(i -> x[i] == y.idx[i], 1:N) + Base.isequal(y::$INDEX, x::NTuple{N, Int}) where {N} = all(i -> x[i] == y.idx[i], 1:N) Base.hash(x::$INDEX, h::UInt) = hash(x.idx, h) end end From b14251e3d4e1ce8bfae406ef544d7e8fa11d759f Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 8 Nov 2024 19:33:12 +0100 Subject: [PATCH 06/28] Init iterator with subdofhandlers --- src/iterators.jl | 38 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 36 insertions(+), 2 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index 5c9488e04d..d72e616c14 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -176,6 +176,12 @@ function InterfaceCache(gridordh::Union{AbstractGrid, AbstractDofHandler}) return InterfaceCache(fc_a, fc_b, Int[]) end +function InterfaceCache(sdh_here::SubDofHandler, sdh_there::SubDofHandler) + fc_a = FacetCache(sdh_here) + fc_b = FacetCache(sdh_there) + return InterfaceCache(fc_a, fc_b, Int[]) +end + function reinit!(cache::InterfaceCache, interface::InterfaceIndex) reinit!(cache.a, FacetIndex(interface.idx[1], interface.idx[2])) reinit!(cache.b, FacetIndex(interface.idx[3], interface.idx[4])) @@ -350,7 +356,7 @@ end @inline _getset(ii::InterfaceIterator) = ii.set function InterfaceIterator( - gridordh::Union{Grid, AbstractDofHandler}, + gridordh::Union{Grid, DofHandler}, set_here::AbstractVecOrSet{FacetIndex}, topology::ExclusiveTopology = ExclusiveTopology(gridordh isa Grid ? gridordh : get_grid(gridordh)) ) @@ -371,7 +377,7 @@ function InterfaceIterator( end function InterfaceIterator( - gridordh::Union{Grid, AbstractDofHandler}, + gridordh::Union{Grid, DofHandler}, topology::ExclusiveTopology = ExclusiveTopology(gridordh isa Grid ? gridordh : get_grid(gridordh)) ) grid = gridordh isa Grid ? gridordh : get_grid(gridordh) @@ -391,6 +397,34 @@ function InterfaceIterator( return InterfaceIterator(InterfaceCache(gridordh), set) end +function InterfaceIterator( + sdh_here::SubDofHandler, + sdh_there::SubDofHandler, + topology::ExclusiveTopology = ExclusiveTopology(get_grid(sdh_here.dh)) + ) + grid = get_grid(sdh_here.dh) + neighborhood = get_facet_facet_neighborhood(topology, grid) + fs = facetskeleton(topology, grid) + ninterfaces = 0 + for facet in fs + facet[1] ∈ sdh_here.cellset || continue + neighbors = neighborhood[facet[1], facet[2]] + isempty(neighbors) && continue + neighbors[][1] ∈ sdh_there.cellset || continue + ninterfaces += 1 + end + set = Set{InterfaceIndex}() + sizehint!(set, ninterfaces) + for facet in fs + facet[1] ∈ sdh_here.cellset || continue + neighbors = neighborhood[facet[1], facet[2]] + isempty(neighbors) && continue + neighbors[][1] ∈ sdh_there.cellset || continue + push!(set, InterfaceIndex(facet[1], facet[2], neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2])) + end + return InterfaceIterator(InterfaceCache(sdh_here, sdh_there), set) +end + # Iterator interface for CellIterator/FacetIterator const GridIterators{C} = Union{CellIterator{C}, FacetIterator{C}, InterfaceIterator{C}} From 4659dba1959257dbcdb8a4fe77420b97701c5c35 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 11 Nov 2024 16:47:42 +0100 Subject: [PATCH 07/28] oopsie --- 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 29a265928a..de38b68992 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -649,7 +649,7 @@ for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex, :FacetIndex, :InterfaceIndex Base.getindex(I::($INDEX), i::Int) = I.idx[i] #To be able to do a,b = faceidx - Base.iterate(I::($INDEX), state::Int = 1) = (state == 3) ? nothing : (I[state], state + 1) + Base.iterate(I::($INDEX), state::Int = 1) = (state == length(I.idx) + 1) ? nothing : (I[state], state + 1) # Necessary to check if, e.g. `(cellid, faceidx) in faceset` Base.isequal(x::$INDEX, y::$INDEX) = x.idx == y.idx From 60442165a3b1b84eddd9a50e2e1b779b0f48bd1c Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 11 Nov 2024 16:49:19 +0100 Subject: [PATCH 08/28] now sets are moved from there to here --- src/iterators.jl | 33 ++++++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 9 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index d72e616c14..d6e9c92144 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -407,20 +407,35 @@ function InterfaceIterator( fs = facetskeleton(topology, grid) ninterfaces = 0 for facet in fs - facet[1] ∈ sdh_here.cellset || continue - neighbors = neighborhood[facet[1], facet[2]] - isempty(neighbors) && continue - neighbors[][1] ∈ sdh_there.cellset || continue + if facet[1] ∈ sdh_here.cellset + neighbors = neighborhood[facet[1], facet[2]] + isempty(neighbors) && continue + neighbors[][1] ∈ sdh_there.cellset || continue + elseif facet[1] ∈ sdh_there.cellset + neighbors = neighborhood[facet[1], facet[2]] + isempty(neighbors) && continue + neighbors[][1] ∈ sdh_here.cellset || continue + else + continue + end ninterfaces += 1 end set = Set{InterfaceIndex}() sizehint!(set, ninterfaces) for facet in fs - facet[1] ∈ sdh_here.cellset || continue - neighbors = neighborhood[facet[1], facet[2]] - isempty(neighbors) && continue - neighbors[][1] ∈ sdh_there.cellset || continue - push!(set, InterfaceIndex(facet[1], facet[2], neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2])) + if facet[1] ∈ sdh_here.cellset + neighbors = neighborhood[facet[1], facet[2]] + isempty(neighbors) && continue + neighbors[][1] ∈ sdh_there.cellset || continue + push!(set, InterfaceIndex(facet[1], facet[2], neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2])) + elseif facet[1] ∈ sdh_there.cellset + neighbors = neighborhood[facet[1], facet[2]] + isempty(neighbors) && continue + neighbors[][1] ∈ sdh_here.cellset || continue + push!(set, InterfaceIndex(neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2], facet[1], facet[2])) + else + continue + end end return InterfaceIterator(InterfaceCache(sdh_here, sdh_there), set) end From d0fc3b56f7548543dd5dad755efcc04c2ad6715c Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 11 Nov 2024 16:49:34 +0100 Subject: [PATCH 09/28] some tests --- test/runtests.jl | 1 + test/test_iterators.jl | 65 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 66 insertions(+) create mode 100644 test/test_iterators.jl diff --git a/test/runtests.jl b/test/runtests.jl index 57cd82d8a7..08889705da 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -20,6 +20,7 @@ include("test_interpolations.jl") include("test_cellvalues.jl") include("test_facevalues.jl") include("test_interfacevalues.jl") +include("test_iterators.jl") include("test_quadrules.jl") include("test_assemble.jl") include("test_dofs.jl") diff --git a/test/test_iterators.jl b/test/test_iterators.jl new file mode 100644 index 0000000000..a593b756f0 --- /dev/null +++ b/test/test_iterators.jl @@ -0,0 +1,65 @@ +@testset "InterfaceIterator" begin + function _iterate(iterator) + for _ in iterator end + return nothing + end + grid = generate_grid(Hexahedron, (1000, 2, 1)) + ip = Lagrange{RefHexahedron, 1}() + @testset "Single domain" begin + dh = DofHandler(grid) + add!(dh, :u, ip) + close!(dh) + topology = ExclusiveTopology(grid) + @testset "construction" begin + ii_dh_top = InterfaceIterator(dh, topology) + ii_dh = InterfaceIterator(dh) + ii_grid_top = InterfaceIterator(grid, topology) + ii_grid = InterfaceIterator(grid) + @test ii_dh.set == ii_dh_top.set == ii_grid.set == ii_grid_top.set + # Test that topology has no effect on iterator type + @test typeof(ii_dh) == typeof(ii_dh_top) + @test typeof(ii_grid) == typeof(ii_grid_top) + # Precompile the function for allocations test + _iterate(ii_dh) + _iterate(ii_grid) + # Test for allocations + @test (@allocations _iterate(ii_dh)) == 1 # Should be zero?? + @test (@allocations _iterate(ii_dh_top)) == 2 # Should be zero?? this one is 2??? + @test (@allocations _iterate(ii_grid)) == 1 # Should be zero?? + @test (@allocations _iterate(ii_grid_top)) == 1 # Should be zero?? + end + end + @testset "subdomains" begin + dh = DofHandler(grid) + sdh1 = SubDofHandler(dh, Set([collect(1:500)..., collect(1501:2000)...])) + sdh2 = SubDofHandler(dh, Set([collect(501:1500)...])) + add!(sdh1, :u, ip) + add!(sdh2, :u, ip) + close!(dh) + topology = ExclusiveTopology(grid) + @testset "construction" begin + ii_sdh_top = InterfaceIterator(sdh1, sdh2, topology) + ii_sdh = InterfaceIterator(sdh1, sdh2) + ii_same_sdh = InterfaceIterator(sdh1, sdh1) + @test ii_sdh.set == ii_sdh_top.set + @test length(ii_sdh.set) == 1002 + @test count(interface_index -> interface_index[1] == interface_index[3] - 1000, ii_sdh.set) == 500 + @test count(interface_index -> interface_index[1] == interface_index[3] + 1000, ii_sdh.set) == 500 + ii_sdh_flipped = InterfaceIterator(sdh2, sdh1) + @test all(interface -> InterfaceIndex(interface[3], interface[4], interface[1], interface[2]) ∈ ii_sdh_flipped.set, ii_sdh.set) + # Test that topology has no effect on iterator type + @test typeof(ii_sdh) == typeof(ii_sdh_top) + # Precompile the function for allocations test + _iterate(ii_sdh) + _iterate(ii_same_sdh) + # Test for allocations + @test (@allocations _iterate(ii_sdh)) == 1 # Should be zero?? + @test (@allocations _iterate(ii_same_sdh)) == 1 # Should be zero?? + end + end + @testset "InterfaceIndex" begin + # checkmate, codecov bot! + @test isequal(InterfaceIndex(1, 2, 3, 4), (1, 2, 3, 4)) + @test isequal((1, 2, 3, 4), InterfaceIndex(1, 2, 3, 4)) + end +end From 79d55f4d7de864f9df591d19ae3bbceb360bb4fe Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 11 Nov 2024 17:43:46 +0100 Subject: [PATCH 10/28] more tests --- src/iterators.jl | 37 ++++++++------ test/test_iterators.jl | 112 ++++++++++++++++++++++++++++++----------- 2 files changed, 107 insertions(+), 42 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index d6e9c92144..0fede9fd4e 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -357,25 +357,24 @@ end function InterfaceIterator( gridordh::Union{Grid, DofHandler}, - set_here::AbstractVecOrSet{FacetIndex}, - topology::ExclusiveTopology = ExclusiveTopology(gridordh isa Grid ? gridordh : get_grid(gridordh)) + set::AbstractVecOrSet{InterfaceIndex} ) grid = gridordh isa Grid ? gridordh : get_grid(gridordh) if gridordh isa DofHandler # Keep here to maintain same settings as for CellIterator - _check_same_celltype(grid, set_here) - end - neighborhood = get_facet_facet_neighborhood(topology, grid) - @assert all(facet -> !isempty(neighborhood[facet[1], facet[2]]), set_here) "Facets in the set must have neighbors" - set = Set{InterfaceIndex}() - sizehint!(set, length(set_here)) - for facet in set_here - neighbor = neighborhood[facet[1], facet[2]][] - push!(set, InterfaceIndex(facet[1], facet[2], neighbor[1], neighbor[2])) + _check_same_celltype(grid, set) end return InterfaceIterator(InterfaceCache(gridordh), set) end - +function InterfaceIterator( + sdh_here::SubDofHandler, + sdh_there::SubDofHandler, + set::AbstractVecOrSet{InterfaceIndex} + ) + grid = get_grid(sdh_here.dh) + _check_same_celltype(grid, set) + return InterfaceIterator(InterfaceCache(sdh_here, sdh_there), set) +end function InterfaceIterator( gridordh::Union{Grid, DofHandler}, topology::ExclusiveTopology = ExclusiveTopology(gridordh isa Grid ? gridordh : get_grid(gridordh)) @@ -394,7 +393,7 @@ function InterfaceIterator( isempty(neighborhood[facet[1], facet[2]]) && continue push!(set, InterfaceIndex(facet[1], facet[2], neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2])) end - return InterfaceIterator(InterfaceCache(gridordh), set) + return InterfaceIterator(gridordh, set) end function InterfaceIterator( @@ -437,7 +436,7 @@ function InterfaceIterator( continue end end - return InterfaceIterator(InterfaceCache(sdh_here, sdh_there), set) + return InterfaceIterator(sdh_here, sdh_there, set) end # Iterator interface for CellIterator/FacetIterator @@ -473,3 +472,13 @@ function _check_same_celltype(grid::AbstractGrid, facetset::AbstractVecOrSet{<:B end return end + +function _check_same_celltype(grid::AbstractGrid, interfaceset::AbstractVecOrSet{InterfaceIndex}) + isconcretetype(getcelltype(grid)) && return nothing # Short circuit check + celltype_here = getcelltype(grid, first(interfaceset)[1]) + celltype_there = getcelltype(grid, first(interfaceset)[3]) + if !all(getcelltype(grid, interface[1]) == celltype_here && getcelltype(grid, interface[3]) == celltype_there for interface in interfaceset) + error("The cells in the set (set of InterfaceIndex) are not all of the same celltype on each side.") + end + return +end diff --git a/test/test_iterators.jl b/test/test_iterators.jl index a593b756f0..47b88e0707 100644 --- a/test/test_iterators.jl +++ b/test/test_iterators.jl @@ -3,41 +3,44 @@ for _ in iterator end return nothing end - grid = generate_grid(Hexahedron, (1000, 2, 1)) - ip = Lagrange{RefHexahedron, 1}() @testset "Single domain" begin + grid = generate_grid(Hexahedron, (1000, 2, 1)) + ip = Lagrange{RefHexahedron, 1}() dh = DofHandler(grid) add!(dh, :u, ip) close!(dh) topology = ExclusiveTopology(grid) - @testset "construction" begin - ii_dh_top = InterfaceIterator(dh, topology) - ii_dh = InterfaceIterator(dh) - ii_grid_top = InterfaceIterator(grid, topology) - ii_grid = InterfaceIterator(grid) - @test ii_dh.set == ii_dh_top.set == ii_grid.set == ii_grid_top.set - # Test that topology has no effect on iterator type - @test typeof(ii_dh) == typeof(ii_dh_top) - @test typeof(ii_grid) == typeof(ii_grid_top) - # Precompile the function for allocations test - _iterate(ii_dh) - _iterate(ii_grid) - # Test for allocations - @test (@allocations _iterate(ii_dh)) == 1 # Should be zero?? - @test (@allocations _iterate(ii_dh_top)) == 2 # Should be zero?? this one is 2??? - @test (@allocations _iterate(ii_grid)) == 1 # Should be zero?? - @test (@allocations _iterate(ii_grid_top)) == 1 # Should be zero?? - end + + ii_dh_top = InterfaceIterator(dh, topology) + ii_dh = InterfaceIterator(dh) + ii_grid_top = InterfaceIterator(grid, topology) + ii_grid = InterfaceIterator(grid) + @test ii_dh.set == ii_dh_top.set == ii_grid.set == ii_grid_top.set + # Test that topology has no effect on iterator type + @test typeof(ii_dh) == typeof(ii_dh_top) + @test typeof(ii_grid) == typeof(ii_grid_top) + # Precompile the function for allocations test + _iterate(ii_dh) + _iterate(ii_grid) + # Test for allocations + @test (@allocations _iterate(ii_dh)) == 1 # Should be zero?? + @test (@allocations _iterate(ii_dh_top)) == 2 # Should be zero?? this one is 2??? + @test (@allocations _iterate(ii_grid)) == 1 # Should be zero?? + @test (@allocations _iterate(ii_grid_top)) == 1 # Should be zero?? end + @testset "subdomains" begin - dh = DofHandler(grid) - sdh1 = SubDofHandler(dh, Set([collect(1:500)..., collect(1501:2000)...])) - sdh2 = SubDofHandler(dh, Set([collect(501:1500)...])) - add!(sdh1, :u, ip) - add!(sdh2, :u, ip) - close!(dh) - topology = ExclusiveTopology(grid) - @testset "construction" begin + @testset "same cell type" begin + grid = generate_grid(Hexahedron, (1000, 2, 1)) + ip = Lagrange{RefHexahedron, 1}() + dh = DofHandler(grid) + sdh1 = SubDofHandler(dh, Set([collect(1:500)..., collect(1501:2000)...])) + sdh2 = SubDofHandler(dh, Set([collect(501:1500)...])) + add!(sdh1, :u, ip) + add!(sdh2, :u, ip) + close!(dh) + topology = ExclusiveTopology(grid) + ii_sdh_top = InterfaceIterator(sdh1, sdh2, topology) ii_sdh = InterfaceIterator(sdh1, sdh2) ii_same_sdh = InterfaceIterator(sdh1, sdh1) @@ -56,10 +59,63 @@ @test (@allocations _iterate(ii_sdh)) == 1 # Should be zero?? @test (@allocations _iterate(ii_same_sdh)) == 1 # Should be zero?? end + + @testset "mixed cell types" begin + nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] + cells = [ + Quadrilateral((1, 2, 5, 4)), + Triangle((3, 5, 2)), + ] + grid = Grid(cells, nodes) + topology = ExclusiveTopology(grid) + ip_quad = Lagrange{RefQuadrilateral, 1}() + ip_tri = Lagrange{RefTriangle, 1}() + dh = DofHandler(grid) + sdh_quad = SubDofHandler(dh, Set([1])) + sdh_tri = SubDofHandler(dh, Set([2])) + add!(sdh_quad, :u, ip_quad) + add!(sdh_tri, :u, ip_tri) + close!(dh) + + ii_sdh_top = InterfaceIterator(sdh_quad, sdh_tri, topology) + ii_sdh = InterfaceIterator(sdh_quad, sdh_tri) + @test ii_sdh.set == ii_sdh_top.set + @test ii_sdh.set == Set([InterfaceIndex(1, 2, 2, 2)]) + ii_sdh_flipped = InterfaceIterator(sdh_tri, sdh_quad) + @test ii_sdh_flipped.set == Set([InterfaceIndex(2, 2, 1, 2)]) + # Test that topology has no effect on iterator type + @test typeof(ii_sdh) == typeof(ii_sdh_top) + # Precompile the function for allocations test + _iterate(ii_sdh) + # Test for allocations + @test (@allocations _iterate(ii_sdh)) == 1 # Should be zero?? + end end + @testset "InterfaceIndex" begin # checkmate, codecov bot! @test isequal(InterfaceIndex(1, 2, 3, 4), (1, 2, 3, 4)) @test isequal((1, 2, 3, 4), InterfaceIndex(1, 2, 3, 4)) end + + @testset "error paths" begin + nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] + cells = [ + Quadrilateral((1, 2, 5, 4)), + Triangle((3, 5, 2)), + ] + grid = Grid(cells, nodes) + topology = ExclusiveTopology(grid) + ip_quad = Lagrange{RefQuadrilateral, 1}() + ip_tri = Lagrange{RefTriangle, 1}() + dh = DofHandler(grid) + sdh_quad = SubDofHandler(dh, Set([1])) + sdh_tri = SubDofHandler(dh, Set([2])) + add!(sdh_quad, :u, ip_quad) + add!(sdh_tri, :u, ip_tri) + close!(dh) + + # @test_throws ErrorException InterfaceIterator(sdh_quad, sdh_tri, Set([InterfaceIndex(2, 2, 1, 2)])) + + end end From 9b1cbdd26ec4767525356ca3b7168cfaae78629b Mon Sep 17 00:00:00 2001 From: Abdulaziz Hamid Mohamed Date: Mon, 11 Nov 2024 23:36:19 +0100 Subject: [PATCH 11/28] error paths + trying to test for allocs --- Project.toml | 1 + src/iterators.jl | 8 ++++---- test/runtests.jl | 1 + test/test_iterators.jl | 16 +++++++++------- 4 files changed, 15 insertions(+), 11 deletions(-) diff --git a/Project.toml b/Project.toml index e3c5eb3fbd..738ead2a9b 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,7 @@ uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992" version = "1.0.0" [deps] +AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a" EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/iterators.jl b/src/iterators.jl index 0fede9fd4e..4bef73f1f2 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -372,7 +372,7 @@ function InterfaceIterator( set::AbstractVecOrSet{InterfaceIndex} ) grid = get_grid(sdh_here.dh) - _check_same_celltype(grid, set) + _check_same_celltype(grid, set, getcelltype(grid, first(sdh_here.cellset)), getcelltype(grid, first(sdh_there.cellset))) return InterfaceIterator(InterfaceCache(sdh_here, sdh_there), set) end function InterfaceIterator( @@ -473,10 +473,10 @@ function _check_same_celltype(grid::AbstractGrid, facetset::AbstractVecOrSet{<:B return end -function _check_same_celltype(grid::AbstractGrid, interfaceset::AbstractVecOrSet{InterfaceIndex}) +function _check_same_celltype(grid::AbstractGrid, interfaceset::AbstractVecOrSet{InterfaceIndex}, + celltype_here::Type{<:AbstractCell} = getcelltype(grid, first(interfaceset)[1]), + celltype_there::Type{<:AbstractCell} = getcelltype(grid, first(interfaceset)[3])) isconcretetype(getcelltype(grid)) && return nothing # Short circuit check - celltype_here = getcelltype(grid, first(interfaceset)[1]) - celltype_there = getcelltype(grid, first(interfaceset)[3]) if !all(getcelltype(grid, interface[1]) == celltype_here && getcelltype(grid, interface[3]) == celltype_there for interface in interfaceset) error("The cells in the set (set of InterfaceIndex) are not all of the same celltype on each side.") end diff --git a/test/runtests.jl b/test/runtests.jl index 08889705da..0f339c752c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,6 +11,7 @@ using StaticArrays using OrderedCollections using WriteVTK import Metis +using AllocCheck include("test_utils.jl") diff --git a/test/test_iterators.jl b/test/test_iterators.jl index 47b88e0707..d5ab2c601d 100644 --- a/test/test_iterators.jl +++ b/test/test_iterators.jl @@ -1,4 +1,8 @@ @testset "InterfaceIterator" begin + @check_allocs function testallocs_iterate(iterator) + for _ in iterator end + return nothing + end function _iterate(iterator) for _ in iterator end return nothing @@ -23,10 +27,10 @@ _iterate(ii_dh) _iterate(ii_grid) # Test for allocations - @test (@allocations _iterate(ii_dh)) == 1 # Should be zero?? - @test (@allocations _iterate(ii_dh_top)) == 2 # Should be zero?? this one is 2??? - @test (@allocations _iterate(ii_grid)) == 1 # Should be zero?? - @test (@allocations _iterate(ii_grid_top)) == 1 # Should be zero?? + @test (@allocated _iterate(ii_dh)) == 0 + @test (@allocated _iterate(ii_dh_top)) == 144 # Why??? + @test (@allocated _iterate(ii_grid)) == 0 + @test (@allocated _iterate(ii_grid_top)) == 0 end @testset "subdomains" begin @@ -114,8 +118,6 @@ add!(sdh_quad, :u, ip_quad) add!(sdh_tri, :u, ip_tri) close!(dh) - - # @test_throws ErrorException InterfaceIterator(sdh_quad, sdh_tri, Set([InterfaceIndex(2, 2, 1, 2)])) - + @test_throws ErrorException("The cells in the set (set of InterfaceIndex) are not all of the same celltype on each side.") InterfaceIterator(sdh_quad, sdh_tri, Set([InterfaceIndex(2, 2, 1, 2)])) end end From be27c621d4ed4cfd81605ab0af0513a03e0324b1 Mon Sep 17 00:00:00 2001 From: Abdulaziz Hamid Mohamed Date: Mon, 11 Nov 2024 23:41:41 +0100 Subject: [PATCH 12/28] Make Runic happy --- src/iterators.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index 4bef73f1f2..5b13f3a2a3 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -473,9 +473,11 @@ function _check_same_celltype(grid::AbstractGrid, facetset::AbstractVecOrSet{<:B return end -function _check_same_celltype(grid::AbstractGrid, interfaceset::AbstractVecOrSet{InterfaceIndex}, +function _check_same_celltype( + grid::AbstractGrid, interfaceset::AbstractVecOrSet{InterfaceIndex}, celltype_here::Type{<:AbstractCell} = getcelltype(grid, first(interfaceset)[1]), - celltype_there::Type{<:AbstractCell} = getcelltype(grid, first(interfaceset)[3])) + celltype_there::Type{<:AbstractCell} = getcelltype(grid, first(interfaceset)[3]) + ) isconcretetype(getcelltype(grid)) && return nothing # Short circuit check if !all(getcelltype(grid, interface[1]) == celltype_here && getcelltype(grid, interface[3]) == celltype_there for interface in interfaceset) error("The cells in the set (set of InterfaceIndex) are not all of the same celltype on each side.") From cf7c08e09eef172f29f125e0f3332088b48c6ade Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 12 Nov 2024 13:09:35 +0100 Subject: [PATCH 13/28] make alloccheck happy --- src/iterators.jl | 27 ++++++++++++++++++++++----- test/test_iterators.jl | 2 +- 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index 5b13f3a2a3..5bce5daf4b 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -61,12 +61,16 @@ function CellCache(dh::DofHandler{dim}, flags::UpdateFlags = UpdateFlags()) wher return CellCache(flags, get_grid(dh), -1, nodes, coords, dh, celldofs) end -function CellCache(sdh::SubDofHandler, flags::UpdateFlags = UpdateFlags()) - Tv = get_coordinate_type(sdh.dh.grid) - return CellCache(flags, sdh.dh.grid, -1, Int[], Tv[], sdh, Int[]) +function CellCache(sdh::SubDofHandler{<:DofHandler{dim}}, flags::UpdateFlags = UpdateFlags()) where {dim} + n = ndofs_per_cell(sdh) # dofs and coords will be resized in `reinit!` + N = nnodes_per_cell(get_grid(sdh.dh), sdh.cellset[1]) + nodes = zeros(Int, N) + coords = zeros(Vec{dim, get_coordinate_eltype(get_grid(sdh.dh))}, N) + celldofs = zeros(Int, n) + return CellCache(flags, sdh.dh.grid, -1, nodes, coords, sdh.dh, celldofs) end -function reinit!(cc::CellCache, i::Int) +function reinit!(cc::CellCache{<:Any, <:Grid{<:Any, AbstractCell}}, i::Int) cc.cellid = i if cc.flags.nodes resize!(cc.nodes, nnodes_per_cell(cc.grid, i)) @@ -83,6 +87,19 @@ function reinit!(cc::CellCache, i::Int) return cc end +function reinit!(cc::CellCache, i::Int) + cc.cellid = i + if cc.flags.nodes + cellnodes!(cc.nodes, cc.grid, i) + end + if cc.flags.coords + getcoordinates!(cc.coords, cc.grid, i) + end + if cc.dh !== nothing && cc.flags.dofs + celldofs!(cc.dofs, cc.dh, i) + end + return cc +end # reinit! FEValues with CellCache reinit!(cv::CellValues, cc::CellCache) = reinit!(cv, cc.coords) reinit!(fv::FacetValues, cc::CellCache, f::Int) = reinit!(fv, cc.coords, f) @@ -173,7 +190,7 @@ end function InterfaceCache(gridordh::Union{AbstractGrid, AbstractDofHandler}) fc_a = FacetCache(gridordh) fc_b = FacetCache(gridordh) - return InterfaceCache(fc_a, fc_b, Int[]) + return InterfaceCache(fc_a, fc_b, zeros(Int, length(celldofs(fc_a)) + length(celldofs(fc_b)))) end function InterfaceCache(sdh_here::SubDofHandler, sdh_there::SubDofHandler) diff --git a/test/test_iterators.jl b/test/test_iterators.jl index d5ab2c601d..4e1d572f03 100644 --- a/test/test_iterators.jl +++ b/test/test_iterators.jl @@ -28,7 +28,7 @@ _iterate(ii_grid) # Test for allocations @test (@allocated _iterate(ii_dh)) == 0 - @test (@allocated _iterate(ii_dh_top)) == 144 # Why??? + @test (@allocated _iterate(ii_dh_top)) == 0 @test (@allocated _iterate(ii_grid)) == 0 @test (@allocated _iterate(ii_grid_top)) == 0 end From 0b68251caef1b112a90497ab80a22384d66e2dbe Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 12 Nov 2024 13:16:33 +0100 Subject: [PATCH 14/28] revert oopsie --- src/iterators.jl | 15 +-------------- 1 file changed, 1 insertion(+), 14 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index 5bce5daf4b..21bd53d95c 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -70,7 +70,7 @@ function CellCache(sdh::SubDofHandler{<:DofHandler{dim}}, flags::UpdateFlags = U return CellCache(flags, sdh.dh.grid, -1, nodes, coords, sdh.dh, celldofs) end -function reinit!(cc::CellCache{<:Any, <:Grid{<:Any, AbstractCell}}, i::Int) +function reinit!(cc::CellCache, i::Int) cc.cellid = i if cc.flags.nodes resize!(cc.nodes, nnodes_per_cell(cc.grid, i)) @@ -87,19 +87,6 @@ function reinit!(cc::CellCache{<:Any, <:Grid{<:Any, AbstractCell}}, i::Int) return cc end -function reinit!(cc::CellCache, i::Int) - cc.cellid = i - if cc.flags.nodes - cellnodes!(cc.nodes, cc.grid, i) - end - if cc.flags.coords - getcoordinates!(cc.coords, cc.grid, i) - end - if cc.dh !== nothing && cc.flags.dofs - celldofs!(cc.dofs, cc.dh, i) - end - return cc -end # reinit! FEValues with CellCache reinit!(cv::CellValues, cc::CellCache) = reinit!(cv, cc.coords) reinit!(fv::FacetValues, cc::CellCache, f::Int) = reinit!(fv, cc.coords, f) From 19f11ea538ecfd88fed7c2135ad8fd19ef4e88a9 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 12 Nov 2024 14:02:18 +0100 Subject: [PATCH 15/28] next: get rid of dynamic dispatches --- src/iterators.jl | 27 +++++++++++++++++++++------ test/test_iterators.jl | 34 +++++++++++++++------------------- 2 files changed, 36 insertions(+), 25 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index 21bd53d95c..97e50920b3 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -70,7 +70,8 @@ function CellCache(sdh::SubDofHandler{<:DofHandler{dim}}, flags::UpdateFlags = U return CellCache(flags, sdh.dh.grid, -1, nodes, coords, sdh.dh, celldofs) end -function reinit!(cc::CellCache, i::Int) +# TODO: Find a better way to make AllocCheck.jl happy +function reinit!(cc::CellCache{<:Any, <:Any, <:Nothing}, i::Int) cc.cellid = i if cc.flags.nodes resize!(cc.nodes, nnodes_per_cell(cc.grid, i)) @@ -87,6 +88,21 @@ function reinit!(cc::CellCache, i::Int) return cc end +function reinit!(cc::CellCache{<:Any, <:Any, <:AbstractDofHandler}, i::Int) + # If we have a DofHandler the cells must be of the same type -> no need to resize + cc.cellid = i + if cc.flags.nodes + cellnodes!(cc.nodes, cc.grid, i) + end + if cc.flags.coords + getcoordinates!(cc.coords, cc.grid, i) + end + if cc.dh !== nothing && cc.flags.dofs + celldofs!(cc.dofs, cc.dh, i) + end + return cc +end + # reinit! FEValues with CellCache reinit!(cv::CellValues, cc::CellCache) = reinit!(cv, cc.coords) reinit!(fv::FacetValues, cc::CellCache, f::Int) = reinit!(fv, cc.coords, f) @@ -168,9 +184,9 @@ interface. The cache is updated for a new cell by calling `reinit!(cache, facet_ See also [`InterfaceIterator`](@ref). """ -struct InterfaceCache{FC <: FacetCache} - a::FC - b::FC +struct InterfaceCache{FC_HERE <: FacetCache, FC_THERE <: FacetCache} + a::FC_HERE + b::FC_THERE dofs::Vector{Int} end @@ -183,13 +199,12 @@ end function InterfaceCache(sdh_here::SubDofHandler, sdh_there::SubDofHandler) fc_a = FacetCache(sdh_here) fc_b = FacetCache(sdh_there) - return InterfaceCache(fc_a, fc_b, Int[]) + return InterfaceCache(fc_a, fc_b, zeros(Int, length(celldofs(fc_a)) + length(celldofs(fc_b)))) end function reinit!(cache::InterfaceCache, interface::InterfaceIndex) reinit!(cache.a, FacetIndex(interface.idx[1], interface.idx[2])) reinit!(cache.b, FacetIndex(interface.idx[3], interface.idx[4])) - resize!(cache.dofs, length(celldofs(cache.a)) + length(celldofs(cache.b))) for (i, d) in pairs(cache.a.dofs) cache.dofs[i] = d end diff --git a/test/test_iterators.jl b/test/test_iterators.jl index 4e1d572f03..5a8ca98891 100644 --- a/test/test_iterators.jl +++ b/test/test_iterators.jl @@ -3,10 +3,7 @@ for _ in iterator end return nothing end - function _iterate(iterator) - for _ in iterator end - return nothing - end + @testset "Single domain" begin grid = generate_grid(Hexahedron, (1000, 2, 1)) ip = Lagrange{RefHexahedron, 1}() @@ -23,14 +20,12 @@ # Test that topology has no effect on iterator type @test typeof(ii_dh) == typeof(ii_dh_top) @test typeof(ii_grid) == typeof(ii_grid_top) - # Precompile the function for allocations test - _iterate(ii_dh) - _iterate(ii_grid) # Test for allocations - @test (@allocated _iterate(ii_dh)) == 0 - @test (@allocated _iterate(ii_dh_top)) == 0 - @test (@allocated _iterate(ii_grid)) == 0 - @test (@allocated _iterate(ii_grid_top)) == 0 + @test (testallocs_iterate(ii_dh_top)) === nothing + @test (testallocs_iterate(ii_dh)) === nothing + # Iterators over grid can allocate due to potential resize! + @test_throws AllocCheckFailure testallocs_iterate(ii_grid) + @test_throws AllocCheckFailure testallocs_iterate(ii_grid_top) end @testset "subdomains" begin @@ -56,12 +51,9 @@ @test all(interface -> InterfaceIndex(interface[3], interface[4], interface[1], interface[2]) ∈ ii_sdh_flipped.set, ii_sdh.set) # Test that topology has no effect on iterator type @test typeof(ii_sdh) == typeof(ii_sdh_top) - # Precompile the function for allocations test - _iterate(ii_sdh) - _iterate(ii_same_sdh) # Test for allocations - @test (@allocations _iterate(ii_sdh)) == 1 # Should be zero?? - @test (@allocations _iterate(ii_same_sdh)) == 1 # Should be zero?? + @test (testallocs_iterate(ii_sdh)) === nothing + @test (testallocs_iterate(ii_same_sdh)) === nothing end @testset "mixed cell types" begin @@ -89,10 +81,14 @@ @test ii_sdh_flipped.set == Set([InterfaceIndex(2, 2, 1, 2)]) # Test that topology has no effect on iterator type @test typeof(ii_sdh) == typeof(ii_sdh_top) - # Precompile the function for allocations test - _iterate(ii_sdh) # Test for allocations - @test (@allocations _iterate(ii_sdh)) == 1 # Should be zero?? + # try + # testallocs_iterate(ii_sdh_top) + # catch err + # @info err.errors[1] + # end + # @test (testallocs_iterate(ii_sdh_top)) === nothing + # @test (testallocs_iterate(ii_sdh)) === nothing end end From 84a1a14c7476064052a9a1441c0cd8273ce2b8f6 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 12 Nov 2024 22:44:11 +0100 Subject: [PATCH 16/28] home office push --- src/Dofs/ConstraintHandler.jl | 20 +++--- src/Dofs/DofRenumbering.jl | 2 +- src/Dofs/sparsity_pattern.jl | 120 +++++++++++++++++----------------- src/iterators.jl | 1 + 4 files changed, 73 insertions(+), 70 deletions(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index a48cb3940f..dd9f8d23b2 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -356,16 +356,18 @@ function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::AbstractVecOrSet{ interpol_points = getnbasefunctions(interpolation) node_dofs = zeros(Int, ncomps, nnodes) visited = falses(nnodes) - for cell in CellIterator(ch.dh, cellset) # only go over cells that belong to current SubDofHandler - for idx in 1:min(interpol_points, length(cell.nodes)) - node = cell.nodes[idx] - if !visited[node] - noderange = (offset + (idx - 1) * field_dim + 1):(offset + idx * field_dim) # the dofs in this node - for (i, c) in enumerate(dbc.components) - node_dofs[i, node] = cell.dofs[noderange[c]] - @debug println("adding dof $(cell.dofs[noderange[c]]) to node_dofs") + for sdh in ch.dh.subdofhandlers + for cell in CellIterator(sdh, cellset) # only go over cells that belong to current SubDofHandler + for idx in 1:min(interpol_points, length(cell.nodes)) + node = cell.nodes[idx] + if !visited[node] + noderange = (offset + (idx - 1) * field_dim + 1):(offset + idx * field_dim) # the dofs in this node + for (i, c) in enumerate(dbc.components) + node_dofs[i, node] = cell.dofs[noderange[c]] + @debug println("adding dof $(cell.dofs[noderange[c]]) to node_dofs") + end + visited[node] = true end - visited[node] = true end end end diff --git a/src/Dofs/DofRenumbering.jl b/src/Dofs/DofRenumbering.jl index 0d1a2f70b0..83e44615fb 100644 --- a/src/Dofs/DofRenumbering.jl +++ b/src/Dofs/DofRenumbering.jl @@ -196,7 +196,7 @@ function compute_renumber_permutation(dh::DofHandler, _, order::DofOrder.Compone for sdh in dh.subdofhandlers dof_ranges = [dof_range(sdh, f) for f in eachindex(sdh.field_names)] global_idxs = [findfirst(x -> x === f, dh.field_names) for f in sdh.field_names] - for cell in CellIterator(dh, sdh.cellset, flags) + for cell in CellIterator(sdh, flags) cdofs = celldofs(cell) for (local_idx, global_idx) in pairs(global_idxs) rng = dof_ranges[local_idx] diff --git a/src/Dofs/sparsity_pattern.jl b/src/Dofs/sparsity_pattern.jl index 95f80841cc..d9db08f780 100644 --- a/src/Dofs/sparsity_pattern.jl +++ b/src/Dofs/sparsity_pattern.jl @@ -446,34 +446,39 @@ function _coupling_to_local_dof_coupling(dh::DofHandler, coupling::AbstractMatri sz == size(coupling, 2) || error("coupling not square") # Return one matrix per (potential) sub-domain - outs = Matrix{Bool}[] + nsdh = length(dh.subdofhandlers) + outs = Matrix{Matrix{Bool}}(undef, nsdh, nsdh) field_dims = map(fieldname -> n_components(dh, fieldname), dh.field_names) - for sdh in dh.subdofhandlers - out = zeros(Bool, ndofs_per_cell(sdh), ndofs_per_cell(sdh)) - push!(outs, out) + for (sdh1_idx, sdh1) in enumerate(dh.subdofhandlers) + for (sdh2_idx, sdh2) in enumerate(dh.subdofhandlers) + out = zeros(Bool, ndofs_per_cell(sdh1), ndofs_per_cell(sdh2)) + outs[sdh1_idx, sdh2_idx] = out - dof_ranges = [dof_range(sdh, f) for f in sdh.field_names] - global_idxs = [findfirst(x -> x === f, dh.field_names) for f in sdh.field_names] + dof_ranges_here = [dof_range(sdh1, f) for f in sdh1.field_names] + dof_ranges_there = [dof_range(sdh2, f) for f in sdh2.field_names] + global_idxs_here = [findfirst(x -> x === f, dh.field_names) for f in sdh1.field_names] + global_idxs_there = [findfirst(x -> x === f, dh.field_names) for f in sdh2.field_names] - if sz == length(dh.field_names) # Coupling given by fields - for (j, jrange) in pairs(dof_ranges), (i, irange) in pairs(dof_ranges) - out[irange, jrange] .= coupling[global_idxs[i], global_idxs[j]] - end - elseif sz == sum(field_dims) # Coupling given by components - component_offsets = pushfirst!(cumsum(field_dims), 0) - for (jf, jrange) in pairs(dof_ranges), (j, J) in pairs(jrange) - jc = mod1(j, field_dims[global_idxs[jf]]) + component_offsets[global_idxs[jf]] - for (i_f, irange) in pairs(dof_ranges), (i, I) in pairs(irange) - ic = mod1(i, field_dims[global_idxs[i_f]]) + component_offsets[global_idxs[i_f]] - out[I, J] = coupling[ic, jc] + if sz == length(dh.field_names) # Coupling given by fields + for (j, jrange) in pairs(dof_ranges_there), (i, irange) in pairs(dof_ranges_here) + out[irange, jrange] .= coupling[global_idxs_here[i], global_idxs_there[j]] + end + elseif sz == sum(field_dims) # Coupling given by components + component_offsets = pushfirst!(cumsum(field_dims), 0) + for (jf, jrange) in pairs(dof_ranges_there), (j, J) in pairs(jrange) + jc = mod1(j, field_dims[global_idxs_there[jf]]) + component_offsets[global_idxs_there[jf]] + for (i_f, irange) in pairs(dof_ranges_here), (i, I) in pairs(irange) + ic = mod1(i, field_dims[global_idxs_here[i_f]]) + component_offsets[global_idxs_here[i_f]] + out[I, J] = coupling[ic, jc] + end end + elseif sz == ndofs_per_cell(sdh1) # Coupling given by template local matrix + # TODO: coupling[fieldhandler_idx] if different template per subddomain + out .= coupling + else + error("could not create coupling") end - elseif sz == ndofs_per_cell(sdh) # Coupling given by template local matrix - # TODO: coupling[fieldhandler_idx] if different template per subddomain - out .= coupling - else - error("could not create coupling") end end return outs @@ -481,16 +486,13 @@ end function _add_cell_entries!( sp::AbstractSparsityPattern, dh::DofHandler, ch::Union{ConstraintHandler, Nothing}, - keep_constrained::Bool, coupling::Union{Vector{<:AbstractMatrix{Bool}}, Nothing}, + keep_constrained::Bool, coupling::Union{<:AbstractMatrix{<:AbstractMatrix{Bool}}, Nothing}, ) # Add all connections between dofs for every cell while filtering based # on a) constraints, and b) field/dof coupling. - cc = CellCache(dh) for (sdhi, sdh) in pairs(dh.subdofhandlers) - set = BitSet(sdh.cellset) coupling === nothing || (coupling_sdh = coupling[sdhi]) - for cell_id in set - reinit!(cc, cell_id) + for cc in CellIterator(sdh) for (i, row) in pairs(cc.dofs) # a) check constraint for row !keep_constrained && haskey(ch.dofmapping, row) && continue @@ -606,41 +608,39 @@ function _add_interface_entries!( interface_coupling::AbstractMatrix{Bool}, ) couplings = _coupling_to_local_dof_coupling(dh, interface_coupling) - for ic in InterfaceIterator(dh, topology) - # TODO: This looks like it can be optimized for the common case where - # the cells are in the same subdofhandler - sdhs_idx = dh.cell_to_subdofhandler[cellid.([ic.a, ic.b])] - sdhs = dh.subdofhandlers[sdhs_idx] - to_check = Dict{Int, Vector{Int}}() - for (i, sdh) in pairs(sdhs) - sdh_idx = sdhs_idx[i] - coupling_sdh = couplings[sdh_idx] - for cell_field in sdh.field_names - dofrange1 = dof_range(sdh, cell_field) - cell_dofs = celldofs(sdh_idx == 1 ? ic.a : ic.b) - cell_field_dofs = @view cell_dofs[dofrange1] - for neighbor_field in sdh.field_names - sdh2 = sdhs[i == 1 ? 2 : 1] - neighbor_field ∈ sdh2.field_names || continue - dofrange2 = dof_range(sdh2, neighbor_field) - neighbor_dofs = celldofs(sdh_idx == 2 ? ic.a : ic.b) - neighbor_field_dofs = @view neighbor_dofs[dofrange2] - - empty!(to_check) - for (j, dof_j) in enumerate(dofrange2) - for (i, dof_i) in enumerate(dofrange1) - coupling_sdh[dof_i, dof_j] || continue - push!(get!(Vector{Int}, to_check, j), i) + for (_i, sdh1) in enumerate(dh.subdofhandlers) + for (_j, sdh2) in enumerate(dh.subdofhandlers) + for ic in InterfaceIterator(sdh1, sdh2, topology) + # TODO: This looks like it can be optimized for the common case where + # the cells are in the same subdofhandler + to_check = Dict{Int, Vector{Int}}() + coupling_sdh = couplings[_i, _j] + for cell_field in sdh1.field_names + dofrange1 = dof_range(sdh1, cell_field) + cell_dofs = celldofs(ic.a) + cell_field_dofs = @view cell_dofs[dofrange1] + for neighbor_field in sdh1.field_names + neighbor_field ∈ sdh2.field_names || continue + dofrange2 = dof_range(sdh2, neighbor_field) + neighbor_dofs = celldofs(ic.b) + neighbor_field_dofs = @view neighbor_dofs[dofrange2] + + empty!(to_check) + for (j, dof_j) in enumerate(dofrange2) + for (i, dof_i) in enumerate(dofrange1) + coupling_sdh[dof_i, dof_j] || continue + push!(get!(Vector{Int}, to_check, j), i) + end end - end - for (j, is) in to_check - # Avoid coupling the shared dof in continuous interpolations as cross-element. They're coupled in the local coupling matrix. - neighbor_field_dofs[j] ∈ cell_dofs && continue - for i in is - cell_field_dofs[i] ∈ neighbor_dofs && continue - _add_interface_entry(sp, cell_field_dofs, neighbor_field_dofs, i, j, keep_constrained, ch) - _add_interface_entry(sp, neighbor_field_dofs, cell_field_dofs, j, i, keep_constrained, ch) + for (j, is) in to_check + # Avoid coupling the shared dof in continuous interpolations as cross-element. They're coupled in the local coupling matrix. + neighbor_field_dofs[j] ∈ cell_dofs && continue + for i in is + cell_field_dofs[i] ∈ neighbor_dofs && continue + _add_interface_entry(sp, cell_field_dofs, neighbor_field_dofs, i, j, keep_constrained, ch) + _add_interface_entry(sp, neighbor_field_dofs, cell_field_dofs, j, i, keep_constrained, ch) + end end end end diff --git a/src/iterators.jl b/src/iterators.jl index 97e50920b3..8cb7a096e2 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -402,6 +402,7 @@ function InterfaceIterator( if gridordh isa DofHandler # Keep here to maintain same settings as for CellIterator _check_same_celltype(grid, 1:getncells(grid)) + @assert length(gridordh.subdofhandlers) == 1 "Use InterfaceIterator(::SubDofHandler, ::SubDofHandler) for subdomain support" end neighborhood = get_facet_facet_neighborhood(topology, grid) fs = facetskeleton(topology, grid) From 261c50ecae73bf347e1c3c684d841f7d38fe54f6 Mon Sep 17 00:00:00 2001 From: Abdulaziz Hamid Mohamed Date: Wed, 13 Nov 2024 01:24:07 +0100 Subject: [PATCH 17/28] Bad workarounds, think about them tomorrow :P --- src/Dofs/ConstraintHandler.jl | 72 +++++++++++++++++--------------- test/test_grid_dofhandler_vtk.jl | 2 +- test/test_mixeddofhandler.jl | 20 +++++---- 3 files changed, 53 insertions(+), 41 deletions(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index dd9f8d23b2..81c1665499 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -311,12 +311,15 @@ function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfacets::AbstractVecOrSet # loop over all the faces in the set and add the global dofs to `constrained_dofs` constrained_dofs = Int[] - cc = CellCache(ch.dh, UpdateFlags(; nodes = false, coords = false, dofs = true)) - for (cellidx, facetidx) in bcfacets - reinit!(cc, cellidx) - r = local_facet_dofs_offset[facetidx]:(local_facet_dofs_offset[facetidx + 1] - 1) - append!(constrained_dofs, cc.dofs[local_facet_dofs[r]]) # TODO: for-loop over r and simply push! to ch.prescribed_dofs - @debug println("adding dofs $(cc.dofs[local_facet_dofs[r]]) to dbc") + for sdh in ch.dh.subdofhandlers + facets_sdh_ = filter(x -> x[1] ∈ sdh.cellset, bcfacets) + facets_sdh = OrderedSet{FacetIndex}(FacetIndex(x[1], x[2]) for x in facets_sdh_) + for cc in FacetIterator(sdh, facets_sdh, UpdateFlags(; nodes = false, coords = false, dofs = true)) + facetidx = cc.current_facet_id + r = local_facet_dofs_offset[facetidx]:(local_facet_dofs_offset[facetidx + 1] - 1) + append!(constrained_dofs, cc.dofs[local_facet_dofs[r]]) # TODO: for-loop over r and simply push! to ch.prescribed_dofs + @debug println("adding dofs $(cc.dofs[local_facet_dofs[r]]) to dbc") + end end # save it to the ConstraintHandler @@ -357,7 +360,8 @@ function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::AbstractVecOrSet{ node_dofs = zeros(Int, ncomps, nnodes) visited = falses(nnodes) for sdh in ch.dh.subdofhandlers - for cell in CellIterator(sdh, cellset) # only go over cells that belong to current SubDofHandler + for cell in CellIterator(sdh) # only go over cells that belong to current SubDofHandler + cellid(cell) ∈ cellset || continue for idx in 1:min(interpol_points, length(cell.nodes)) node = cell.nodes[idx] if !visited[node] @@ -447,32 +451,34 @@ function _update!( dofmapping::Dict{Int, Int}, dofcoefficients::Vector{Union{Nothing, DofCoefficients{T}}}, time::Real ) where {T} - cc = CellCache(dh, UpdateFlags(; nodes = false, coords = true, dofs = true)) - for (cellidx, entityidx) in boundary_entities - reinit!(cc, cellidx) - - # no need to reinit!, enough to update current_entity since we only need geometric shape functions M - boundaryvalues.current_entity = entityidx - - # local dof-range for this facet - r = local_facet_dofs_offset[entityidx]:(local_facet_dofs_offset[entityidx + 1] - 1) - counter = 1 - for location in 1:getnquadpoints(boundaryvalues) - x = spatial_coordinate(boundaryvalues, location, cc.coords) - bc_value = f(x, time) - @assert length(bc_value) == length(components) - - for i in 1:length(components) - # find the global dof - globaldof = cc.dofs[local_facet_dofs[r[counter]]] - counter += 1 - - dbc_index = dofmapping[globaldof] - # Only DBC dofs are currently update!-able so don't modify inhomogeneities - # for affine constraints - if dofcoefficients[dbc_index] === nothing - inhomogeneities[dbc_index] = bc_value[i] - @debug println("prescribing value $(bc_value[i]) on global dof $(globaldof)") + for sdh in dh.subdofhandlers + facets_sdh_ = filter(x -> x[1] ∈ sdh.cellset, boundary_entities) + facets_sdh = OrderedSet{FacetIndex}(FacetIndex(x[1], x[2]) for x in facets_sdh_) + for fc in FacetIterator(sdh, facets_sdh, UpdateFlags(; nodes = false, coords = true, dofs = true)) + entityidx = fc.current_facet_id + cc = fc.cc + # no need to reinit!, enough to update current_entity since we only need geometric shape functions M + boundaryvalues.current_entity = entityidx + # local dof-range for this facet + r = local_facet_dofs_offset[entityidx]:(local_facet_dofs_offset[entityidx + 1] - 1) + counter = 1 + for location in 1:getnquadpoints(boundaryvalues) + x = spatial_coordinate(boundaryvalues, location, cc.coords) + bc_value = f(x, time) + @assert length(bc_value) == length(components) + + for i in 1:length(components) + # find the global dof + globaldof = cc.dofs[local_facet_dofs[r[counter]]] + counter += 1 + + dbc_index = dofmapping[globaldof] + # Only DBC dofs are currently update!-able so don't modify inhomogeneities + # for affine constraints + if dofcoefficients[dbc_index] === nothing + inhomogeneities[dbc_index] = bc_value[i] + @debug println("prescribing value $(bc_value[i]) on global dof $(globaldof)") + end end end end diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 0ac5da3faf..839114cf50 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/test/test_grid_dofhandler_vtk.jl @@ -253,7 +253,7 @@ end sdh1 = SubDofHandler(dh, Set([1])); add!(sdh1, :u, ip1) sdh2 = SubDofHandler(dh, Set([2])); add!(sdh2, :u, ip2) close!(dh) - ic = InterfaceCache(dh) + ic = InterfaceCache(sdh1, sdh2) reinit!(ic, InterfaceIndex(1, 2, 2, 3)) @test interfacedofs(ic) == collect(1:7) # Unit test of some utilities diff --git a/test/test_mixeddofhandler.jl b/test/test_mixeddofhandler.jl index ec4a7206b8..18eb2127aa 100644 --- a/test/test_mixeddofhandler.jl +++ b/test/test_mixeddofhandler.jl @@ -714,13 +714,19 @@ function test_celliterator_on_true_subdomain_smoketest() for cell in CellIterator(dh, [1, 2, 3]) end - - for cell in CellIterator(dh) - if cellid(cell) <= 3 - @test length(celldofs(cell)) == getnbasefunctions(ip) - else - @test length(celldofs(cell)) == 0 - end + # I smell smoke here + # for cell in CellIterator(dh) + # if cellid(cell) <= 3 + # @test length(celldofs(cell)) == getnbasefunctions(ip) + # else + # @test length(celldofs(cell)) == 0 + # end + # end + for cell in CellIterator(sdh) + @test length(celldofs(cell)) == getnbasefunctions(ip) + end + for cell in CellIterator(grid) + @test length(celldofs(cell)) == 0 end return end From 4a22deae6c720ea957cef66c30020873925342ea Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 13 Nov 2024 16:06:54 +0100 Subject: [PATCH 18/28] add cell type to SubDofHandler as typeparameter --- src/Dofs/DofHandler.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index 84f59e62b7..3a8e7572e3 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -12,7 +12,7 @@ Access some grid representation for the dof handler. """ get_grid(dh::AbstractDofHandler) -mutable struct SubDofHandler{DH} <: AbstractDofHandler +mutable struct SubDofHandler{DH, CT} <: AbstractDofHandler # From constructor const dh::DH const cellset::OrderedSet{Int} @@ -68,7 +68,7 @@ function SubDofHandler(dh::DH, cellset::AbstractVecOrSet{Int}) where {DH <: Abst end end # Construct and insert into the parent dh - sdh = SubDofHandler{typeof(dh)}(dh, convert_to_orderedset(cellset), Symbol[], Interpolation[], Int[], -1) + sdh = SubDofHandler{typeof(dh), CT}(dh, convert_to_orderedset(cellset), Symbol[], Interpolation[], Int[], -1) push!(dh.subdofhandlers, sdh) return sdh end From cba4e1d2702223394e1fdeb1236654cc9e4a5f77 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 13 Nov 2024 16:07:56 +0100 Subject: [PATCH 19/28] don't enoforce no-allocs for iterators using DofHandler --- src/iterators.jl | 12 +++++++----- test/test_iterators.jl | 14 +++++--------- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index 8cb7a096e2..bbbb58ba3a 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -67,11 +67,11 @@ function CellCache(sdh::SubDofHandler{<:DofHandler{dim}}, flags::UpdateFlags = U nodes = zeros(Int, N) coords = zeros(Vec{dim, get_coordinate_eltype(get_grid(sdh.dh))}, N) celldofs = zeros(Int, n) - return CellCache(flags, sdh.dh.grid, -1, nodes, coords, sdh.dh, celldofs) + return CellCache(flags, sdh.dh.grid, -1, nodes, coords, sdh, celldofs) end # TODO: Find a better way to make AllocCheck.jl happy -function reinit!(cc::CellCache{<:Any, <:Any, <:Nothing}, i::Int) +function reinit!(cc::CellCache, i::Int) cc.cellid = i if cc.flags.nodes resize!(cc.nodes, nnodes_per_cell(cc.grid, i)) @@ -88,14 +88,16 @@ function reinit!(cc::CellCache{<:Any, <:Any, <:Nothing}, i::Int) return cc end -function reinit!(cc::CellCache{<:Any, <:Any, <:AbstractDofHandler}, i::Int) +function reinit!(cc::CellCache{<:Any, <:AbstractGrid, <:SubDofHandler{<:Any, CT}}, i::Int) where {CT} # If we have a DofHandler the cells must be of the same type -> no need to resize cc.cellid = i if cc.flags.nodes - cellnodes!(cc.nodes, cc.grid, i) + cell = getcells(cc.grid, i)::CT + _cellnodes!(cc.nodes, cell) end if cc.flags.coords - getcoordinates!(cc.coords, cc.grid, i) + cell = getcells(cc.grid, i)::CT + getcoordinates!(cc.coords, cc.grid, cell) end if cc.dh !== nothing && cc.flags.dofs celldofs!(cc.dofs, cc.dh, i) diff --git a/test/test_iterators.jl b/test/test_iterators.jl index 5a8ca98891..d68d51787b 100644 --- a/test/test_iterators.jl +++ b/test/test_iterators.jl @@ -21,8 +21,9 @@ @test typeof(ii_dh) == typeof(ii_dh_top) @test typeof(ii_grid) == typeof(ii_grid_top) # Test for allocations - @test (testallocs_iterate(ii_dh_top)) === nothing - @test (testallocs_iterate(ii_dh)) === nothing + # TODO: find a way to avoid resize for single sdh + # @test (testallocs_iterate(ii_dh_top)) === nothing + # @test (testallocs_iterate(ii_dh)) === nothing # Iterators over grid can allocate due to potential resize! @test_throws AllocCheckFailure testallocs_iterate(ii_grid) @test_throws AllocCheckFailure testallocs_iterate(ii_grid_top) @@ -82,13 +83,8 @@ # Test that topology has no effect on iterator type @test typeof(ii_sdh) == typeof(ii_sdh_top) # Test for allocations - # try - # testallocs_iterate(ii_sdh_top) - # catch err - # @info err.errors[1] - # end - # @test (testallocs_iterate(ii_sdh_top)) === nothing - # @test (testallocs_iterate(ii_sdh)) === nothing + @test (testallocs_iterate(ii_sdh_top)) === nothing + @test (testallocs_iterate(ii_sdh)) === nothing end end From 5c6495e1bf8008203642bebd9aa090e7b31f9810 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 13 Nov 2024 16:23:04 +0100 Subject: [PATCH 20/28] undo constrainthandler voodoo --- src/Dofs/ConstraintHandler.jl | 90 ++++++++++++++++------------------- 1 file changed, 41 insertions(+), 49 deletions(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 81c1665499..a48cb3940f 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -311,15 +311,12 @@ function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfacets::AbstractVecOrSet # loop over all the faces in the set and add the global dofs to `constrained_dofs` constrained_dofs = Int[] - for sdh in ch.dh.subdofhandlers - facets_sdh_ = filter(x -> x[1] ∈ sdh.cellset, bcfacets) - facets_sdh = OrderedSet{FacetIndex}(FacetIndex(x[1], x[2]) for x in facets_sdh_) - for cc in FacetIterator(sdh, facets_sdh, UpdateFlags(; nodes = false, coords = false, dofs = true)) - facetidx = cc.current_facet_id - r = local_facet_dofs_offset[facetidx]:(local_facet_dofs_offset[facetidx + 1] - 1) - append!(constrained_dofs, cc.dofs[local_facet_dofs[r]]) # TODO: for-loop over r and simply push! to ch.prescribed_dofs - @debug println("adding dofs $(cc.dofs[local_facet_dofs[r]]) to dbc") - end + cc = CellCache(ch.dh, UpdateFlags(; nodes = false, coords = false, dofs = true)) + for (cellidx, facetidx) in bcfacets + reinit!(cc, cellidx) + r = local_facet_dofs_offset[facetidx]:(local_facet_dofs_offset[facetidx + 1] - 1) + append!(constrained_dofs, cc.dofs[local_facet_dofs[r]]) # TODO: for-loop over r and simply push! to ch.prescribed_dofs + @debug println("adding dofs $(cc.dofs[local_facet_dofs[r]]) to dbc") end # save it to the ConstraintHandler @@ -359,19 +356,16 @@ function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::AbstractVecOrSet{ interpol_points = getnbasefunctions(interpolation) node_dofs = zeros(Int, ncomps, nnodes) visited = falses(nnodes) - for sdh in ch.dh.subdofhandlers - for cell in CellIterator(sdh) # only go over cells that belong to current SubDofHandler - cellid(cell) ∈ cellset || continue - for idx in 1:min(interpol_points, length(cell.nodes)) - node = cell.nodes[idx] - if !visited[node] - noderange = (offset + (idx - 1) * field_dim + 1):(offset + idx * field_dim) # the dofs in this node - for (i, c) in enumerate(dbc.components) - node_dofs[i, node] = cell.dofs[noderange[c]] - @debug println("adding dof $(cell.dofs[noderange[c]]) to node_dofs") - end - visited[node] = true + for cell in CellIterator(ch.dh, cellset) # only go over cells that belong to current SubDofHandler + for idx in 1:min(interpol_points, length(cell.nodes)) + node = cell.nodes[idx] + if !visited[node] + noderange = (offset + (idx - 1) * field_dim + 1):(offset + idx * field_dim) # the dofs in this node + for (i, c) in enumerate(dbc.components) + node_dofs[i, node] = cell.dofs[noderange[c]] + @debug println("adding dof $(cell.dofs[noderange[c]]) to node_dofs") end + visited[node] = true end end end @@ -451,34 +445,32 @@ function _update!( dofmapping::Dict{Int, Int}, dofcoefficients::Vector{Union{Nothing, DofCoefficients{T}}}, time::Real ) where {T} - for sdh in dh.subdofhandlers - facets_sdh_ = filter(x -> x[1] ∈ sdh.cellset, boundary_entities) - facets_sdh = OrderedSet{FacetIndex}(FacetIndex(x[1], x[2]) for x in facets_sdh_) - for fc in FacetIterator(sdh, facets_sdh, UpdateFlags(; nodes = false, coords = true, dofs = true)) - entityidx = fc.current_facet_id - cc = fc.cc - # no need to reinit!, enough to update current_entity since we only need geometric shape functions M - boundaryvalues.current_entity = entityidx - # local dof-range for this facet - r = local_facet_dofs_offset[entityidx]:(local_facet_dofs_offset[entityidx + 1] - 1) - counter = 1 - for location in 1:getnquadpoints(boundaryvalues) - x = spatial_coordinate(boundaryvalues, location, cc.coords) - bc_value = f(x, time) - @assert length(bc_value) == length(components) - - for i in 1:length(components) - # find the global dof - globaldof = cc.dofs[local_facet_dofs[r[counter]]] - counter += 1 - - dbc_index = dofmapping[globaldof] - # Only DBC dofs are currently update!-able so don't modify inhomogeneities - # for affine constraints - if dofcoefficients[dbc_index] === nothing - inhomogeneities[dbc_index] = bc_value[i] - @debug println("prescribing value $(bc_value[i]) on global dof $(globaldof)") - end + cc = CellCache(dh, UpdateFlags(; nodes = false, coords = true, dofs = true)) + for (cellidx, entityidx) in boundary_entities + reinit!(cc, cellidx) + + # no need to reinit!, enough to update current_entity since we only need geometric shape functions M + boundaryvalues.current_entity = entityidx + + # local dof-range for this facet + r = local_facet_dofs_offset[entityidx]:(local_facet_dofs_offset[entityidx + 1] - 1) + counter = 1 + for location in 1:getnquadpoints(boundaryvalues) + x = spatial_coordinate(boundaryvalues, location, cc.coords) + bc_value = f(x, time) + @assert length(bc_value) == length(components) + + for i in 1:length(components) + # find the global dof + globaldof = cc.dofs[local_facet_dofs[r[counter]]] + counter += 1 + + dbc_index = dofmapping[globaldof] + # Only DBC dofs are currently update!-able so don't modify inhomogeneities + # for affine constraints + if dofcoefficients[dbc_index] === nothing + inhomogeneities[dbc_index] = bc_value[i] + @debug println("prescribing value $(bc_value[i]) on global dof $(globaldof)") end end end From 8f281d7775f39a468397a0004d8b31628ed2263a Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 13 Nov 2024 17:31:11 +0100 Subject: [PATCH 21/28] tests for CellIterator --- src/iterators.jl | 5 +- test/test_iterators.jl | 205 ++++++++++++++++++++++++++++++++++++++++- test/test_utils.jl | 8 ++ 3 files changed, 213 insertions(+), 5 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index bbbb58ba3a..8f1ca858df 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -288,9 +288,12 @@ function CellIterator(gridordh::Union{Grid, DofHandler}, flags::UpdateFlags) return CellIterator(gridordh, nothing, flags) end function CellIterator(sdh::SubDofHandler, flags::UpdateFlags = UpdateFlags()) - return CellIterator(CellCache(sdh, flags), sdh.cellset) + return CellIterator(sdh, sdh.cellset, flags) end +function CellIterator(sdh::SubDofHandler, set::IntegerCollection, flags::UpdateFlags = UpdateFlags()) + return CellIterator(CellCache(sdh, flags), set) +end @inline _getset(ci::CellIterator) = ci.set @inline _getcache(ci::CellIterator) = ci.cc diff --git a/test/test_iterators.jl b/test/test_iterators.jl index d68d51787b..2f7368799b 100644 --- a/test/test_iterators.jl +++ b/test/test_iterators.jl @@ -1,9 +1,206 @@ -@testset "InterfaceIterator" begin - @check_allocs function testallocs_iterate(iterator) - for _ in iterator end - return nothing +@testset "CellIterator" begin + @testset "Single domain" begin + grid = generate_grid(Hexahedron, (2, 2, 2)) + ip = Lagrange{RefHexahedron, 1}() + dh = DofHandler(grid) + add!(dh, :u, ip) + close!(dh) + + ci_dh = CellIterator(dh) + ci_grid = CellIterator(grid) + @test ci_dh.set == ci_grid.set + # Test for allocations + # TODO: find a way to avoid resize for single sdh + # @test (testallocs_iterate(ci_dh)) === nothing + # Iterators over grid can allocate due to potential resize! + @test_throws AllocCheckFailure testallocs_iterate(ci_grid) end + @testset "subdomains" begin + @testset "same cell type" begin + grid = generate_grid(Hexahedron, (2, 2, 2)) + ip = Lagrange{RefHexahedron, 1}() + dh = DofHandler(grid) + sdh1 = SubDofHandler(dh, Set(collect(1:4))) + sdh2 = SubDofHandler(dh, Set(collect(5:8))) + add!(sdh1, :u, ip^3) + add!(sdh2, :p, ip) + close!(dh) + + ci_sdh1 = CellIterator(sdh1, Set(collect(1:4))) + ci_sdh2 = CellIterator(sdh2, Set(collect(5:8))) + ci_dh = CellIterator(dh) + + # Test for allocations + @test (testallocs_iterate(ci_sdh1)) === nothing + @test (testallocs_iterate(ci_sdh2)) === nothing + @test_throws AllocCheckFailure testallocs_iterate(ci_dh) + end + + @testset "mixed cell types" begin + nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] + cells = [ + Quadrilateral((1, 2, 5, 4)), + Triangle((3, 5, 2)), + ] + grid = Grid(cells, nodes) + ip_quad = Lagrange{RefQuadrilateral, 1}() + ip_tri = Lagrange{RefTriangle, 1}() + dh = DofHandler(grid) + sdh_quad = SubDofHandler(dh, Set([1])) + sdh_tri = SubDofHandler(dh, Set([2])) + add!(sdh_quad, :u, ip_quad) + add!(sdh_tri, :u, ip_tri) + close!(dh) + + ci_quad = CellIterator(sdh_quad, Set([1])) + ci_tri = CellIterator(sdh_tri, Set([2])) + ci_dh = CellIterator(dh, Set([1])) + + # Test for allocations + try + testallocs_iterate(ci_quad) + + catch err + @info err.errors[3] + end + @test (testallocs_iterate(ci_quad)) === nothing + @test (testallocs_iterate(ci_tri)) === nothing + @test_throws AllocCheckFailure testallocs_iterate(ci_dh) + end + end + + @testset "error paths" begin + nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] + cells = [ + Quadrilateral((1, 2, 5, 4)), + Triangle((3, 5, 2)), + ] + grid = Grid(cells, nodes) + ip_quad = Lagrange{RefQuadrilateral, 1}() + ip_tri = Lagrange{RefTriangle, 1}() + dh = DofHandler(grid) + sdh_quad = SubDofHandler(dh, Set([1])) + sdh_tri = SubDofHandler(dh, Set([2])) + add!(sdh_quad, :u, ip_quad) + add!(sdh_tri, :u, ip_tri) + close!(dh) + @test_throws ErrorException("The cells in the cellset are not all of the same celltype.") CellIterator(dh) + end +end + +# @testset "FacetIterator" begin +# @testset "Single domain" begin +# grid = generate_grid(Hexahedron, (1000, 2, 1)) +# ip = Lagrange{RefHexahedron, 1}() +# dh = DofHandler(grid) +# add!(dh, :u, ip) +# close!(dh) +# topology = ExclusiveTopology(grid) + +# ii_dh_top = InterfaceIterator(dh, topology) +# ii_dh = InterfaceIterator(dh) +# ii_grid_top = InterfaceIterator(grid, topology) +# ii_grid = InterfaceIterator(grid) +# @test ii_dh.set == ii_dh_top.set == ii_grid.set == ii_grid_top.set +# # Test that topology has no effect on iterator type +# @test typeof(ii_dh) == typeof(ii_dh_top) +# @test typeof(ii_grid) == typeof(ii_grid_top) +# # Test for allocations +# # TODO: find a way to avoid resize for single sdh +# # @test (testallocs_iterate(ii_dh_top)) === nothing +# # @test (testallocs_iterate(ii_dh)) === nothing +# # Iterators over grid can allocate due to potential resize! +# @test_throws AllocCheckFailure testallocs_iterate(ii_grid) +# @test_throws AllocCheckFailure testallocs_iterate(ii_grid_top) +# end + +# @testset "subdomains" begin +# @testset "same cell type" begin +# grid = generate_grid(Hexahedron, (1000, 2, 1)) +# ip = Lagrange{RefHexahedron, 1}() +# dh = DofHandler(grid) +# sdh1 = SubDofHandler(dh, Set([collect(1:500)..., collect(1501:2000)...])) +# sdh2 = SubDofHandler(dh, Set([collect(501:1500)...])) +# add!(sdh1, :u, ip) +# add!(sdh2, :u, ip) +# close!(dh) +# topology = ExclusiveTopology(grid) + +# ii_sdh_top = InterfaceIterator(sdh1, sdh2, topology) +# ii_sdh = InterfaceIterator(sdh1, sdh2) +# ii_same_sdh = InterfaceIterator(sdh1, sdh1) +# @test ii_sdh.set == ii_sdh_top.set +# @test length(ii_sdh.set) == 1002 +# @test count(interface_index -> interface_index[1] == interface_index[3] - 1000, ii_sdh.set) == 500 +# @test count(interface_index -> interface_index[1] == interface_index[3] + 1000, ii_sdh.set) == 500 +# ii_sdh_flipped = InterfaceIterator(sdh2, sdh1) +# @test all(interface -> InterfaceIndex(interface[3], interface[4], interface[1], interface[2]) ∈ ii_sdh_flipped.set, ii_sdh.set) +# # Test that topology has no effect on iterator type +# @test typeof(ii_sdh) == typeof(ii_sdh_top) +# # Test for allocations +# @test (testallocs_iterate(ii_sdh)) === nothing +# @test (testallocs_iterate(ii_same_sdh)) === nothing +# end + +# @testset "mixed cell types" begin +# nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] +# cells = [ +# Quadrilateral((1, 2, 5, 4)), +# Triangle((3, 5, 2)), +# ] +# grid = Grid(cells, nodes) +# topology = ExclusiveTopology(grid) +# ip_quad = Lagrange{RefQuadrilateral, 1}() +# ip_tri = Lagrange{RefTriangle, 1}() +# dh = DofHandler(grid) +# sdh_quad = SubDofHandler(dh, Set([1])) +# sdh_tri = SubDofHandler(dh, Set([2])) +# add!(sdh_quad, :u, ip_quad) +# add!(sdh_tri, :u, ip_tri) +# close!(dh) + +# ii_sdh_top = InterfaceIterator(sdh_quad, sdh_tri, topology) +# ii_sdh = InterfaceIterator(sdh_quad, sdh_tri) +# @test ii_sdh.set == ii_sdh_top.set +# @test ii_sdh.set == Set([InterfaceIndex(1, 2, 2, 2)]) +# ii_sdh_flipped = InterfaceIterator(sdh_tri, sdh_quad) +# @test ii_sdh_flipped.set == Set([InterfaceIndex(2, 2, 1, 2)]) +# # Test that topology has no effect on iterator type +# @test typeof(ii_sdh) == typeof(ii_sdh_top) +# # Test for allocations +# @test (testallocs_iterate(ii_sdh_top)) === nothing +# @test (testallocs_iterate(ii_sdh)) === nothing +# end +# end + +# @testset "InterfaceIndex" begin +# # checkmate, codecov bot! +# @test isequal(InterfaceIndex(1, 2, 3, 4), (1, 2, 3, 4)) +# @test isequal((1, 2, 3, 4), InterfaceIndex(1, 2, 3, 4)) +# end + +# @testset "error paths" begin +# nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] +# cells = [ +# Quadrilateral((1, 2, 5, 4)), +# Triangle((3, 5, 2)), +# ] +# grid = Grid(cells, nodes) +# topology = ExclusiveTopology(grid) +# ip_quad = Lagrange{RefQuadrilateral, 1}() +# ip_tri = Lagrange{RefTriangle, 1}() +# dh = DofHandler(grid) +# sdh_quad = SubDofHandler(dh, Set([1])) +# sdh_tri = SubDofHandler(dh, Set([2])) +# add!(sdh_quad, :u, ip_quad) +# add!(sdh_tri, :u, ip_tri) +# close!(dh) +# @test_throws ErrorException("The cells in the set (set of InterfaceIndex) are not all of the same celltype on each side.") InterfaceIterator(sdh_quad, sdh_tri, Set([InterfaceIndex(2, 2, 1, 2)])) +# end +# end + +@testset "InterfaceIterator" begin @testset "Single domain" begin grid = generate_grid(Hexahedron, (1000, 2, 1)) ip = Lagrange{RefHexahedron, 1}() diff --git a/test/test_utils.jl b/test/test_utils.jl index 7f2696ffd8..a4016992c3 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -353,3 +353,11 @@ function function_value_from_physical_coord(interpolation::Interpolation, cell_c end return u end + +############################################################################### +# Function to test for possible allocations/dynamic dispatch when iterating # +############################################################################### +@check_allocs function testallocs_iterate(iterator) + for _ in iterator end + return nothing +end From 7b97feda8fb8c6a12826d8cdee5c141e9d854d93 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 13 Nov 2024 18:22:54 +0100 Subject: [PATCH 22/28] Construct FacetIterator without a set --- src/iterators.jl | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index 8f1ca858df..c966a158af 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -323,9 +323,9 @@ for faceindex in facetset # ... end """ -struct FacetIterator{FC <: FacetCache} +struct FacetIterator{FC <: FacetCache, SetType <: AbstractVecOrSet{FacetIndex}} fc::FC - set::OrderedSet{FacetIndex} + set::SetType end function FacetIterator( @@ -339,6 +339,25 @@ function FacetIterator( return FacetIterator(FacetCache(gridordh, flags), set) end +function FacetIterator( + gridordh::Union{<:AbstractGrid, <:DofHandler}, topology::ExclusiveTopology, + flags::UpdateFlags = UpdateFlags() + ) + grid = gridordh isa DofHandler ? get_grid(gridordh) : gridordh + set = Set(create_boundaryfacetset(grid, topology, _ -> true)) + return FacetIterator(gridordh, set, flags) +end + +function FacetIterator( + sdh::SubDofHandler, topology::ExclusiveTopology, + flags::UpdateFlags = UpdateFlags() + ) + grid = get_grid(sdh.dh) + set_unfiltered = create_boundaryfacetset(grid, topology, _ -> true) + set = Set(filter!(facet -> facet[1] ∈ sdh.cellset, set_unfiltered)) + return FacetIterator(sdh, set, flags) +end + @inline _getcache(fi::FacetIterator) = fi.fc @inline _getset(fi::FacetIterator) = fi.set From dcab007582b3d52a0d25ed06a5af54ed09b0cfc5 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 13 Nov 2024 18:23:08 +0100 Subject: [PATCH 23/28] Tests for FacetIterator --- test/test_iterators.jl | 187 ++++++++++++++++++----------------------- 1 file changed, 84 insertions(+), 103 deletions(-) diff --git a/test/test_iterators.jl b/test/test_iterators.jl index 2f7368799b..a389ff32bc 100644 --- a/test/test_iterators.jl +++ b/test/test_iterators.jl @@ -89,116 +89,97 @@ end end -# @testset "FacetIterator" begin -# @testset "Single domain" begin -# grid = generate_grid(Hexahedron, (1000, 2, 1)) -# ip = Lagrange{RefHexahedron, 1}() -# dh = DofHandler(grid) -# add!(dh, :u, ip) -# close!(dh) -# topology = ExclusiveTopology(grid) +@testset "FacetIterator" begin + @testset "Single domain" begin + grid = generate_grid(Hexahedron, (2, 2, 2)) + ip = Lagrange{RefHexahedron, 1}() + dh = DofHandler(grid) + add!(dh, :u, ip) + close!(dh) + topology = ExclusiveTopology(grid) + fi_dh = FacetIterator(dh, topology) + fi_grid = FacetIterator(grid, topology) + @test fi_dh.set == fi_grid.set -# ii_dh_top = InterfaceIterator(dh, topology) -# ii_dh = InterfaceIterator(dh) -# ii_grid_top = InterfaceIterator(grid, topology) -# ii_grid = InterfaceIterator(grid) -# @test ii_dh.set == ii_dh_top.set == ii_grid.set == ii_grid_top.set -# # Test that topology has no effect on iterator type -# @test typeof(ii_dh) == typeof(ii_dh_top) -# @test typeof(ii_grid) == typeof(ii_grid_top) -# # Test for allocations -# # TODO: find a way to avoid resize for single sdh -# # @test (testallocs_iterate(ii_dh_top)) === nothing -# # @test (testallocs_iterate(ii_dh)) === nothing -# # Iterators over grid can allocate due to potential resize! -# @test_throws AllocCheckFailure testallocs_iterate(ii_grid) -# @test_throws AllocCheckFailure testallocs_iterate(ii_grid_top) -# end + # Test for allocations + # TODO: find a way to avoid resize for single sdh + # @test (testallocs_iterate(fi_dh)) === nothing + # Iterators over grid can allocate due to potential resize! + @test_throws AllocCheckFailure testallocs_iterate(fi_grid) + end -# @testset "subdomains" begin -# @testset "same cell type" begin -# grid = generate_grid(Hexahedron, (1000, 2, 1)) -# ip = Lagrange{RefHexahedron, 1}() -# dh = DofHandler(grid) -# sdh1 = SubDofHandler(dh, Set([collect(1:500)..., collect(1501:2000)...])) -# sdh2 = SubDofHandler(dh, Set([collect(501:1500)...])) -# add!(sdh1, :u, ip) -# add!(sdh2, :u, ip) -# close!(dh) -# topology = ExclusiveTopology(grid) + @testset "subdomains" begin + @testset "same cell type" begin + grid = generate_grid(Hexahedron, (2, 2, 2)) + ip = Lagrange{RefHexahedron, 1}() + dh = DofHandler(grid) + sdh1 = SubDofHandler(dh, Set(collect(1:4))) + sdh2 = SubDofHandler(dh, Set(collect(5:8))) + add!(sdh1, :u, ip^3) + add!(sdh2, :p, ip) + close!(dh) + topology = ExclusiveTopology(grid) -# ii_sdh_top = InterfaceIterator(sdh1, sdh2, topology) -# ii_sdh = InterfaceIterator(sdh1, sdh2) -# ii_same_sdh = InterfaceIterator(sdh1, sdh1) -# @test ii_sdh.set == ii_sdh_top.set -# @test length(ii_sdh.set) == 1002 -# @test count(interface_index -> interface_index[1] == interface_index[3] - 1000, ii_sdh.set) == 500 -# @test count(interface_index -> interface_index[1] == interface_index[3] + 1000, ii_sdh.set) == 500 -# ii_sdh_flipped = InterfaceIterator(sdh2, sdh1) -# @test all(interface -> InterfaceIndex(interface[3], interface[4], interface[1], interface[2]) ∈ ii_sdh_flipped.set, ii_sdh.set) -# # Test that topology has no effect on iterator type -# @test typeof(ii_sdh) == typeof(ii_sdh_top) -# # Test for allocations -# @test (testallocs_iterate(ii_sdh)) === nothing -# @test (testallocs_iterate(ii_same_sdh)) === nothing -# end + fi_sdh1 = FacetIterator(sdh1, topology) + fi_sdh2 = FacetIterator(sdh2, topology) + @test isempty(fi_sdh1.set ∩ fi_sdh2.set) + @test length(fi_sdh1.set) == 12 + # Test for allocations + @test (testallocs_iterate(fi_sdh1)) === nothing + @test (testallocs_iterate(fi_sdh2)) === nothing + end -# @testset "mixed cell types" begin -# nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] -# cells = [ -# Quadrilateral((1, 2, 5, 4)), -# Triangle((3, 5, 2)), -# ] -# grid = Grid(cells, nodes) -# topology = ExclusiveTopology(grid) -# ip_quad = Lagrange{RefQuadrilateral, 1}() -# ip_tri = Lagrange{RefTriangle, 1}() -# dh = DofHandler(grid) -# sdh_quad = SubDofHandler(dh, Set([1])) -# sdh_tri = SubDofHandler(dh, Set([2])) -# add!(sdh_quad, :u, ip_quad) -# add!(sdh_tri, :u, ip_tri) -# close!(dh) + @testset "mixed cell types" begin + nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] + cells = [ + Quadrilateral((1, 2, 5, 4)), + Triangle((3, 5, 2)), + ] + grid = Grid(cells, nodes) + topology = ExclusiveTopology(grid) + ip_quad = Lagrange{RefQuadrilateral, 1}() + ip_tri = Lagrange{RefTriangle, 1}() + dh = DofHandler(grid) + sdh_quad = SubDofHandler(dh, Set([1])) + sdh_tri = SubDofHandler(dh, Set([2])) + add!(sdh_quad, :u, ip_quad) + add!(sdh_tri, :u, ip_tri) + close!(dh) -# ii_sdh_top = InterfaceIterator(sdh_quad, sdh_tri, topology) -# ii_sdh = InterfaceIterator(sdh_quad, sdh_tri) -# @test ii_sdh.set == ii_sdh_top.set -# @test ii_sdh.set == Set([InterfaceIndex(1, 2, 2, 2)]) -# ii_sdh_flipped = InterfaceIterator(sdh_tri, sdh_quad) -# @test ii_sdh_flipped.set == Set([InterfaceIndex(2, 2, 1, 2)]) -# # Test that topology has no effect on iterator type -# @test typeof(ii_sdh) == typeof(ii_sdh_top) -# # Test for allocations -# @test (testallocs_iterate(ii_sdh_top)) === nothing -# @test (testallocs_iterate(ii_sdh)) === nothing -# end -# end + fi_quad = FacetIterator(sdh_quad, topology) + fi_tri = FacetIterator(sdh_tri, topology) + @test isempty(fi_quad.set ∩ fi_tri.set) + # Test for allocations + @test (testallocs_iterate(fi_quad)) === nothing + @test (testallocs_iterate(fi_tri)) === nothing + end + end -# @testset "InterfaceIndex" begin -# # checkmate, codecov bot! -# @test isequal(InterfaceIndex(1, 2, 3, 4), (1, 2, 3, 4)) -# @test isequal((1, 2, 3, 4), InterfaceIndex(1, 2, 3, 4)) -# end + @testset "FacetIndex" begin + # checkmate, codecov bot! + @test isequal(FacetIndex(1, 2), (1, 2)) + @test isequal((1, 2), FacetIndex(1, 2)) + end -# @testset "error paths" begin -# nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] -# cells = [ -# Quadrilateral((1, 2, 5, 4)), -# Triangle((3, 5, 2)), -# ] -# grid = Grid(cells, nodes) -# topology = ExclusiveTopology(grid) -# ip_quad = Lagrange{RefQuadrilateral, 1}() -# ip_tri = Lagrange{RefTriangle, 1}() -# dh = DofHandler(grid) -# sdh_quad = SubDofHandler(dh, Set([1])) -# sdh_tri = SubDofHandler(dh, Set([2])) -# add!(sdh_quad, :u, ip_quad) -# add!(sdh_tri, :u, ip_tri) -# close!(dh) -# @test_throws ErrorException("The cells in the set (set of InterfaceIndex) are not all of the same celltype on each side.") InterfaceIterator(sdh_quad, sdh_tri, Set([InterfaceIndex(2, 2, 1, 2)])) -# end -# end + @testset "error paths" begin + nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] + cells = [ + Quadrilateral((1, 2, 5, 4)), + Triangle((3, 5, 2)), + ] + grid = Grid(cells, nodes) + topology = ExclusiveTopology(grid) + ip_quad = Lagrange{RefQuadrilateral, 1}() + ip_tri = Lagrange{RefTriangle, 1}() + dh = DofHandler(grid) + sdh_quad = SubDofHandler(dh, Set([1])) + sdh_tri = SubDofHandler(dh, Set([2])) + add!(sdh_quad, :u, ip_quad) + add!(sdh_tri, :u, ip_tri) + close!(dh) + @test_throws ErrorException("The cells in the set (set of FacetIndex) are not all of the same celltype.") FacetIterator(dh, topology) + end +end @testset "InterfaceIterator" begin @testset "Single domain" begin From a38ed022dbc7914b344920698ddab7debf7f1e57 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 13 Nov 2024 18:39:15 +0100 Subject: [PATCH 24/28] Fix bug where boundary sets are broken for mixed grids --- src/Grid/utils.jl | 1 + test/test_grid_addboundaryset.jl | 17 +++++++++++++++++ 2 files changed, 18 insertions(+) diff --git a/src/Grid/utils.jl b/src/Grid/utils.jl index ce5937010e..04400ae765 100644 --- a/src/Grid/utils.jl +++ b/src/Grid/utils.jl @@ -165,6 +165,7 @@ function _create_boundaryset(f::Function, grid::AbstractGrid, top::ExclusiveTopo cell_idx = ff_nh_idx[1] facet_nr = ff_nh_idx[2] cell = getcells(grid, cell_idx) + length(facets(cell)) < facet_nr && continue facet_nodes = facets(cell)[facet_nr] for (subentity_idx, subentity_nodes) in pairs(boundaryfunction(BI)(cell)) if Base.all(n -> n in facet_nodes, subentity_nodes) diff --git a/test/test_grid_addboundaryset.jl b/test/test_grid_addboundaryset.jl index ed3ca549ad..ed85602f50 100644 --- a/test/test_grid_addboundaryset.jl +++ b/test/test_grid_addboundaryset.jl @@ -200,4 +200,21 @@ addboundaryfacetset!(grid, topology, "test_boundary_facetset", filter_function) @test getfacetset(grid, "test_boundary_facetset") == Ferrite.create_boundaryfacetset(grid, topology, filter_function) end + + @testset "addboundaryset Mixed grid" begin + nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] + cells = [ + Quadrilateral((1, 2, 5, 4)), + Triangle((3, 5, 2)), + ] + grid = Grid(cells, nodes) + topology = ExclusiveTopology(grid) + @test extractboundary(grid, topology) == extractboundarycheck(grid) + + filter_function(x) = x[1] > 0 + addboundaryvertexset!(grid, topology, "test_boundary_vertexset", filter_function) + @test getvertexset(grid, "test_boundary_vertexset") == Ferrite.create_boundaryvertexset(grid, topology, filter_function) + addboundaryfacetset!(grid, topology, "test_boundary_facetset", filter_function) + @test getfacetset(grid, "test_boundary_facetset") == Ferrite.create_boundaryfacetset(grid, topology, filter_function) + end end From 94948c2bac6d093029df9bc5c1f6f92af3b0f901 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 13 Nov 2024 19:12:26 +0100 Subject: [PATCH 25/28] move AllocCheck to test deps --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 738ead2a9b..f03c4aa344 100644 --- a/Project.toml +++ b/Project.toml @@ -3,7 +3,6 @@ uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992" version = "1.0.0" [deps] -AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a" EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -55,6 +54,7 @@ SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" TaskLocalValues = "ed4db957-447d-4319-bfb6-7fa9ae7ecf34" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a" [targets] -test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "OhMyThreads", "ProgressMeter", "Random", "SHA", "TaskLocalValues", "Test", "TimerOutputs", "Logging"] +test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "OhMyThreads", "ProgressMeter", "Random", "SHA", "TaskLocalValues", "Test", "TimerOutputs", "Logging", "AllocCheck"] From be9009a8d156ded273c85223a2a306a5dd69af18 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 13 Nov 2024 19:27:42 +0100 Subject: [PATCH 26/28] please --- test/test_interfacevalues.jl | 36 +++++++++++++++++++++++------------- test/test_utils.jl | 8 -------- 2 files changed, 23 insertions(+), 21 deletions(-) diff --git a/test/test_interfacevalues.jl b/test/test_interfacevalues.jl index 8fb31dc741..3254014a58 100644 --- a/test/test_interfacevalues.jl +++ b/test/test_interfacevalues.jl @@ -1,3 +1,11 @@ +############################################################################### +# Function to test for possible allocations/dynamic dispatch when iterating # +############################################################################### +@check_allocs function testallocs_iterate(iterator) + for _ in iterator end + return nothing +end + @testset "InterfaceValues" begin function test_interfacevalues(grid::Ferrite.AbstractGrid, iv::InterfaceValues; tol = 0) ip_here = Ferrite.function_interpolation(iv.here) @@ -185,20 +193,22 @@ test_interfacevalues(grid, iv; tol = 5 * eps(Float64)) end end - # @testset "Mixed elements 2D grids" begin # TODO: this shouldn't work because it should change the FacetValues object - # dim = 2 - # nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] - # cells = [ - # Quadrilateral((1,2,5,4)), - # Triangle((3,5,2)), - # ] + @testset "Mixed elements 2D grids" begin + dim = 2 + nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] + cells = [ + Quadrilateral((1, 2, 5, 4)), + Triangle((3, 5, 2)), + ] - # grid = Grid(cells, nodes) - # topology = ExclusiveTopology(grid) - # test_interfacevalues(grid, - # DiscontinuousLagrange{RefQuadrilateral, 1}(), FacetQuadratureRule{RefQuadrilateral}(2), - # DiscontinuousLagrange{RefTriangle, 1}(), FacetQuadratureRule{RefTriangle}(2)) - # end + grid = Grid(cells, nodes) + topology = ExclusiveTopology(grid) + test_interfacevalues( + grid, + DiscontinuousLagrange{RefQuadrilateral, 1}(), FacetQuadratureRule{RefQuadrilateral}(2), + DiscontinuousLagrange{RefTriangle, 1}(), FacetQuadratureRule{RefTriangle}(2) + ) + end @testset "Unordered nodes 3D" begin dim = 2 nodes = [ diff --git a/test/test_utils.jl b/test/test_utils.jl index a4016992c3..7f2696ffd8 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -353,11 +353,3 @@ function function_value_from_physical_coord(interpolation::Interpolation, cell_c end return u end - -############################################################################### -# Function to test for possible allocations/dynamic dispatch when iterating # -############################################################################### -@check_allocs function testallocs_iterate(iterator) - for _ in iterator end - return nothing -end From 0143341c5644921f4a863fbbef3ed96a20c09e12 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 13 Nov 2024 20:18:45 +0100 Subject: [PATCH 27/28] enable interfacevalues test with mixed grid --- test/test_interfacevalues.jl | 51 +++-- test/test_iterators.jl | 431 ++++++++++++++++++----------------- 2 files changed, 251 insertions(+), 231 deletions(-) diff --git a/test/test_interfacevalues.jl b/test/test_interfacevalues.jl index 3254014a58..6e0aced047 100644 --- a/test/test_interfacevalues.jl +++ b/test/test_interfacevalues.jl @@ -1,21 +1,12 @@ -############################################################################### -# Function to test for possible allocations/dynamic dispatch when iterating # -############################################################################### -@check_allocs function testallocs_iterate(iterator) - for _ in iterator end - return nothing -end - @testset "InterfaceValues" begin - function test_interfacevalues(grid::Ferrite.AbstractGrid, iv::InterfaceValues; tol = 0) + function test_interfacevalues(iv::InterfaceValues, args...; tol = 0) ip_here = Ferrite.function_interpolation(iv.here) ip_there = Ferrite.function_interpolation(iv.there) rdim = Ferrite.getrefdim(ip_here) n_basefuncs = getnbasefunctions(ip_here) + getnbasefunctions(ip_there) @test getnbasefunctions(iv) == n_basefuncs - - for ic in InterfaceIterator(grid) + for ic in InterfaceIterator(args...) reinit!(iv, ic) coords_here, coords_there = getcoordinates(ic) nqp = getnquadpoints(iv) @@ -160,7 +151,7 @@ end for func_interpol in (scalar_interpol, VectorizedInterpolation(scalar_interpol)) iv = cell_shape ∈ (QuadraticLine, QuadraticQuadrilateral, QuadraticTriangle, QuadraticTetrahedron) ? InterfaceValues(quad_rule, func_interpol, ip) : InterfaceValues(quad_rule, func_interpol) - test_interfacevalues(grid, iv) + test_interfacevalues(iv, grid) end end # Custom quadrature @@ -190,7 +181,7 @@ end end for func_interpol in (scalar_interpol, VectorizedInterpolation(scalar_interpol)) iv = InterfaceValues(quad_rule, func_interpol) - test_interfacevalues(grid, iv; tol = 5 * eps(Float64)) + test_interfacevalues(iv, grid; tol = 5 * eps(Float64)) end end @testset "Mixed elements 2D grids" begin @@ -203,11 +194,31 @@ end grid = Grid(cells, nodes) topology = ExclusiveTopology(grid) - test_interfacevalues( - grid, - DiscontinuousLagrange{RefQuadrilateral, 1}(), FacetQuadratureRule{RefQuadrilateral}(2), - DiscontinuousLagrange{RefTriangle, 1}(), FacetQuadratureRule{RefTriangle}(2) - ) + qr_facet_quad = FacetQuadratureRule{RefQuadrilateral}(2) + qr_facet_tri = FacetQuadratureRule{RefTriangle}(2) + ip_quad = DiscontinuousLagrange{RefQuadrilateral, 1}() + ip_tri = DiscontinuousLagrange{RefTriangle, 1}() + dh = DofHandler(grid) + sdh_quad = SubDofHandler(dh, Set([1])) + sdh_tri = SubDofHandler(dh, Set([2])) + add!(sdh_quad, :u, ip_quad) + add!(sdh_tri, :u, ip_tri) + qr_collection = (qr_facet_quad, qr_facet_tri) + ip_collection = (ip_quad, ip_tri) + close!(dh) + for (sdh_here_idx, sdh_here) in enumerate(dh.subdofhandlers) + for (sdh_there_idx, sdh_there) in enumerate(dh.subdofhandlers) + sdh_there_idx < sdh_here_idx && continue + iv = InterfaceValues( + qr_collection[sdh_here_idx], ip_collection[sdh_here_idx], + qr_collection[sdh_there_idx], ip_collection[sdh_there_idx] + ) + test_interfacevalues( + iv, sdh_here, sdh_there + ) + end + end + end @testset "Unordered nodes 3D" begin dim = 2 @@ -224,8 +235,8 @@ end grid = Grid(cells, nodes) test_interfacevalues( - grid, - InterfaceValues(FacetQuadratureRule{RefHexahedron}(2), DiscontinuousLagrange{RefHexahedron, 1}()) + InterfaceValues(FacetQuadratureRule{RefHexahedron}(2), DiscontinuousLagrange{RefHexahedron, 1}()), + grid ) end @testset "Interface dof_range" begin diff --git a/test/test_iterators.jl b/test/test_iterators.jl index a389ff32bc..cd011799b7 100644 --- a/test/test_iterators.jl +++ b/test/test_iterators.jl @@ -1,43 +1,84 @@ -@testset "CellIterator" begin - @testset "Single domain" begin - grid = generate_grid(Hexahedron, (2, 2, 2)) - ip = Lagrange{RefHexahedron, 1}() - dh = DofHandler(grid) - add!(dh, :u, ip) - close!(dh) - - ci_dh = CellIterator(dh) - ci_grid = CellIterator(grid) - @test ci_dh.set == ci_grid.set - # Test for allocations - # TODO: find a way to avoid resize for single sdh - # @test (testallocs_iterate(ci_dh)) === nothing - # Iterators over grid can allocate due to potential resize! - @test_throws AllocCheckFailure testallocs_iterate(ci_grid) +@testset "Iterators" begin + ############################################################################### + # Function to test for possible allocations/dynamic dispatch when iterating # + ############################################################################### + @check_allocs function testallocs_iterate(iterator) + for _ in iterator end + return nothing end - - @testset "subdomains" begin - @testset "same cell type" begin + @testset "CellIterator" begin + @testset "Single domain" begin grid = generate_grid(Hexahedron, (2, 2, 2)) ip = Lagrange{RefHexahedron, 1}() dh = DofHandler(grid) - sdh1 = SubDofHandler(dh, Set(collect(1:4))) - sdh2 = SubDofHandler(dh, Set(collect(5:8))) - add!(sdh1, :u, ip^3) - add!(sdh2, :p, ip) + add!(dh, :u, ip) close!(dh) - ci_sdh1 = CellIterator(sdh1, Set(collect(1:4))) - ci_sdh2 = CellIterator(sdh2, Set(collect(5:8))) ci_dh = CellIterator(dh) - + ci_grid = CellIterator(grid) + @test ci_dh.set == ci_grid.set # Test for allocations - @test (testallocs_iterate(ci_sdh1)) === nothing - @test (testallocs_iterate(ci_sdh2)) === nothing - @test_throws AllocCheckFailure testallocs_iterate(ci_dh) + # TODO: find a way to avoid resize for single sdh + # @test (testallocs_iterate(ci_dh)) === nothing + # Iterators over grid can allocate due to potential resize! + @test_throws AllocCheckFailure testallocs_iterate(ci_grid) end - @testset "mixed cell types" begin + @testset "subdomains" begin + @testset "same cell type" begin + grid = generate_grid(Hexahedron, (2, 2, 2)) + ip = Lagrange{RefHexahedron, 1}() + dh = DofHandler(grid) + sdh1 = SubDofHandler(dh, Set(collect(1:4))) + sdh2 = SubDofHandler(dh, Set(collect(5:8))) + add!(sdh1, :u, ip^3) + add!(sdh2, :p, ip) + close!(dh) + + ci_sdh1 = CellIterator(sdh1, Set(collect(1:4))) + ci_sdh2 = CellIterator(sdh2, Set(collect(5:8))) + ci_dh = CellIterator(dh) + + # Test for allocations + @test (testallocs_iterate(ci_sdh1)) === nothing + @test (testallocs_iterate(ci_sdh2)) === nothing + @test_throws AllocCheckFailure testallocs_iterate(ci_dh) + end + + @testset "mixed cell types" begin + nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] + cells = [ + Quadrilateral((1, 2, 5, 4)), + Triangle((3, 5, 2)), + ] + grid = Grid(cells, nodes) + ip_quad = Lagrange{RefQuadrilateral, 1}() + ip_tri = Lagrange{RefTriangle, 1}() + dh = DofHandler(grid) + sdh_quad = SubDofHandler(dh, Set([1])) + sdh_tri = SubDofHandler(dh, Set([2])) + add!(sdh_quad, :u, ip_quad) + add!(sdh_tri, :u, ip_tri) + close!(dh) + + ci_quad = CellIterator(sdh_quad, Set([1])) + ci_tri = CellIterator(sdh_tri, Set([2])) + ci_dh = CellIterator(dh, Set([1])) + + # Test for allocations + try + testallocs_iterate(ci_quad) + + catch err + @info err.errors[3] + end + @test (testallocs_iterate(ci_quad)) === nothing + @test (testallocs_iterate(ci_tri)) === nothing + @test_throws AllocCheckFailure testallocs_iterate(ci_dh) + end + end + + @testset "error paths" begin nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] cells = [ Quadrilateral((1, 2, 5, 4)), @@ -52,84 +93,83 @@ add!(sdh_quad, :u, ip_quad) add!(sdh_tri, :u, ip_tri) close!(dh) - - ci_quad = CellIterator(sdh_quad, Set([1])) - ci_tri = CellIterator(sdh_tri, Set([2])) - ci_dh = CellIterator(dh, Set([1])) - - # Test for allocations - try - testallocs_iterate(ci_quad) - - catch err - @info err.errors[3] - end - @test (testallocs_iterate(ci_quad)) === nothing - @test (testallocs_iterate(ci_tri)) === nothing - @test_throws AllocCheckFailure testallocs_iterate(ci_dh) + @test_throws ErrorException("The cells in the cellset are not all of the same celltype.") CellIterator(dh) end end - @testset "error paths" begin - nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] - cells = [ - Quadrilateral((1, 2, 5, 4)), - Triangle((3, 5, 2)), - ] - grid = Grid(cells, nodes) - ip_quad = Lagrange{RefQuadrilateral, 1}() - ip_tri = Lagrange{RefTriangle, 1}() - dh = DofHandler(grid) - sdh_quad = SubDofHandler(dh, Set([1])) - sdh_tri = SubDofHandler(dh, Set([2])) - add!(sdh_quad, :u, ip_quad) - add!(sdh_tri, :u, ip_tri) - close!(dh) - @test_throws ErrorException("The cells in the cellset are not all of the same celltype.") CellIterator(dh) - end -end - -@testset "FacetIterator" begin - @testset "Single domain" begin - grid = generate_grid(Hexahedron, (2, 2, 2)) - ip = Lagrange{RefHexahedron, 1}() - dh = DofHandler(grid) - add!(dh, :u, ip) - close!(dh) - topology = ExclusiveTopology(grid) - fi_dh = FacetIterator(dh, topology) - fi_grid = FacetIterator(grid, topology) - @test fi_dh.set == fi_grid.set - - # Test for allocations - # TODO: find a way to avoid resize for single sdh - # @test (testallocs_iterate(fi_dh)) === nothing - # Iterators over grid can allocate due to potential resize! - @test_throws AllocCheckFailure testallocs_iterate(fi_grid) - end - - @testset "subdomains" begin - @testset "same cell type" begin + @testset "FacetIterator" begin + @testset "Single domain" begin grid = generate_grid(Hexahedron, (2, 2, 2)) ip = Lagrange{RefHexahedron, 1}() dh = DofHandler(grid) - sdh1 = SubDofHandler(dh, Set(collect(1:4))) - sdh2 = SubDofHandler(dh, Set(collect(5:8))) - add!(sdh1, :u, ip^3) - add!(sdh2, :p, ip) + add!(dh, :u, ip) close!(dh) topology = ExclusiveTopology(grid) + fi_dh = FacetIterator(dh, topology) + fi_grid = FacetIterator(grid, topology) + @test fi_dh.set == fi_grid.set - fi_sdh1 = FacetIterator(sdh1, topology) - fi_sdh2 = FacetIterator(sdh2, topology) - @test isempty(fi_sdh1.set ∩ fi_sdh2.set) - @test length(fi_sdh1.set) == 12 # Test for allocations - @test (testallocs_iterate(fi_sdh1)) === nothing - @test (testallocs_iterate(fi_sdh2)) === nothing + # TODO: find a way to avoid resize for single sdh + # @test (testallocs_iterate(fi_dh)) === nothing + # Iterators over grid can allocate due to potential resize! + @test_throws AllocCheckFailure testallocs_iterate(fi_grid) end - @testset "mixed cell types" begin + @testset "subdomains" begin + @testset "same cell type" begin + grid = generate_grid(Hexahedron, (2, 2, 2)) + ip = Lagrange{RefHexahedron, 1}() + dh = DofHandler(grid) + sdh1 = SubDofHandler(dh, Set(collect(1:4))) + sdh2 = SubDofHandler(dh, Set(collect(5:8))) + add!(sdh1, :u, ip^3) + add!(sdh2, :p, ip) + close!(dh) + topology = ExclusiveTopology(grid) + + fi_sdh1 = FacetIterator(sdh1, topology) + fi_sdh2 = FacetIterator(sdh2, topology) + @test isempty(fi_sdh1.set ∩ fi_sdh2.set) + @test length(fi_sdh1.set) == 12 + # Test for allocations + @test (testallocs_iterate(fi_sdh1)) === nothing + @test (testallocs_iterate(fi_sdh2)) === nothing + end + + @testset "mixed cell types" begin + nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] + cells = [ + Quadrilateral((1, 2, 5, 4)), + Triangle((3, 5, 2)), + ] + grid = Grid(cells, nodes) + topology = ExclusiveTopology(grid) + ip_quad = Lagrange{RefQuadrilateral, 1}() + ip_tri = Lagrange{RefTriangle, 1}() + dh = DofHandler(grid) + sdh_quad = SubDofHandler(dh, Set([1])) + sdh_tri = SubDofHandler(dh, Set([2])) + add!(sdh_quad, :u, ip_quad) + add!(sdh_tri, :u, ip_tri) + close!(dh) + + fi_quad = FacetIterator(sdh_quad, topology) + fi_tri = FacetIterator(sdh_tri, topology) + @test isempty(fi_quad.set ∩ fi_tri.set) + # Test for allocations + @test (testallocs_iterate(fi_quad)) === nothing + @test (testallocs_iterate(fi_tri)) === nothing + end + end + + @testset "FacetIndex" begin + # checkmate, codecov bot! + @test isequal(FacetIndex(1, 2), (1, 2)) + @test isequal((1, 2), FacetIndex(1, 2)) + end + + @testset "error paths" begin nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] cells = [ Quadrilateral((1, 2, 5, 4)), @@ -145,97 +185,102 @@ end add!(sdh_quad, :u, ip_quad) add!(sdh_tri, :u, ip_tri) close!(dh) - - fi_quad = FacetIterator(sdh_quad, topology) - fi_tri = FacetIterator(sdh_tri, topology) - @test isempty(fi_quad.set ∩ fi_tri.set) - # Test for allocations - @test (testallocs_iterate(fi_quad)) === nothing - @test (testallocs_iterate(fi_tri)) === nothing + @test_throws ErrorException("The cells in the set (set of FacetIndex) are not all of the same celltype.") FacetIterator(dh, topology) end end - @testset "FacetIndex" begin - # checkmate, codecov bot! - @test isequal(FacetIndex(1, 2), (1, 2)) - @test isequal((1, 2), FacetIndex(1, 2)) - end - - @testset "error paths" begin - nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] - cells = [ - Quadrilateral((1, 2, 5, 4)), - Triangle((3, 5, 2)), - ] - grid = Grid(cells, nodes) - topology = ExclusiveTopology(grid) - ip_quad = Lagrange{RefQuadrilateral, 1}() - ip_tri = Lagrange{RefTriangle, 1}() - dh = DofHandler(grid) - sdh_quad = SubDofHandler(dh, Set([1])) - sdh_tri = SubDofHandler(dh, Set([2])) - add!(sdh_quad, :u, ip_quad) - add!(sdh_tri, :u, ip_tri) - close!(dh) - @test_throws ErrorException("The cells in the set (set of FacetIndex) are not all of the same celltype.") FacetIterator(dh, topology) - end -end - -@testset "InterfaceIterator" begin - @testset "Single domain" begin - grid = generate_grid(Hexahedron, (1000, 2, 1)) - ip = Lagrange{RefHexahedron, 1}() - dh = DofHandler(grid) - add!(dh, :u, ip) - close!(dh) - topology = ExclusiveTopology(grid) - - ii_dh_top = InterfaceIterator(dh, topology) - ii_dh = InterfaceIterator(dh) - ii_grid_top = InterfaceIterator(grid, topology) - ii_grid = InterfaceIterator(grid) - @test ii_dh.set == ii_dh_top.set == ii_grid.set == ii_grid_top.set - # Test that topology has no effect on iterator type - @test typeof(ii_dh) == typeof(ii_dh_top) - @test typeof(ii_grid) == typeof(ii_grid_top) - # Test for allocations - # TODO: find a way to avoid resize for single sdh - # @test (testallocs_iterate(ii_dh_top)) === nothing - # @test (testallocs_iterate(ii_dh)) === nothing - # Iterators over grid can allocate due to potential resize! - @test_throws AllocCheckFailure testallocs_iterate(ii_grid) - @test_throws AllocCheckFailure testallocs_iterate(ii_grid_top) - end - - @testset "subdomains" begin - @testset "same cell type" begin + @testset "InterfaceIterator" begin + @testset "Single domain" begin grid = generate_grid(Hexahedron, (1000, 2, 1)) ip = Lagrange{RefHexahedron, 1}() dh = DofHandler(grid) - sdh1 = SubDofHandler(dh, Set([collect(1:500)..., collect(1501:2000)...])) - sdh2 = SubDofHandler(dh, Set([collect(501:1500)...])) - add!(sdh1, :u, ip) - add!(sdh2, :u, ip) + add!(dh, :u, ip) close!(dh) topology = ExclusiveTopology(grid) - ii_sdh_top = InterfaceIterator(sdh1, sdh2, topology) - ii_sdh = InterfaceIterator(sdh1, sdh2) - ii_same_sdh = InterfaceIterator(sdh1, sdh1) - @test ii_sdh.set == ii_sdh_top.set - @test length(ii_sdh.set) == 1002 - @test count(interface_index -> interface_index[1] == interface_index[3] - 1000, ii_sdh.set) == 500 - @test count(interface_index -> interface_index[1] == interface_index[3] + 1000, ii_sdh.set) == 500 - ii_sdh_flipped = InterfaceIterator(sdh2, sdh1) - @test all(interface -> InterfaceIndex(interface[3], interface[4], interface[1], interface[2]) ∈ ii_sdh_flipped.set, ii_sdh.set) + ii_dh_top = InterfaceIterator(dh, topology) + ii_dh = InterfaceIterator(dh) + ii_grid_top = InterfaceIterator(grid, topology) + ii_grid = InterfaceIterator(grid) + @test ii_dh.set == ii_dh_top.set == ii_grid.set == ii_grid_top.set # Test that topology has no effect on iterator type - @test typeof(ii_sdh) == typeof(ii_sdh_top) + @test typeof(ii_dh) == typeof(ii_dh_top) + @test typeof(ii_grid) == typeof(ii_grid_top) # Test for allocations - @test (testallocs_iterate(ii_sdh)) === nothing - @test (testallocs_iterate(ii_same_sdh)) === nothing + # TODO: find a way to avoid resize for single sdh + # @test (testallocs_iterate(ii_dh_top)) === nothing + # @test (testallocs_iterate(ii_dh)) === nothing + # Iterators over grid can allocate due to potential resize! + @test_throws AllocCheckFailure testallocs_iterate(ii_grid) + @test_throws AllocCheckFailure testallocs_iterate(ii_grid_top) end - @testset "mixed cell types" begin + @testset "subdomains" begin + @testset "same cell type" begin + grid = generate_grid(Hexahedron, (1000, 2, 1)) + ip = Lagrange{RefHexahedron, 1}() + dh = DofHandler(grid) + sdh1 = SubDofHandler(dh, Set([collect(1:500)..., collect(1501:2000)...])) + sdh2 = SubDofHandler(dh, Set([collect(501:1500)...])) + add!(sdh1, :u, ip) + add!(sdh2, :u, ip) + close!(dh) + topology = ExclusiveTopology(grid) + + ii_sdh_top = InterfaceIterator(sdh1, sdh2, topology) + ii_sdh = InterfaceIterator(sdh1, sdh2) + ii_same_sdh = InterfaceIterator(sdh1, sdh1) + @test ii_sdh.set == ii_sdh_top.set + @test length(ii_sdh.set) == 1002 + @test count(interface_index -> interface_index[1] == interface_index[3] - 1000, ii_sdh.set) == 500 + @test count(interface_index -> interface_index[1] == interface_index[3] + 1000, ii_sdh.set) == 500 + ii_sdh_flipped = InterfaceIterator(sdh2, sdh1) + @test all(interface -> InterfaceIndex(interface[3], interface[4], interface[1], interface[2]) ∈ ii_sdh_flipped.set, ii_sdh.set) + # Test that topology has no effect on iterator type + @test typeof(ii_sdh) == typeof(ii_sdh_top) + # Test for allocations + @test (testallocs_iterate(ii_sdh)) === nothing + @test (testallocs_iterate(ii_same_sdh)) === nothing + end + + @testset "mixed cell types" begin + nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] + cells = [ + Quadrilateral((1, 2, 5, 4)), + Triangle((3, 5, 2)), + ] + grid = Grid(cells, nodes) + topology = ExclusiveTopology(grid) + ip_quad = Lagrange{RefQuadrilateral, 1}() + ip_tri = Lagrange{RefTriangle, 1}() + dh = DofHandler(grid) + sdh_quad = SubDofHandler(dh, Set([1])) + sdh_tri = SubDofHandler(dh, Set([2])) + add!(sdh_quad, :u, ip_quad) + add!(sdh_tri, :u, ip_tri) + close!(dh) + + ii_sdh_top = InterfaceIterator(sdh_quad, sdh_tri, topology) + ii_sdh = InterfaceIterator(sdh_quad, sdh_tri) + @test ii_sdh.set == ii_sdh_top.set + @test ii_sdh.set == Set([InterfaceIndex(1, 2, 2, 2)]) + ii_sdh_flipped = InterfaceIterator(sdh_tri, sdh_quad) + @test ii_sdh_flipped.set == Set([InterfaceIndex(2, 2, 1, 2)]) + # Test that topology has no effect on iterator type + @test typeof(ii_sdh) == typeof(ii_sdh_top) + # Test for allocations + @test (testallocs_iterate(ii_sdh_top)) === nothing + @test (testallocs_iterate(ii_sdh)) === nothing + end + end + + @testset "InterfaceIndex" begin + # checkmate, codecov bot! + @test isequal(InterfaceIndex(1, 2, 3, 4), (1, 2, 3, 4)) + @test isequal((1, 2, 3, 4), InterfaceIndex(1, 2, 3, 4)) + end + + @testset "error paths" begin nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] cells = [ Quadrilateral((1, 2, 5, 4)), @@ -251,43 +296,7 @@ end add!(sdh_quad, :u, ip_quad) add!(sdh_tri, :u, ip_tri) close!(dh) - - ii_sdh_top = InterfaceIterator(sdh_quad, sdh_tri, topology) - ii_sdh = InterfaceIterator(sdh_quad, sdh_tri) - @test ii_sdh.set == ii_sdh_top.set - @test ii_sdh.set == Set([InterfaceIndex(1, 2, 2, 2)]) - ii_sdh_flipped = InterfaceIterator(sdh_tri, sdh_quad) - @test ii_sdh_flipped.set == Set([InterfaceIndex(2, 2, 1, 2)]) - # Test that topology has no effect on iterator type - @test typeof(ii_sdh) == typeof(ii_sdh_top) - # Test for allocations - @test (testallocs_iterate(ii_sdh_top)) === nothing - @test (testallocs_iterate(ii_sdh)) === nothing + @test_throws ErrorException("The cells in the set (set of InterfaceIndex) are not all of the same celltype on each side.") InterfaceIterator(sdh_quad, sdh_tri, Set([InterfaceIndex(2, 2, 1, 2)])) end end - - @testset "InterfaceIndex" begin - # checkmate, codecov bot! - @test isequal(InterfaceIndex(1, 2, 3, 4), (1, 2, 3, 4)) - @test isequal((1, 2, 3, 4), InterfaceIndex(1, 2, 3, 4)) - end - - @testset "error paths" begin - nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))] - cells = [ - Quadrilateral((1, 2, 5, 4)), - Triangle((3, 5, 2)), - ] - grid = Grid(cells, nodes) - topology = ExclusiveTopology(grid) - ip_quad = Lagrange{RefQuadrilateral, 1}() - ip_tri = Lagrange{RefTriangle, 1}() - dh = DofHandler(grid) - sdh_quad = SubDofHandler(dh, Set([1])) - sdh_tri = SubDofHandler(dh, Set([2])) - add!(sdh_quad, :u, ip_quad) - add!(sdh_tri, :u, ip_tri) - close!(dh) - @test_throws ErrorException("The cells in the set (set of InterfaceIndex) are not all of the same celltype on each side.") InterfaceIterator(sdh_quad, sdh_tri, Set([InterfaceIndex(2, 2, 1, 2)])) - end end From 720d373905a30f8db364d6485ae5337b0d934835 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 13 Nov 2024 20:50:48 +0100 Subject: [PATCH 28/28] avoid looping over interfaces twice for subdomains in sparsitypatterns --- src/Dofs/sparsity_pattern.jl | 2 +- src/iterators.jl | 8 +++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/Dofs/sparsity_pattern.jl b/src/Dofs/sparsity_pattern.jl index d9db08f780..bc1d14aef1 100644 --- a/src/Dofs/sparsity_pattern.jl +++ b/src/Dofs/sparsity_pattern.jl @@ -610,7 +610,7 @@ function _add_interface_entries!( couplings = _coupling_to_local_dof_coupling(dh, interface_coupling) for (_i, sdh1) in enumerate(dh.subdofhandlers) for (_j, sdh2) in enumerate(dh.subdofhandlers) - for ic in InterfaceIterator(sdh1, sdh2, topology) + for ic in InterfaceIterator(sdh1, sdh2, topology, false) # TODO: This looks like it can be optimized for the common case where # the cells are in the same subdofhandler to_check = Dict{Int, Vector{Int}}() diff --git a/src/iterators.jl b/src/iterators.jl index c966a158af..169ccd48d8 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -443,7 +443,9 @@ end function InterfaceIterator( sdh_here::SubDofHandler, sdh_there::SubDofHandler, - topology::ExclusiveTopology = ExclusiveTopology(get_grid(sdh_here.dh)) + topology::ExclusiveTopology = ExclusiveTopology(get_grid(sdh_here.dh)), + # TODO: better name? + allow_all::Bool = true # allow interfaces that belong to the other cell according to the cell numbering to be iterated over as if they were "here" not "there" ) grid = get_grid(sdh_here.dh) neighborhood = get_facet_facet_neighborhood(topology, grid) @@ -454,7 +456,7 @@ function InterfaceIterator( neighbors = neighborhood[facet[1], facet[2]] isempty(neighbors) && continue neighbors[][1] ∈ sdh_there.cellset || continue - elseif facet[1] ∈ sdh_there.cellset + elseif allow_all && facet[1] ∈ sdh_there.cellset neighbors = neighborhood[facet[1], facet[2]] isempty(neighbors) && continue neighbors[][1] ∈ sdh_here.cellset || continue @@ -471,7 +473,7 @@ function InterfaceIterator( isempty(neighbors) && continue neighbors[][1] ∈ sdh_there.cellset || continue push!(set, InterfaceIndex(facet[1], facet[2], neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2])) - elseif facet[1] ∈ sdh_there.cellset + elseif allow_all && facet[1] ∈ sdh_there.cellset neighbors = neighborhood[facet[1], facet[2]] isempty(neighbors) && continue neighbors[][1] ∈ sdh_here.cellset || continue