Subdomain support for Interfaces + AllocCheck tests for iterators + Sparsitypattern bugfix

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"

src/Ferrite.jl

@@ -107,6 +107,13 @@ 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} # cell - side - cell - side
+end
+
 const AbstractVecOrSet{T} = Union{AbstractSet{T}, AbstractVector{T}}
 const IntegerCollection = AbstractVecOrSet{<:Integer}
 

src/Grid/grid.jl

@@ -641,20 +641,20 @@ 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]
 
         #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
-        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

src/exports.jl

@@ -77,6 +77,7 @@ export
     EdgeIndex,
     VertexIndex,
     FacetIndex,
+    InterfaceIndex,
     geometric_interpolation,
     ExclusiveTopology,
     getneighborhood,

src/iterators.jl

@@ -61,12 +61,17 @@ 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)
+# 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))
@@ -83,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)
@@ -164,22 +184,27 @@ 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
 
 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 reinit!(cache::InterfaceCache, facet_a::BoundaryIndex, facet_b::BoundaryIndex)
-    reinit!(cache.a, facet_a)
-    reinit!(cache.b, facet_b)
-    resize!(cache.dofs, length(celldofs(cache.a)) + length(celldofs(cache.b)))
+function InterfaceCache(sdh_here::SubDofHandler, sdh_there::SubDofHandler)
+    fc_a = FacetCache(sdh_here)
+    fc_b = FacetCache(sdh_there)
+    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]))
     for (i, d) in pairs(cache.a.dofs)
         cache.dofs[i] = d
     end
@@ -341,40 +366,98 @@ 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, SetType <: AbstractSet{InterfaceIndex}}
     cache::IC
-    grid::G
-    topology::ExclusiveTopology
+    set::SetType
 end
 
+@inline _getcache(ii::InterfaceIterator) = ii.cache
+@inline _getset(ii::InterfaceIterator) = ii.set
+
 function InterfaceIterator(
-        gridordh::Union{Grid, AbstractDofHandler},
+        gridordh::Union{Grid, DofHandler},
+        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)
+    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, getcelltype(grid, first(sdh_here.cellset)), getcelltype(grid, first(sdh_there.cellset)))
+    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))
     )
     grid = gridordh isa Grid ? gridordh : get_grid(gridordh)
-    return InterfaceIterator(InterfaceCache(gridordh), grid, topology)
+    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 = Set{InterfaceIndex}()
+    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]))
+    end
+    return InterfaceIterator(gridordh, set)
 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]])
+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
+        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
-        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)
+        ninterfaces += 1
     end
-    return
+    set = Set{InterfaceIndex}()
+    sizehint!(set, ninterfaces)
+    for facet in fs
+        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(sdh_here, sdh_there, set)
 end
 
-
 # Iterator interface for CellIterator/FacetIterator
 const GridIterators{C} = Union{CellIterator{C}, FacetIterator{C}, InterfaceIterator{C}}
 
@@ -408,3 +491,15 @@ function _check_same_celltype(grid::AbstractGrid, facetset::AbstractVecOrSet{<:B
     end
     return
 end
+
+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
+    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

test/runtests.jl

@@ -11,6 +11,7 @@ using StaticArrays
 using OrderedCollections
 using WriteVTK
 import Metis
+using AllocCheck
 
 include("test_utils.jl")
 
@@ -20,6 +21,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")

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(