21 std::vector<double>&& absolute_epsilons)
22 : _absolute_epsilons(std::move(absolute_epsilons))
26 OGS_FATAL(
"No values for the absolute epsilons have been given.");
32 const std::vector<double>& local_x_data,
33 const std::vector<double>& local_xdot_data,
const double dxdot_dx,
34 const double dx_dx, std::vector<double>& local_M_data,
35 std::vector<double>& local_K_data, std::vector<double>& local_b_data,
36 std::vector<double>& local_Jac_data)
42 "The number of specified epsilons ({:d}) and the number of local "
43 "d.o.f.s ({:d}) do not match, i.e., the latter is not divisible by "
49 static_cast<Eigen::MatrixXd::Index
>(local_x_data.size());
52 MathLib::toVector<Eigen::VectorXd>(local_x_data, num_r_c);
53 auto const local_xdot =
54 MathLib::toVector<Eigen::VectorXd>(local_xdot_data, num_r_c);
61 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;
81 local_K_data, local_b_data);
92 if (!local_M_data.empty())
94 auto const local_M_p =
96 auto const local_M_m =
98 local_Jac.col(i).noalias() +=
100 (local_M_p - local_M_m) * local_xdot / (2.0 * eps);
101 local_M_data.clear();
104 if (!local_K_data.empty())
106 auto const local_K_p =
108 auto const local_K_m =
110 local_Jac.col(i).noalias() +=
112 (local_K_p - local_K_m) * local_x / (2.0 * eps);
113 local_K_data.clear();
116 if (!local_b_data.empty())
118 auto const local_b_p =
119 MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
120 auto const local_b_m =
122 local_Jac.col(i).noalias() -=
124 (local_b_p - local_b_m) / (2.0 * eps);
125 local_b_data.clear();
131 local_assembler.
assemble(t, dt, local_x_data, local_xdot_data, local_M_data,
132 local_K_data, local_b_data);
135 if (dxdot_dx != 0.0 && !local_M_data.empty())
138 local_Jac.noalias() += local_M * dxdot_dx;
140 if (dx_dx != 0.0 && !local_K_data.empty())
143 local_Jac.noalias() += local_K * dx_dx;
147 std::unique_ptr<CentralDifferencesJacobianAssembler>
156 "relative_epsilons");
159 "component_magnitudes");
161 if (!!rel_eps != !!comp_mag)
164 "Either both or none of <relative_epsilons> and "
165 "<component_magnitudes> have to be specified.");
168 std::vector<double> abs_eps;
172 if (rel_eps->size() != comp_mag->size())
175 "The numbers of components of <relative_epsilons> and "
176 "<component_magnitudes> do not match.");
179 abs_eps.resize(rel_eps->size());
180 for (std::size_t i = 0; i < rel_eps->size(); ++i)
182 abs_eps[i] = (*rel_eps)[i] * (*comp_mag)[i];
189 abs_eps.emplace_back(1e-8);
192 return std::make_unique<CentralDifferencesJacobianAssembler>(
void checkConfigParameter(std::string const ¶m, T const &value) const
std::optional< T > getConfigParameterOptional(std::string const ¶m) const
std::vector< double > _local_b_data
std::vector< double > _local_K_data
CentralDifferencesJacobianAssembler(std::vector< double > &&absolute_epsilons)
std::vector< double > _local_x_perturbed_data
std::vector< double > const _absolute_epsilons
std::vector< double > _local_xdot_perturbed_data
void assembleWithJacobian(LocalAssemblerInterface &local_assembler, double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_xdot, const double dxdot_dx, const double dx_dx, 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
std::vector< double > _local_M_data
virtual void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_xdot, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
std::unique_ptr< CentralDifferencesJacobianAssembler > createCentralDifferencesJacobianAssembler(BaseLib::ConfigTree const &config)