Cross-elements coupling for DiscontinuousLagrange sparsity patterns.

docs/src/devdocs/interpolations.md

@@ -31,6 +31,8 @@ Ferrite.edgedof_interior_indices(::Interpolation)
 Ferrite.celldof_interior_indices(::Interpolation)
 Ferrite.getnbasefunctions(::Interpolation)
 Ferrite.reference_coordinates(::Interpolation)
+Ferrite.IsDiscontinuous(::Interpolation)
+Ferrite.get_continuous_interpolation(::Interpolation)
 ```
 
 for all entities which exist on that reference element. The dof functions default to having no

src/Dofs/sparsity_pattern.jl

@@ -1,5 +1,5 @@
 """
-    create_sparsity_pattern(dh::DofHandler; coupling)
+    create_sparsity_pattern(dh::DofHandler; coupling, topology::Union{Nothing, AbstractTopology} = nothing)
 
 Create the sparsity pattern corresponding to the degree of freedom
 numbering in the [`DofHandler`](@ref). Return a `SparseMatrixCSC`
@@ -12,44 +12,82 @@ in the DofHandler with `true` if fields are coupled and `false` if
 not. By default full coupling is assumed.
 
 See the [Sparsity Pattern](@ref) section of the manual.
+
+```julia-repl
+julia> using Ferrite
+
+julia> grid = generate_grid(Line, (3,));
+
+julia> topology = ExclusiveTopology(grid);
+
+julia> ip = DiscontinuousLagrange{1, RefCube, 2}();
+
+julia> ipc = Lagrange{1, RefCube, 1}();
+
+julia> dh = DofHandler(grid);
+
+julia> add!(dh, :u, 1,ip);
+
+julia> add!(dh, :v, 1,ipc);
+
+julia> close!(dh);;
+
+julia> K = create_sparsity_pattern(dh, topology = topology)
+13×13 SparseMatrixCSC{Float64, Int64} with 109 stored entries:
+ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   ⋅    ⋅    ⋅    ⋅    ⋅
+ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   ⋅    ⋅    ⋅    ⋅    ⋅
+ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   ⋅    ⋅    ⋅    ⋅    ⋅
+ 0.0  0.0  0.0  0.0  0.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
+ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   ⋅    ⋅    ⋅    ⋅
+ 0.0  0.0  0.0   ⋅   0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   ⋅
+ 0.0  0.0  0.0   ⋅   0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   ⋅
+ 0.0  0.0  0.0   ⋅   0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   ⋅
+  ⋅    ⋅    ⋅    ⋅   0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
+  ⋅    ⋅    ⋅    ⋅    ⋅   0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
+  ⋅    ⋅    ⋅    ⋅    ⋅   0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
+  ⋅    ⋅    ⋅    ⋅    ⋅   0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
+  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   0.0  0.0  0.0  0.0  0.0
+```
 """
-function create_sparsity_pattern(dh::AbstractDofHandler; coupling=nothing)
-    return _create_sparsity_pattern(dh, nothing, false, true, coupling)
+function create_sparsity_pattern(dh::AbstractDofHandler; coupling=nothing,
+    topology::Union{Nothing, AbstractTopology} = nothing)
+    return _create_sparsity_pattern(dh, nothing, false, true, coupling; topology)
 end
 
 """
-    create_symmetric_sparsity_pattern(dh::DofHandler; coupling)
+    create_symmetric_sparsity_pattern(dh::DofHandler; coupling, topology::Union{Nothing, AbstractTopology} = nothing)
 
 Create the symmetric sparsity pattern corresponding to the degree of freedom
 numbering in the [`DofHandler`](@ref) by only considering the upper
 triangle of the matrix. Return a `Symmetric{SparseMatrixCSC}`.
 
 See the [Sparsity Pattern](@ref) section of the manual.
 """
-function create_symmetric_sparsity_pattern(dh::AbstractDofHandler; coupling=nothing)
-    return Symmetric(_create_sparsity_pattern(dh, nothing, true, true, coupling), :U)
+function create_symmetric_sparsity_pattern(dh::AbstractDofHandler; coupling=nothing,
+    topology::Union{Nothing, AbstractTopology} = nothing)
+    return Symmetric(_create_sparsity_pattern(dh, nothing, true, true, coupling; topology), :U)
 end
 
 """
-    create_symmetric_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler, coupling)
+    create_symmetric_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler, coupling, topology::Union{Nothing, AbstractTopology} = nothing)
 
 Create a symmetric sparsity pattern accounting for affine constraints in `ch`. See
 the Affine Constraints section of the manual for further details.
 """
 function create_symmetric_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler;
-        keep_constrained::Bool=true, coupling=nothing)
-    return Symmetric(_create_sparsity_pattern(dh, ch, true, keep_constrained, coupling), :U)
+        keep_constrained::Bool=true, coupling=nothing, topology::Union{Nothing, AbstractTopology} = nothing)
+    return Symmetric(_create_sparsity_pattern(dh, ch, true, keep_constrained, coupling; topology), :U)
 end
 
 """
-    create_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler; coupling)
+    create_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler; coupling, topology::Union{Nothing, AbstractTopology} = nothing)
 
 Create a sparsity pattern accounting for affine constraints in `ch`. See
 the Affine Constraints section of the manual for further details.
 """
 function create_sparsity_pattern(dh::AbstractDofHandler, ch::ConstraintHandler;
-        keep_constrained::Bool=true, coupling=nothing)
-    return _create_sparsity_pattern(dh, ch, false, keep_constrained, coupling)
+        keep_constrained::Bool=true, coupling=nothing, topology::Union{Nothing, AbstractTopology} = nothing)
+    return _create_sparsity_pattern(dh, ch, false, keep_constrained, coupling; topology)
 end
 
 # Compute a coupling matrix of size (ndofs_per_cell × ndofs_per_cell) based on the input
@@ -96,21 +134,78 @@ function _coupling_to_local_dof_coupling(dh::DofHandler, coupling::AbstractMatri
     return outs
 end
 
-function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{ConstraintHandler, Nothing}=#, sym::Bool, keep_constrained::Bool, coupling::Union{AbstractMatrix{Bool},Nothing})
+"""
+    cross_element_coupling!(dh::AbstractDofHandler, topology::ExclusiveTopology, sym::Bool, keep_constrained::Bool, couplings::Union{AbstractVector{<:AbstractMatrix{Bool}},Nothing}, cnt::Int, I::Vector{Int}, J::Vector{Int})
+Calculates `I, J` for cross-element coupling
+Returns the updated value of `cnt`
+Used internally for sparsity patterns with discontinuous interpolations.
+"""
+function cross_element_coupling!(dh::AbstractDofHandler, topology::ExclusiveTopology, sym::Bool, keep_constrained::Bool, couplings::Union{AbstractVector{<:AbstractMatrix{Bool}},Nothing}, cnt::Int, I::Vector{Int}, J::Vector{Int})
+    for (fhi, fh) in pairs(dh.fieldhandlers)
+        ip_infos = InterpolationInfo[]
+        nbasefunctions = Int[]
+        for cell_field in fh.fields
+            ip_info = InterpolationInfo(cell_field.interpolation)
+            push!(ip_infos, ip_info)
+            push!(nbasefunctions, getnbasefunctions(cell_field.interpolation))
+        end
+        element_dof_start = 0
+        isnothing(couplings) || (coupling_fh = couplings[fhi])
+        for (cell_field_i, cell_field) in enumerate(fh.fields)
+            fii = ip_infos[cell_field_i]
+            if(!fii.is_discontinuous)
+                element_dof_start += nbasefunctions[cell_field_i]
+                continue
+            end
+            for cell_idx in eachindex(getcells(dh.grid))
+                cell_field ∈ dh.fieldhandlers[dh.cell_to_fieldhandler[cell_idx]].fields || continue
+                cell_dofs = @view dh.cell_dofs[dh.cell_dofs_offset[cell_idx] : dh.cell_dofs_offset[cell_idx] + ndofs_per_cell(dh, cell_idx) - 1]
+                cell_field_dofs = @view cell_dofs[element_dof_start + 1 : element_dof_start + nbasefunctions[cell_field_i]]
+                for neighbor_cell in (getdim(dh.grid.cells[cell_idx]) >1 ? topology.cell_face_neighbor[cell_idx] : topology.cell_neighbor[cell_idx])
+                    neighbour_dof_start = 0
+                    for (neighbor_field_i, neighbor_field) in enumerate(fh.fields)
+                        fii2 = ip_infos[neighbor_field_i]
+                        neighbor_field ∈ dh.fieldhandlers[dh.cell_to_fieldhandler[neighbor_cell.idx]].fields || continue
+                        if(!fii2.is_discontinuous)
+                            neighbour_dof_start += nbasefunctions[neighbor_field_i]
+                            continue
+                        end
+                        neighbor_dofs = @view dh.cell_dofs[dh.cell_dofs_offset[neighbor_cell.idx] : dh.cell_dofs_offset[neighbor_cell.idx] + ndofs_per_cell(dh, neighbor_cell.idx) - 1]
+                        neighbor_field_dofs = @view neighbor_dofs[neighbour_dof_start + 1 : neighbour_dof_start + nbasefunctions[neighbor_field_i]]
+                        for j in eachindex(neighbor_field_dofs), i in eachindex(cell_field_dofs)
+                            isnothing(couplings) || coupling_fh[i+element_dof_start,j+neighbour_dof_start] || continue
+                            dofi = cell_field_dofs[i]
+                            dofj = neighbor_field_dofs[j]
+                            sym && (dofi > dofj && continue)
+                            !keep_constrained && (haskey(ch.dofmapping, dofi) || haskey(ch.dofmapping, dofj)) && continue
+                            cnt += 1
+                            _add_or_grow(cnt, I, J, dofi, dofj)
+                        end
+                        neighbour_dof_start += nbasefunctions[neighbor_field_i]
+                    end
+                end
+            end
+            element_dof_start += nbasefunctions[cell_field_i]
+        end
+    end
+    return cnt
+end
+
+function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{ConstraintHandler, Nothing}=#, sym::Bool, keep_constrained::Bool, coupling::Union{AbstractMatrix{Bool},Nothing};
+    topology::Union{Nothing, AbstractTopology} = nothing)
     @assert isclosed(dh)
     if !keep_constrained
         @assert ch !== nothing && isclosed(ch)
     end
-    if coupling !== nothing
-        # Extend coupling to be of size (ndofs_per_cell × ndofs_per_cell)
-        couplings = _coupling_to_local_dof_coupling(dh, coupling, sym)
-    end
+
+    couplings = isnothing(coupling) ? nothing : _coupling_to_local_dof_coupling(dh, coupling, sym)
 
     # Allocate buffers. Compute an upper bound for the buffer length and allocate it all up
     # front since they will become large and expensive to re-allocate. The bound is exact
     # when keeping constrained dofs (default) and if not it only over-estimates with number
     # of entries eliminated by constraints.
     max_buffer_length = ndofs(dh) # diagonal elements
+    has_discontinuous_ip = false
     for (fhi, fh) in pairs(dh.fieldhandlers)
         set = fh.cellset
         n = ndofs_per_cell(fh)
@@ -120,6 +215,7 @@ function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{Constraint
             coupling_fh = couplings[fhi]
             count(coupling_fh[i, j] for i in 1:n for j in (sym ? i : 1):n)
         end
+        has_discontinuous_ip = any(ip -> IsDiscontinuous(ip),fh.field_interpolations)
         max_buffer_length += entries_per_cell * length(set)
     end
     I = Vector{Int}(undef, max_buffer_length)
@@ -147,7 +243,7 @@ function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{Constraint
             end
         end
     end
-
+    has_discontinuous_ip && (cnt = cross_element_coupling!(dh,topology,sym, keep_constrained, couplings, cnt, I, J))
     # Always add diagonal entries
     resize!(I, cnt + ndofs(dh))
     resize!(J, cnt + ndofs(dh))