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

Project.toml

@@ -54,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"]

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

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]

src/Dofs/sparsity_pattern.jl

@@ -446,51 +446,53 @@ 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
 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

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/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)

src/exports.jl

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