diff --git a/docs/src/literate/computational_homogenization.jl b/docs/src/literate/computational_homogenization.jl
index 09872f7b45..baa486f09c 100644
--- a/docs/src/literate/computational_homogenization.jl
+++ b/docs/src/literate/computational_homogenization.jl
@@ -368,7 +368,6 @@ rhsdata = (
 )
 
 apply!(K.dirichlet, ch.dirichlet)
-Kp = copy(K.periodic) #hide
 apply!(K.periodic,  ch.periodic)
 
 # We can now solve the problem(s). Note that we only use `apply_rhs!` in the loops below.
@@ -388,14 +387,9 @@ for i in 1:size(rhs.dirichlet, 2)
     push!(u.dirichlet, u_i)                            # Save the solution vector
 end
 
-rhs_p = copy(rhs.periodic) #hide
 for i in 1:size(rhs.periodic, 2)
     rhs_i = @view rhs.periodic[:, i]                   # Extract this RHS
     apply_rhs!(rhsdata.periodic, rhs_i, ch.periodic)   # Apply BC
-    rhs_i = @view rhs_p[:, i] #hide
-    Kpp = copy(Kp) #hide
-    apply!(Kpp, rhs_i, ch.periodic) #hide
-    copy!(K.periodic, Kpp) #hide
     u_i = cholesky(Symmetric(K.periodic)) \ rhs_i      # Solve
     apply!(u_i, ch.periodic)                           # Apply BC on the solution
     push!(u.periodic, u_i)                             # Save the solution vector
diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl
index eeae5fb52b..fc941e1b42 100644
--- a/src/Dofs/ConstraintHandler.jl
+++ b/src/Dofs/ConstraintHandler.jl
@@ -111,22 +111,31 @@ See also: [`get_rhs_data`](@ref).
 """
 function apply_rhs!(data::RHSData, f::AbstractVector, ch::ConstraintHandler, applyzero::Bool=false)
     K = data.constrained_columns
-    @assert length(f) == 0 || length(f) == size(K, 1)
-    @boundscheck length(f) == 0 || checkbounds(f, ch.prescribed_dofs)
-
+    @assert length(f) == size(K, 1)
+    @boundscheck checkbounds(f, ch.prescribed_dofs)
     m = data.m
+
+    # TODO: Can the loops be combined or does the order matter?
     @inbounds for i in 1:length(ch.inhomogeneities)
-        d = ch.prescribed_dofs[i]
         v = ch.inhomogeneities[i]
         if !applyzero && v != 0
             for j in nzrange(K, i)
                 f[K.rowval[j]] -= v * K.nzval[j]
             end
         end
-        if length(f) != 0
-            vz = applyzero ? zero(eltype(f)) : v
-            f[d] = vz * m
+    end
+    @inbounds for ac in ch.acs
+        global_dof = ac.constrained_dof
+        for (d, v) in ac.entries
+            f[d] += f[global_dof] * v
         end
+        f[global_dof] = 0
+    end
+    @inbounds for i in 1:length(ch.inhomogeneities)
+        d = ch.prescribed_dofs[i]
+        v = ch.inhomogeneities[i]
+        vz = applyzero ? zero(eltype(f)) : v
+        f[d] = vz * m
     end
 end
 
@@ -679,6 +688,7 @@ function _condense!(K::SparseMatrixCSC, f::AbstractVector, acs::Vector{AffineCon
                 for (d,v) in ac.entries
                     f[d] += f[col] * v
                 end
+                @assert ac.constrained_dof == col
                 f[ac.constrained_dof] = 0.0
             end
         end
diff --git a/test/test_constraints.jl b/test/test_constraints.jl
index 3151cecdfd..079f6c3f12 100644
--- a/test/test_constraints.jl
+++ b/test/test_constraints.jl
@@ -150,6 +150,7 @@ end
         for cell in CellIterator(dh)
             K[celldofs(cell), celldofs(cell)] += 2.0 * [1 -1; -1 1]
         end
+        copies = (K = copy(K), f1 = copy(f), f2 = copy(f))
 
         # Solve by actually condensing the matrix
         ff  = C' * (f - K * g)
@@ -162,7 +163,15 @@ end
         a = K \ f
         apply!(a, ch)
 
-        @test a ≈ aa
+        # Solve with extracted RHS data
+        rhs = get_rhs_data(ch, copies.K)
+        apply!(copies.K, ch)
+        apply_rhs!(rhs, copies.f1, ch)
+        a_rhs1 = apply!(copies.K \ copies.f1, ch)
+        apply_rhs!(rhs, copies.f2, ch)
+        a_rhs2 = apply!(copies.K \ copies.f2, ch)
+
+        @test a ≈ aa ≈ a_rhs1 ≈ a_rhs2
     end
 
 end