diff --git a/src/core/fiber.cpp b/src/core/fiber.cpp
index ea28f6e..4faece7 100644
--- a/src/core/fiber.cpp
+++ b/src/core/fiber.cpp
@@ -233,16 +233,7 @@ void Fiber::update_RHS(double dt, MatrixRef &flow, MatrixRef &f_external) {
         // FIXME: Does this imply always having velocity boundary conditions?
         xs_vT(bc_start_i + 3) = flow.col(minus_node).dot(xs_.col(minus_node));
 
-        RHS_ += vT_in - xs_vT;
-
-        RHS_.segment(0 * np, np) += flow.row(0);
-        RHS_.segment(1 * np, np) += flow.row(1);
-        RHS_.segment(2 * np, np) += flow.row(2);
-
-        RHS_.segment(3 * np, np) += (xs_x.transpose() * (flow.row(0) * D_1.matrix()).array() +
-                                     xs_y.transpose() * (flow.row(1) * D_1.matrix()).array() +
-                                     xs_z.transpose() * (flow.row(2) * D_1.matrix()).array())
-                                        .matrix();
+        RHS_ += -vT + vT_in - xs_vT;
     }
     if (f_external.size()) {
         ArrayXXd fs = f_external * D_1;
diff --git a/src/core/system.cpp b/src/core/system.cpp
index 6124e76..32384fb 100644
--- a/src/core/system.cpp
+++ b/src/core/system.cpp
@@ -274,12 +274,13 @@ void prep_state_for_solver() {
     CVectorMap shell_velocities_flat(v_all.data() + 3 * fib_node_count, 3 * shell_node_count);
     shell_->solution_vec_ = shell_->M_inv_ * shell_velocities_flat;
     v_fibers_ = shell_->flow(fc_.get_local_node_positions(), shell_->solution_vec_, params_.eta);
-    v_fibers_ += v_all.block(0, 0, 3, fib_node_count);
 
     MatrixXd external_force_fibers = MatrixXd::Zero(3, fib_node_count);
     fc_.update_RHS(properties.dt, v_fibers_, external_force_fibers);
     fc_.update_boundary_conditions(*shell_);
     fc_.apply_bc_rectangular(properties.dt, v_fibers_, external_force_fibers);
+
+    v_fibers_ += v_all.block(0, 0, 3, fib_node_count);
 }