Implement Christian's brilliant finding

Author: pjaap

src/common_operators/combinedofs.jl

@@ -87,71 +87,59 @@ function build_assembler!(CD::CombineDofs{UT, CT}, FE::Array{<:FEVectorBlock, 1}
         function assemble!(A::AbstractSparseArray{T}, b::AbstractVector{T}, assemble_matrix::Bool, assemble_rhs::Bool, kwargs...) where {T}
             if assemble_matrix
 
-                # put FE adjacency information into this matrix row
-                for dof_i in 1:size(coupling_matrix, 2)
+                # transpose the matrix once for efficient row access
+                transposed_coupling_matrix = sparse(transpose(coupling_matrix))
+
+                # go through each coupled dof and update the FE adjacency info
+                # from the constrained dofs here
+
+                for dof_i in 1:size(transposed_coupling_matrix, 2)
                     # this col-view is efficient
-                    coupling_i = @views coupling_matrix[:, dof_i]
-                    # do nothing if no coupling for dof_i
+                    coupling_i = @views transposed_coupling_matrix[:, dof_i]
+                    # do nothing if dof_k is not coupled to any constrained dof
                     if nnz(coupling_i) == 0
                         continue
                     end
 
+                    # write the FE adjacency of the constrained dofs into this row
                     targetrow = dof_i + offsetX
 
+                    # extract the constrained dofs and the weights
                     coupled_dofs_i, weights_i = findnz(coupling_i)
 
-                    for dof_j in 1:size(coupling_matrix, 2)
+                    # parse through all cols and update the entries
+                    for dof_j in 1:size(transposed_coupling_matrix, 2)
                         # this col-view is efficient
-                        coupling_j = @views coupling_matrix[:, dof_j]
+                        coupling_j = @views transposed_coupling_matrix[:, dof_j]
 
+                        # if both dof_i and dof_j are coupled to a constrained dof, then
+                        # the FE adjacency A_ij is not updated: this is covered by the linear combinations
+                        # expressed in the rows of the constrained dofs_on_boundary
+                        # Hence, check that dof_j is not coupled to anything
                         if nnz(coupling_j) == 0
-                            # not both i and j on the boundary
-
                             targetcol = dof_j + offsetY
-                            for (dof_k, weight_k) in zip(coupled_dofs_i, weights_i)
+                            for (dof_k, weight_ik) in zip(coupled_dofs_i, weights_i)
                                 sourcerow = dof_k + offsetX
                                 sourcecol = targetcol
                                 val = A[sourcerow, sourcecol]
-                                _addnz(A, targetrow, targetcol, val, weight_k)
-                            end
-
-
-                        else
-                            # i and j on the same boundary
-
-                            coupled_dofs_j, weights_j = findnz(coupling_j)
-                            targetcol = dof_j + offsetY
-
-                            for (dof_k, weight_k) in zip(coupled_dofs_i, weights_i)
-                                for (dof_l, weight_l) in zip(coupled_dofs_j, weights_j)
-
-                                    sourcerow = dof_k + offsetX
-                                    sourcecol = dof_l + offsetY
-
-                                    val = A[sourcerow, sourcecol]
-                                    # @show sourcerow, sourcecol, targetrow, targetcol, val, weight_k, weight_l
-                                    _addnz(A, targetrow, targetcol, val, weight_k * weight_l)
-
-                                end
+                                _addnz(A, targetrow, targetcol, val, weight_ik)
                             end
                         end
                     end
                 end
 
-                # replace the geometric coupling rows
+                # replace the geometric coupling rows based
+                # on the original coupling matrix
                 for dof_i in 1:size(coupling_matrix, 2)
-                    # go through the rows of the matrix; use transposed matrix implicitly
 
-                    # this row-view is not efficient :(
-                    coupling_i = coupling_matrix[dof_i, :]
+                    coupling_i = coupling_matrix[:, dof_i]
                     # do nothing if no coupling for dof_i
                     if nnz(coupling_i) == 0
                         continue
                     end
 
-
                     # get the coupled dofs of dof_i and the corresponding weights
-                    coupled_dofs, weights = findnz(coupling_i)
+                    coupled_dofs_i, weights_i = findnz(coupling_i)
 
                     sourcerow = dof_i + offsetX
 
@@ -162,9 +150,9 @@ function build_assembler!(CD::CombineDofs{UT, CT}, FE::Array{<:FEVectorBlock, 1}
 
                     # replace sourcerow with coupling linear combination
                     _addnz(A, sourcerow, sourcerow, -1.0, 1)
-                    for (dof_j, weight) in zip(coupled_dofs, weights)
+                    for (dof_j, weight_ij) in zip(coupled_dofs_i, weights_i)
                         # weights for ∑ⱼ wⱼdofⱼ - dofᵢ = 0
-                        _addnz(A, sourcerow, dof_j + offsetY, weight, 1)
+                        _addnz(A, sourcerow, dof_j + offsetY, weight_ij, 1)
                     end
 
                 end
@@ -173,9 +161,9 @@ function build_assembler!(CD::CombineDofs{UT, CT}, FE::Array{<:FEVectorBlock, 1}
 
             if assemble_rhs
 
-                for dof_i in 1:size(coupling_matrix, 2)
+                for dof_i in 1:size(transposed_coupling_matrix, 2)
                     # this col-view is efficient
-                    coupling_i = @views coupling_matrix[:, dof_i]
+                    coupling_i = @views transposed_coupling_matrix[:, dof_i]
                     # do nothing if no coupling for dof_i
                     if nnz(coupling_i) == 0
                         continue
@@ -193,16 +181,15 @@ function build_assembler!(CD::CombineDofs{UT, CT}, FE::Array{<:FEVectorBlock, 1}
                 end
 
 
+                # now set the rows of the constrained dofs to zero to enforce the linear combination
                 for dof_i in 1:size(coupling_matrix, 2)
                     coupling_i = coupling_matrix[dof_i, :]
                     # do nothing if no coupling for dof_i
                     if nnz(coupling_i) == 0
                         continue
                     end
 
-                    sourcerow = dof_i + offsetX
-
-                    b[sourcerow] = 0.0
+                    b[dof_i + offsetX] = 0.0
 
                 end
             end