32 const std::vector<double>& local_x_data,
33 const std::vector<double>& local_x_prev_data,
34 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
35 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
41 "The number of specified epsilons ({:d}) and the number of local "
42 "d.o.f.s ({:d}) do not match, i.e., the latter is not divisible by "
48 static_cast<Eigen::MatrixXd::Index
>(local_x_data.size());
51 MathLib::toVector<Eigen::VectorXd>(local_x_data, num_r_c);
52 auto const local_x_prev =
53 MathLib::toVector<Eigen::VectorXd>(local_x_prev_data, num_r_c);
54 Eigen::VectorXd
const local_xdot = (local_x - local_x_prev) / dt;
60 auto const num_dofs_per_component =
71 for (Eigen::MatrixXd::Index i = 0; i < num_r_c; ++i)
74 auto const component = i / num_dofs_per_component;
79 local_x_prev_data, local_M_data, local_K_data,
89 if (!local_M_data.empty())
91 auto const local_M_p =
93 auto const local_M_m =
95 local_Jac.col(i).noalias() +=
97 (local_M_p - local_M_m) * local_xdot / (2.0 * eps);
101 if (!local_K_data.empty())
103 auto const local_K_p =
105 auto const local_K_m =
107 local_Jac.col(i).noalias() +=
109 (local_K_p - local_K_m) * local_x / (2.0 * eps);
110 local_K_data.clear();
113 if (!local_b_data.empty())
115 auto const local_b_p =
116 MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
117 auto const local_b_m =
119 local_Jac.col(i).noalias() -=
121 (local_b_p - local_b_m) / (2.0 * eps);
122 local_b_data.clear();
128 local_assembler.
assemble(t, dt, local_x_data, local_x_prev_data,
129 local_M_data, local_K_data, local_b_data);
132 if (!local_M_data.empty())
135 local_Jac.noalias() += local_M / dt;
137 if (!local_K_data.empty())
140 local_Jac.noalias() += local_K;
148 if (!local_b_data.empty())
150 return MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
152 return MathLib::createZeroedVector<Eigen::VectorXd>(local_b_data,
156 if (!local_M_data.empty())
160 local_M_data.clear();
162 if (!local_K_data.empty())
166 local_K_data.clear();
179 "relative_epsilons");
182 "component_magnitudes");
184 if (!!rel_eps != !!comp_mag)
187 "Either both or none of <relative_epsilons> and "
188 "<component_magnitudes> have to be specified.");
191 std::vector<double> abs_eps;
195 if (rel_eps->size() != comp_mag->size())
198 "The numbers of components of <relative_epsilons> and "
199 "<component_magnitudes> do not match.");
202 abs_eps.resize(rel_eps->size());
203 for (std::size_t i = 0; i < rel_eps->size(); ++i)
205 abs_eps[i] = (*rel_eps)[i] * (*comp_mag)[i];
212 abs_eps.emplace_back(1e-8);
215 return std::make_unique<CentralDifferencesJacobianAssembler>(
void assembleWithJacobian(LocalAssemblerInterface &local_assembler, double const t, double const dt, std::vector< double > const &local_x_data, std::vector< double > const &local_x_prev_data, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
virtual void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)