Merge pull request #85 from WIAS-PDELib/fix/residual

Fix fixed_dofs of CoupledDofsRestriction

Author: pjaap

examples/Example212_PeriodicBoundary2D.jl

@@ -163,7 +163,8 @@ function main(;
         end
     end
 
-    sol = solve(PD, FES)
+    sol, SC = solve(PD, FES, return_config = true)
+    residual(SC) < 1.0e-10 || error("Residual is not zero!")
 
     plt = GridVisualizer(; Plotter, size = (1300, 800))
 

examples/Example312_PeriodicBoundary3D.jl

@@ -179,7 +179,8 @@ function main(;
     end
 
     ## solve
-    sol = solve(PD, FES; kwargs...)
+    sol, SC = solve(PD, FES; return_config = true, kwargs...)
+    residual(SC) < 1.0e-10 || error("Residual is not zero!")
 
     displace_mesh!(xgrid, sol[u])
     plt = plot([grid(u)], sol; Plotter, do_vector_plots = false, width = 1200, height = 800, title = "displaced mesh", scene3d = :LScene)

src/common_restrictions/coupled_dofs_restriction.jl

@@ -33,7 +33,10 @@ function assemble!(R::CoupledDofsRestriction, sol, SC; kwargs...)
 
     R.parameters[:matrix] = B
     R.parameters[:rhs] = Zeros(length(unique_cols))
-    R.parameters[:fixed_dofs] = unique_cols
+
+    # fixed dofs are all active rows of B
+    I, _, _ = findnz(B)
+    R.parameters[:fixed_dofs] = unique(I)
 
     return nothing
 end