diff --git a/include/fiber.hpp b/include/fiber.hpp index a708b5b..82672a1 100644 --- a/include/fiber.hpp +++ b/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; diff --git a/src/core/fiber.cpp b/src/core/fiber.cpp index 4faece7..875e088 100644 --- a/src/core/fiber.cpp +++ b/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 /// diff --git a/src/core/fiber_container.cpp b/src/core/fiber_container.cpp index dec0324..bb322c6 100644 --- a/src/core/fiber_container.cpp +++ b/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; diff --git a/src/core/system.cpp b/src/core/system.cpp index 32384fb..22040a7 100644 --- a/src/core/system.cpp +++ b/src/core/system.cpp @@ -32,7 +32,6 @@ BackgroundSource bs_; ///< Background flow std::unique_ptr 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, external_force_fibers); }