core: remove mixing due to reapplication of BC

include/fiber.hpp

@@ -135,7 +135,7 @@ class Fiber {
     }
 
 
-    Eigen::VectorXd matvec(VectorRef x, MatrixRef v) const;
+    Eigen::VectorXd matvec(VectorRef x) const;
     void update_preconditioner();
     void update_force_operator();
     void update_boundary_conditions(Periphery &shell);
@@ -217,7 +217,7 @@ class FiberContainer {
     const Eigen::MatrixXd &get_local_node_positions() const { return r_fib_local_; };
     Eigen::VectorXd get_RHS() const;
     Eigen::MatrixXd flow(const MatrixRef &r_trg, const MatrixRef &forces, double eta, bool subtract_self) const;
-    Eigen::VectorXd matvec(VectorRef &x_all, MatrixRef &v_fib) const;
+    Eigen::VectorXd matvec(VectorRef &x_all) const;
     Eigen::MatrixXd apply_fiber_force(VectorRef &x_all) const;
     Eigen::VectorXd apply_preconditioner(VectorRef &x_all) const;
 

src/core/fiber.cpp

@@ -209,31 +209,14 @@ void Fiber::update_RHS(double dt, MatrixRef &flow, MatrixRef &f_external) {
     RHS_.segment(3 * np, np) = -penalty_param_ * VectorXd::Ones(np);
 
     if (flow.size()) {
-        const int bc_start_i = 4 * np - 14;
-        MatrixXd xsDs = (D_1.array().colwise() * xs_.row(0).transpose().array()).transpose();
-        MatrixXd ysDs = (D_1.array().colwise() * xs_.row(1).transpose().array()).transpose();
-        MatrixXd zsDs = (D_1.array().colwise() * xs_.row(2).transpose().array()).transpose();
-
-        VectorXd vT = VectorXd(np * 4);
-        VectorXd v_x = flow.row(0).transpose();
-        VectorXd v_y = flow.row(1).transpose();
-        VectorXd v_z = flow.row(2).transpose();
-
-        vT.segment(0 * np, np) = v_x;
-        vT.segment(1 * np, np) = v_y;
-        vT.segment(2 * np, np) = v_z;
-        vT.segment(3 * np, np) = xsDs * v_x + ysDs * v_y + zsDs * v_z;
-
-        VectorXd vT_in = VectorXd::Zero(4 * np);
-        vT_in.segment(0, bc_start_i) = mats.P_downsample_bc * vT;
-
-        VectorXd xs_vT = VectorXd::Zero(4 * np); // from attachments
-        const int minus_node = 0;
-
-        // 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 + 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();
     }
     if (f_external.size()) {
         ArrayXXd fs = f_external * D_1;
@@ -278,9 +261,7 @@ void Fiber::update_RHS(double dt, MatrixRef &flow, MatrixRef &f_external) {
     }
 }
 
-VectorXd Fiber::matvec(VectorRef x, MatrixRef v) const {
-    return A_ * x;
-}
+VectorXd Fiber::matvec(VectorRef x) const { return A_ * x; }
 
 /// @brief Calculate the force operator cache variable
 ///

src/core/fiber_container.cpp

@@ -71,14 +71,14 @@ VectorXd FiberContainer::apply_preconditioner(VectorRef &x_all) const {
     return y;
 }
 
-VectorXd FiberContainer::matvec(VectorRef &x_all, MatrixRef &v_fib) const {
+VectorXd FiberContainer::matvec(VectorRef &x_all) const {
     VectorXd res = VectorXd::Zero(get_local_solution_size());
 
     size_t offset = 0;
     int i_fib = 0;
     for (const auto &fib : *this) {
         const int np = fib.n_nodes_;
-        res.segment(offset, 4 * np) = fib.matvec(x_all.segment(offset, 4 * np), v_fib.block(0, i_fib, 3, np));
+        res.segment(offset, 4 * np) = fib.matvec(x_all.segment(offset, 4 * np));
 
         i_fib++;
         offset += 4 * np;

src/core/system.cpp

@@ -32,7 +32,6 @@ BackgroundSource bs_;      ///< Background flow
 
 std::unique_ptr<Periphery> shell_; ///< Periphery
 Eigen::VectorXd curr_solution_;    ///< Current MPI-rank local solution vector
-Eigen::MatrixXd v_fibers_;    ///< Current MPI-rank local velocity on fibers
 
 FiberContainer fc_bak_;   ///< Copy of fibers for timestep reversion
 int rank_;                ///< MPI rank
@@ -196,7 +195,7 @@ Eigen::VectorXd apply_matvec(VectorRef &x) {
     const FiberContainer &fc = fc_;
 
     auto [x_fibers] = get_solution_maps(x.data());
-    return fc.matvec(x_fibers, v_fibers_);
+    return fc.matvec(x_fibers);
 }
 
 /// @brief Evaluate the velocity at a list of target points
@@ -271,16 +270,20 @@ void prep_state_for_solver() {
     v_all += psc_.flow(r_all, params_.eta, properties.time);
     v_all += bs_.flow(r_all, params_.eta);
 
-    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);
+    MatrixXd v_shell2fib;
+    if (shell_->is_active()) {
+        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_shell2fib = shell_->flow(fc_.get_local_node_positions(), shell_->solution_vec_, params_.eta);
+    }
+    else {
+        v_shell2fib = MatrixXd::Zero(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_RHS(properties.dt, v_shell2fib, 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);
+    fc_.apply_bc_rectangular(properties.dt, v_shell2fib, force_fibers);
 }