Merge pull request #86 from WIAS-PDELib/feature/multi-coupling

`CoupledDofsRestriction`: allow multiple matrices, reduce cols if redundant

Author: pjaap

CHANGELOG.md

@@ -1,5 +1,15 @@
 # CHANGES
 
+## v1.6.0
+
+### Added
+
+  - `CoupledDofsRestriction` accepts now an array of multiple coupling matrices.
+
+### Fixed
+
+  - Periodic coupling with multiple redundant coupling matrices are now reduced to full rank to make the system uniquely solvable.
+
 ## v1.5.0
 
 ### Added

Project.toml

@@ -1,6 +1,6 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
-version = "1.5.0"
+version = "1.6.0"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
 
 [deps]

src/ExtendableFEM.jl

@@ -70,7 +70,7 @@ using LinearSolve: LinearSolve, LinearProblem, UMFPACKFactorization, deleteat!,
     init, solve
 using Printf: Printf, @printf, @sprintf
 using SparseArrays: SparseArrays, AbstractSparseArray, SparseMatrixCSC, findnz, nnz,
-    nzrange, rowvals, sparse, SparseVector, spzeros
+    nzrange, rowvals, sparse, SparseVector, spzeros, qr, rank
 using StaticArrays: @MArray
 using SparseDiffTools: SparseDiffTools, ForwardColorJacCache,
     forwarddiff_color_jacobian!, matrix_colors

src/common_restrictions/coupled_dofs_restriction.jl

@@ -1,13 +1,13 @@
 struct CoupledDofsRestriction{TM} <: AbstractRestriction
-    coupling_matrix::TM
+    coupling_matrices::Vector{TM}
     parameters::Dict{Symbol, Any}
 end
 
 
 """
     CoupledDofsRestriction(matrix::AbstractMatrix)
 
-    Creates an restriction that couples dofs together.
+    Creates a restriction that couples dofs together.
 
     The coupling is stored in a CSC Matrix `matrix`, s.t.,
 
@@ -16,27 +16,58 @@ end
     The matrix can be obtained from, e.g., `get_periodic_coupling_matrix`.
 """
 function CoupledDofsRestriction(matrix::AbstractMatrix)
-    return CoupledDofsRestriction(matrix, Dict{Symbol, Any}(:name => "CoupledDofsRestriction"))
+    return CoupledDofsRestriction(
+        [matrix],
+        Dict{Symbol, Any}(
+            :name => "CoupledDofsRestriction",
+            :reduce_col_space => false
+        )
+    )
 end
 
 
-function assemble!(R::CoupledDofsRestriction, sol, SC; kwargs...)
+"""
+    CoupledDofsRestriction(matrices::Vector{AM}) where {AM <: AbstractMatrix}
+
+    Creates a `CoupledDofsRestriction` from multiple given coupling matrices.
+
+    By default, the column space of the matrices is reduced to be of full rank.
+    This can toggled by the `:reduce_col_space` parameter.
+"""
+function CoupledDofsRestriction(matrices::Vector{AM}) where {AM <: AbstractMatrix}
+    return CoupledDofsRestriction(
+        matrices,
+        Dict{Symbol, Any}(
+            :name => "CoupledDofsRestriction",
+            :reduce_col_space => true
+        )
+    )
+end
 
-    # extract all col indices
-    _, J, _ = findnz(R.coupling_matrix)
 
-    # remove duplicates
-    unique_cols = unique(J)
+function assemble!(R::CoupledDofsRestriction, sol, SC; kwargs...)
 
+    # extract all col indices and remove duplicates
     # subtract diagonal and shrink matrix to non-empty cols
-    B = (R.coupling_matrix - LinearAlgebra.I)[:, unique_cols]
+    Bs = [ (matrix - LinearAlgebra.I)[:, unique(findnz(matrix)[2])] for matrix in R.coupling_matrices ]
+
+    # combine all into one matrix
+    B = hcat(Bs...)
+
+    if R.parameters[:reduce_col_space]
+        # eliminate redundant cols by QR:
+        qr_result = qr(B)
+
+        # pick minimal number of cols that are rank preserving
+        cols_of_interest = qr_result.pcol[1:rank(qr_result)]
+        B = B[:, cols_of_interest]
+    end
 
     R.parameters[:matrix] = B
-    R.parameters[:rhs] = Zeros(length(unique_cols))
+    R.parameters[:rhs] = Zeros(size(B, 2))
 
     # fixed dofs are all active rows of B
-    I, _, _ = findnz(B)
-    R.parameters[:fixed_dofs] = unique(I)
+    R.parameters[:fixed_dofs] = unique(findnz(B)[1])
 
     return nothing
 end

src/helper_functions.jl

@@ -371,7 +371,7 @@ function _get_periodic_coupling_matrix(
     # throw error if no search area had been found for a bface
     for source in 1:num_sources(searchareas)
         if num_targets(searchareas, source) == 0
-            throw("bface $source has no valid search area on the opposite side of the grid. Are from/to regions and give_opposite! function correct?")
+            throw("bface $source has no valid search area on the opposite side of the grid. Double check the provided from/to regions and your give_opposite! function")
         end
     end
 
@@ -499,7 +499,7 @@ For each x in the grid, the resulting y has to be in the grid, too: incorporate
 Example: If b_from is at x[1] = 0 and the opposite boundary is at y[1] = 1, then give_opposite!(y,x) = y .= [ 1-x[1], x[2] ]
 
 The return value is a (𝑛 × 𝑛) sparse matrix 𝐴 (𝑛 is the total number of dofs) containing the periodic coupling information.
-The relation ship between the degrees of freedome is  dofᵢ = ∑ⱼ Aⱼᵢ ⋅ dofⱼ.
+The relation ship between the degrees of freedom is  dofᵢ = ∑ⱼ Aⱼᵢ ⋅ dofⱼ.
 It is guaranteed that
     i)  Aᵢⱼ=0 if dofᵢ is 𝑛𝑜𝑡 on the boundary b_from.
     ii) Aᵢⱼ=0 if the opposite of dofᵢ is not in the same grid cell as dofⱼ.