Merge pull request #114 from WIAS-PDELib/chmerdon/faster_combinedofs_operator

faster combinedofs operator assembly

Author: pjaap

CHANGELOG.md

@@ -1,5 +1,12 @@
 # CHANGES
 
+
+## v.1.10.4
+
+### Fixed
+  - assembly of CombineDofs operator now much faster
+
+
 ## v.1.10.3
 
 ### Fixed

Project.toml

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

examples/Example312_PeriodicElasticity3D.jl

@@ -129,7 +129,7 @@ end
 
 function main(;
         order = 1,
-        periodic_coupling = :none, # :restriction, :operator, :high_level_restriction
+        periodic_coupling = :high_level_restriction, # :restriction, :operator, :high_level_restriction
         Plotter = nothing,
         force = 1.0,
         h = 1.0e-4,

src/ExtendableFEM.jl

@@ -61,7 +61,7 @@ using ExtendableGrids: ExtendableGrids, AT_NODES, AbstractElementGeometry,
     unique, update_trafo!, xrefFACE2xrefCELL,
     xrefFACE2xrefOFACE
 using ExtendableSparse: ExtendableSparse, ExtendableSparseMatrix, flush!,
-    MTExtendableSparseMatrixCSC, findindex,
+    MTExtendableSparseMatrixCSC,
     rawupdateindex!
 using ForwardDiff: ForwardDiff
 using GridVisualize: GridVisualize, GridVisualizer, gridplot!, reveal, save,

src/common_operators/combinedofs.jl

@@ -107,29 +107,27 @@ 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
+
+                @views AT = sparse(A.cscmatrix[fixed_dofs .+ offsetX, :]')
+
                 # go through each constrained dof and update the FE adjacency info
                 # of the coupled dofs
-                for dof_i in fixed_dofs
+                for (i, dof_i) in enumerate(fixed_dofs)
                     # this col-view is efficient
                     coupling_i = @views coupling_matrix[:, dof_i]
 
-                    # write the FE adjacency of the constrained dofs into this row
-                    sourcerow = dof_i + offsetX
-
                     # extract the constrained dofs and the weights
                     coupled_dofs_i, weights_i = findnz(coupling_i)
 
                     # parse through sourcerow and add the contents to the coupled dofs
-                    for col in 1:size(A, 2)
-                        r = findindex(A.cscmatrix, sourcerow, col)
-                        if r > 0
-                            val = A.cscmatrix.nzval[r]
-                            if abs(val) > 1.0e-15
-                                for (dof_k, weight_ik) in zip(coupled_dofs_i, weights_i)
-                                    targetrow = dof_k + offsetX
-                                    _addnz(A, targetrow, col, val, weight_ik)
-                                end
-                            end
+                    A_row = @views AT[:, i]
+                    cols_i, vals_i = findnz(A_row)
+
+                    # copy values of sourcerow to the coupled targetrow
+                    for (col, val) in zip(cols_i, vals_i)
+                        for (dof_k, weight_ik) in zip(coupled_dofs_i, weights_i)
+                            targetrow = dof_k + offsetX
+                            _addnz(A, targetrow, col, val, weight_ik)
                         end
                     end
                 end