Cleanup: imports, variable names, timers

Author: pjaap

src/ExtendableFEM.jl

@@ -5,7 +5,7 @@ $(read(joinpath(@__DIR__, "..", "README.md"), String))
 """
 module ExtendableFEM
 
-using BlockArrays: BlockMatrix, BlockVector, Block, BlockedMatrix, blocks, axes
+using BlockArrays: BlockMatrix, BlockVector, Block, blocks, axes
 using CommonSolve: CommonSolve
 using DiffResults: DiffResults
 using DocStringExtensions: DocStringExtensions, TYPEDEF, TYPEDSIGNATURES

src/solvers.jl

@@ -223,7 +223,7 @@ Solves the linear system and updates the solution vector. This includes:
 """
 function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns, freedofs, damping, PD, SC, stats, is_linear, timer, kwargs...)
 
-    @timeit timer "Lagrange restrictions" begin
+    @timeit timer "assembly" begin
         ## assemble restrctions
         if !SC.parameters[:initialized]
             for restriction in PD.restrictions
@@ -255,67 +255,64 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
             linsolve.b = b_unrestricted
         else
 
-            @timeit timer "LM restrictions" begin
-                ## add possible Lagrange restrictions
-                @timeit timer "prepare" begin
-                    restriction_matrices = [length(freedofs) > 0 ? re.parameters[:matrix][freedofs, :] : re.parameters[:matrix] for re in PD.restrictions ]
-                    restriction_rhs = [length(freedofs) > 0 ? re.parameters[:rhs][freedofs] : re.parameters[:rhs] for re in PD.restrictions ]
+            # add possible Lagrange restrictions
+            restriction_matrices = [length(freedofs) > 0 ? re.parameters[:matrix][freedofs, :] : re.parameters[:matrix] for re in PD.restrictions ]
+            restriction_rhs = [length(freedofs) > 0 ? re.parameters[:rhs][freedofs] : re.parameters[:rhs] for re in PD.restrictions ]
 
-                    ## we need to add the (initial) solution to the rhs, since we work with the residual equation
-                    for (B, rhs) in zip(restriction_matrices, restriction_rhs)
-                        rhs .+= B'sol.entries
-                    end
-                end
+            # block sizes for the block matrix
+            block_sizes = [size(A_unrestricted, 2), [ size(B, 2) for B in restriction_matrices ]...]
 
-                @timeit timer "set blocks" begin
-                    # block sizes for the block matrix
-                    block_sizes = [size(A_unrestricted, 2), [ size(B, 2) for B in restriction_matrices ]...]
+            # total number of additional LM dofs
+            nLMs = @views sum(block_sizes[2:end])
 
-                    total_size = sum(block_sizes)
-                    Tv = eltype(A_unrestricted)
+            @timeit timer "LM restrictions (nLMs = $nLMs)" begin
 
-                    ## create block matrix
-                    A_block = BlockMatrix(spzeros(Tv, total_size, total_size), block_sizes, block_sizes)
-                    A_block[Block(1, 1)] = A_unrestricted
+                ## we need to add the (initial) solution to the rhs, since we work with the residual equation
+                for (B, rhs) in zip(restriction_matrices, restriction_rhs)
+                    rhs .+= B'sol.entries
+                end
 
-                    b_block = BlockVector(zeros(Tv, total_size), block_sizes)
-                    b_block[Block(1)] = b_unrestricted
+                total_size = sum(block_sizes)
+                Tv = eltype(A_unrestricted)
 
-                    u_unrestricted = linsolve.u
-                    u_block = BlockVector(zeros(Tv, total_size), block_sizes)
-                    u_block[Block(1)] = u_unrestricted
+                ## create block matrix
+                A_block = BlockMatrix(spzeros(Tv, total_size, total_size), block_sizes, block_sizes)
+                A_block[Block(1, 1)] = A_unrestricted
 
-                    for i in eachindex(PD.restrictions)
-                        A_block[Block(1, i + 1)] = restriction_matrices[i]
-                        A_block[Block(i + 1, 1)] = transpose(restriction_matrices[i])
-                        b_block[Block(i + 1)] = restriction_rhs[i]
+                b_block = BlockVector(zeros(Tv, total_size), block_sizes)
+                b_block[Block(1)] = b_unrestricted
 
-                    end
-                end
+                u_unrestricted = linsolve.u
+                u_block = BlockVector(zeros(Tv, total_size), block_sizes)
+                u_block[Block(1)] = u_unrestricted
+
+                for i in eachindex(PD.restrictions)
+                    A_block[Block(1, i + 1)] = restriction_matrices[i]
+                    A_block[Block(i + 1, 1)] = transpose(restriction_matrices[i])
+                    b_block[Block(i + 1)] = restriction_rhs[i]
 
-                @timeit timer "convert" begin
+                end
 
-                    linsolve.b = Vector(b_block) # convert to dense vector
-                    linsolve.u = Vector(u_block) # convert to dense vector
+                linsolve.b = Vector(b_block) # convert to dense vector
+                linsolve.u = Vector(u_block) # convert to dense vector
 
-                    # linsolve.A = sparse(A_block) # convert to CSC Matrix is very slow https://github.com/JuliaArrays/BlockArrays.jl/issues/78
-                    # do it manually:
-                    A_flat = spzeros(size(A_block))
+                # linsolve.A = sparse(A_block) # convert to CSC Matrix is very slow https://github.com/JuliaArrays/BlockArrays.jl/issues/78
+                # do it manually:
+                A_flat = spzeros(size(A_block))
 
-                    lasts_i = [0, axes(A_block)[1].lasts...] # add zero
-                    lasts_j = [0, axes(A_block)[2].lasts...] # add zero
+                lasts_row = [0, axes(A_block)[1].lasts...] # add leading zero
+                lasts_col = [0, axes(A_block)[2].lasts...] # add leading zero
 
-                    (ni, nj) = size(blocks(A_block))
+                (n_row, n_col) = size(blocks(A_block))
 
-                    for i in 1:ni, j in 1:nj
-                        range_i = (lasts_i[i] + 1):lasts_i[i + 1]
-                        range_j = (lasts_j[j] + 1):lasts_j[j + 1]
+                for i in 1:n_row, j in 1:n_col
+                    range_row = (lasts_row[i] + 1):lasts_row[i + 1]
+                    range_col = (lasts_col[j] + 1):lasts_col[j + 1]
 
-                        # write each block directly in the resulting matrix
-                        A_flat[range_i, range_j] = A_block[Block(i, j)]
-                    end
-                    linsolve.A = A_flat
+                    # write each block directly in the resulting matrix
+                    A_flat[range_row, range_col] = A_block[Block(i, j)]
                 end
+                linsolve.A = A_flat
             end
         end