combinedofs should not reassemble matrix or rhs if that is set constant (#57)

* combinedofs should not reassemble matrix or rhs if that is set constant

* changed test for Example252 to some more stable case (otherwise fails only on my machine and not caused by the changes of this PR)

* apply_penalties! of InterpolateBoundaryData now also recognizes assemblr_rhs and assemble_matrix switches

Author: chmerdon

examples/Example252_NSEPlanarLatticeFlow.jl

@@ -127,8 +127,8 @@ end
 
 generateplots = ExtendableFEM.default_generateplots(Example252_NSEPlanarLatticeFlow, "example252.png") #hide
 function runtests() #hide
-    L2errorU, plt = main(; nrefs = 3) #hide
-    @test L2errorU ≈ 0.11892169556349004 #hide
+    L2errorU, plt = main(; nrefs = 4) #hide
+    @test L2errorU ≈ 0.010616923333947861 #hide
     return nothing #hide
 end #hide
 end # module

src/common_operators/combinedofs.jl

@@ -85,11 +85,11 @@ function build_assembler!(CD::CombineDofs{UT, CT}, FE::Array{<:FEVectorBlock, 1}
             @info ".... coupling $(length(coupling_matrix.nzval)) dofs"
         end
         function assemble!(A::AbstractSparseArray{T}, b::AbstractVector{T}, assemble_matrix::Bool, assemble_rhs::Bool, kwargs...) where {T}
-            if assemble_matrix
 
-                # transpose the matrix once for efficient row access
-                transposed_coupling_matrix = sparse(transpose(coupling_matrix))
+            # transpose the matrix once for efficient row access
+            transposed_coupling_matrix = sparse(transpose(coupling_matrix))
 
+            if assemble_matrix
                 # go through each coupled dof and update the FE adjacency info
                 # from the constrained dofs here
 

src/common_operators/interpolateboundarydata_operator.jl

@@ -188,7 +188,7 @@ apply_penalties!(A, b, sol, O::InterpolateBoundaryData{UT}, SC::SolverConfigurat
 modifies the linear system A|b such that the boundary dofs are penalized and attain the
 correct values from the last assemble! call of O. Also applies the correct values to sol.
 """
-function apply_penalties!(A, b, sol, O::InterpolateBoundaryData{UT}, SC::SolverConfiguration; kwargs...) where {UT}
+function apply_penalties!(A, b, sol, O::InterpolateBoundaryData{UT}, SC::SolverConfiguration; assemble_matrix = true, assemble_rhs = true, assemble_sol = true, kwargs...) where {UT}
     time = @elapsed begin
         if UT <: Integer
             ind = O.u
@@ -201,16 +201,22 @@ function apply_penalties!(A, b, sol, O::InterpolateBoundaryData{UT}, SC::SolverC
         bddata = O.bddata
         bdofs = O.bdofs
         penalty = O.parameters[:penalty]
-        AE = A.entries
-        BE = b.entries
-        for dof in bdofs
-            AE[dof, dof] = penalty
+        if assemble_matrix
+            AE = A.entries
+            for dof in bdofs
+                AE[dof, dof] = penalty
+            end
+            flush!(AE)
+        end
+        if assemble_rhs
+            BE = b.entries
+            for dof in bdofs
+                BE[dof] = penalty * bddata.entries[dof - offset]
+            end
         end
-        flush!(AE)
-        for dof in bdofs
-            BE[dof] = penalty * bddata.entries[dof - offset]
+        if assemble_sol
+            apply!(sol[ind_sol], O; offset = offset)
         end
-        apply!(sol[ind_sol], O; offset = offset)
     end
     return if O.parameters[:verbosity] > 1
         @info "$(O.parameters[:name]) : applying penalties took $time s"

src/plots.jl

@@ -19,7 +19,7 @@ function plot!(p::GridVisualizer, ops, sol; kwargs...)
 Plots the operator evaluations ops of blocks in sol into the GridVisualizer.
 
 """
-function plot!(p::GridVisualizer, ops, sol; rasterpoints = 10, keep = [], ncols = size(p.subplots, 2), do_abs = true, do_vector_plots = true, title_add = "", kwargs...)
+function plot!(p::GridVisualizer, ops, sol; rasterpoints = 10, linewidth = 1, keep = [], ncols = size(p.subplots, 2), do_abs = true, do_vector_plots = true, title_add = "", kwargs...)
     col, row, id = 0, 1, 0
     for op in ops
         col += 1
@@ -35,9 +35,9 @@ function plot!(p::GridVisualizer, ops, sol; rasterpoints = 10, keep = [], ncols
             end
         end
         if op[2] == "grid"
-            gridplot!(p[row, col], sol[op[1]].FES.xgrid; kwargs...)
+            gridplot!(p[row, col], sol[op[1]].FES.xgrid; linewidth = linewidth, kwargs...)
         elseif op[2] == "dofgrid"
-            gridplot!(p[row, col], sol[op[1]].FES.dofgrid; kwargs...)
+            gridplot!(p[row, col], sol[op[1]].FES.dofgrid; linewidth = linewidth, kwargs...)
         else
             ncomponents = get_ncomponents(sol[op[1]])
             edim = size(sol[op[1]].FES.xgrid[Coordinates], 1)
@@ -96,7 +96,7 @@ function plot(ops, sol; add = 0, Plotter = nothing, ncols = min(2, length(ops) +
     if height == 0
         height = width / ncols * nrows
     end
-    p = GridVisualizer(; Plotter = Plotter, layout = (nrows, ncols), clear = true, resolution = (width, height))
+    p = GridVisualizer(; Plotter = Plotter, layout = (nrows, ncols), clear = true, size = (width, height))
     return plot!(p, ops, sol; do_abs = do_abs, kwargs...)
 end
 

src/solvers.jl

@@ -191,7 +191,7 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{
 
                 ## penalize fixed dofs
                 time_assembly += @elapsed for op in PD.operators
-                    allocs_assembly += @allocated apply_penalties!(A, b, sol, op, SC; kwargs...)
+                    allocs_assembly += @allocated apply_penalties!(A, b, sol, op, SC; assemble_matrix = !SC.parameters[:initialized] || !SC.parameters[:constant_matrix], assemble_rhs = !SC.parameters[:initialized] || !SC.parameters[:constant_rhs], kwargs...)
                 end
                 flush!(A.entries)
                 # end