Fix ranges of the solution vector

Author: pjaap

src/solvers.jl

@@ -361,7 +361,19 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{
                 push!(stats[:matrix_nnz], nnz(linsolve.A))
                 Δx = LinearSolve.solve!(linsolve)
 
-                x = sol.entries - Δx.u
+                # x = sol.entries - Δx.u for free dofs or partial solutions
+                if length(freedofs) > 0
+                    x = sol.entries[freedofs] - Δx.u
+                else
+                    x = zero(Δx)
+                    offset = 0
+                    for u in unknowns
+                        ndofs_u = length(view(sol[u]))
+                        x_range = (offset + 1):(offset + ndofs_u)
+                        x[x_range] .= view(sol[u]) .- view(Δx, x_range)
+                        offset += ndofs_u
+                    end
+                end
 
                 ## check linear residual with full matrix
                 if length(freedofs) > 0
@@ -656,7 +668,15 @@ function iterate_until_stationarity(
                         ## solve
                         Δx = LinearSolve.solve!(linsolve)
 
-                        x = sol.entries - Δx.u
+                        # x = sol.entries - Δx.u ... in the entry ranges of the present unknowns
+                        x = zero(Δx.u)
+                        offset = 0
+                        for u in unknowns[p]
+                            ndofs_u = length(view(sol[u]))
+                            x_range = (offset + 1):(offset + ndofs_u)
+                            x[x_range] .= view(sol[u]) .- view(Δx.u, x_range)
+                            offset += ndofs_u
+                        end
 
                         fill!(residual.entries, 0)
                         mul!(residual.entries, A.entries.cscmatrix, x)