OGS
ProcessLib::CentralDifferencesJacobianAssembler Class Referencefinal

Detailed Description

Assembles the Jacobian matrix using central differences.

Definition at line 18 of file CentralDifferencesJacobianAssembler.h.

#include <CentralDifferencesJacobianAssembler.h>

Inheritance diagram for ProcessLib::CentralDifferencesJacobianAssembler:
[legend]
Collaboration diagram for ProcessLib::CentralDifferencesJacobianAssembler:
[legend]

Public Member Functions

 CentralDifferencesJacobianAssembler (std::vector< double > &&absolute_epsilons)
void assembleWithJacobian (std::size_t const mesh_item_id, 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_b_data, std::vector< double > &local_Jac_data) override
std::unique_ptr< AbstractJacobianAssembler > copy () const override
Public Member Functions inherited from ProcessLib::NumericalJacobianAssembler
 NumericalJacobianAssembler (std::vector< double > &&absolute_epsilons)
bool needsPicardAssembly () const override
void checkPerturbationSize (int const max_non_deformation_dofs_per_node) const override
void setNonDeformationComponentIDs (std::vector< int > const &non_deformation_component_ids) override
void setNonDeformationComponentIDsNoSizeCheck (std::vector< int > const &non_deformation_component_ids) override
Public Member Functions inherited from ProcessLib::AbstractJacobianAssembler
virtual void assembleWithJacobianForStaggeredScheme (LocalAssemblerInterface &, double const, double const, Eigen::VectorXd const &, Eigen::VectorXd const &, int const, std::vector< double > &, std::vector< double > &)
virtual ~AbstractJacobianAssembler ()=default
virtual void preIteration (const unsigned)

Private Attributes

std::vector< double > _local_M_data
std::vector< double > _local_K_data
std::vector< double > _local_b_data
std::vector< double > _local_x_perturbed_data

Additional Inherited Members

Protected Member Functions inherited from ProcessLib::NumericalJacobianAssembler
auto getVariableComponentEpsilonsView () const
Protected Attributes inherited from ProcessLib::NumericalJacobianAssembler
std::vector< double > const absolute_epsilons_
std::vector< int > non_deformation_component_ids_

Constructor & Destructor Documentation

◆ CentralDifferencesJacobianAssembler()

ProcessLib::CentralDifferencesJacobianAssembler::CentralDifferencesJacobianAssembler ( std::vector< double > && absolute_epsilons)
explicit

Constructs a new instance.

Parameters
absolute_epsilonsperturbations of the variable components used for evaluating the finite differences.
Note
The size of absolute_epsilons defines the "number of components" of the local solution vector (This is not the number of elements of the vector!). Therefore the size of the local solution vector must be divisible by the size of absolute_epsilons. This is the only consistency check performed. It is not checked whether said "number of components" is sensible. E.g., one could pass one epsilon per node, which would be valid but would not make sense at all.

Definition at line 13 of file CentralDifferencesJacobianAssembler.cpp.

15 : NumericalJacobianAssembler(std::move(absolute_epsilons))
16{
17}
ProcessLib::NumericalJacobianAssembler::NumericalJacobianAssembler
NumericalJacobianAssembler(std::vector< double > &&absolute_epsilons)
Definition NumericalJacobianAssembler.h:17

References ProcessLib::NumericalJacobianAssembler::NumericalJacobianAssembler().

Member Function Documentation

◆ assembleWithJacobian()

void ProcessLib::CentralDifferencesJacobianAssembler::assembleWithJacobian ( std::size_t const mesh_item_id,
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_b_data,
std::vector< double > & local_Jac_data )
overridevirtual

Assembles the Jacobian, the matrices \(M\) and \(K\), and the vector \(b\). For the assembly the assemble() method of the given local_assembler is called several times and the Jacobian is built from finite differences. The number of calls of the assemble() method is \(2N+1\) if \(N\) is the size of local_x_data.

Attention
It is assumed that the local vectors and matrices are ordered by component.

Implements ProcessLib::AbstractJacobianAssembler.

Definition at line 19 of file CentralDifferencesJacobianAssembler.cpp.

25{
26 std::vector<double> local_M_data(local_Jac_data.size());
27 std::vector<double> local_K_data(local_Jac_data.size());
28
29 auto const num_r_c =
30 static_cast<Eigen::MatrixXd::Index>(local_x_data.size());
31
32 auto const local_x =
33 MathLib::toVector<Eigen::VectorXd>(local_x_data, num_r_c);
34 auto const local_x_prev =
35 MathLib::toVector<Eigen::VectorXd>(local_x_prev_data, num_r_c);
36 Eigen::VectorXd const local_xdot = (local_x - local_x_prev) / dt;
37
38 auto local_Jac =
39 MathLib::createZeroedMatrix(local_Jac_data, num_r_c, num_r_c);
40 _local_x_perturbed_data = local_x_data;
41
42 assert(this->non_deformation_component_ids_.size() > 0);
43 auto const num_dofs_per_component =
44 local_x_data.size() / this->non_deformation_component_ids_.size();
45
46 auto const nved = local_assembler.getNumberOfVectorElementsForDeformation();
47
48 // Residual res := M xdot + K x - b
49 // Computing Jac := dres/dx
50 // = M dxdot/dx + dM/dx xdot + K dx/dx + dK/dx x - db/dx
51 // with dxdot/dx = 1/dt and dx/dx = 1
52 // (Note: dM/dx and dK/dx actually have the second and
53 // third index transposed.)
54 // The loop computes the dM/dx, dK/dx and db/dx terms, the rest is computed
55 // afterwards. The loop skips the entries corresponding to the deformation
56 // part of the solution vector if a vector segment size is given by nved.
57 // This is to avoid recomputing the analytic block of the deformation
58 // process.
59 auto const num_purterbated_colums = num_r_c - nved;
60 auto const perturbations = this->getVariableComponentEpsilonsView();
61 for (Eigen::MatrixXd::Index i = 0; i < num_purterbated_colums; ++i)
62 {
63 // assume that local_x_data is ordered by component.
64 auto const component = i / num_dofs_per_component;
65 auto const eps = perturbations[component];
66
67 _local_x_perturbed_data[i] += eps;
68 local_assembler.assemble(t, dt, _local_x_perturbed_data,
69 local_x_prev_data, local_M_data, local_K_data,
70 local_b_data);
71
72 _local_x_perturbed_data[i] = local_x_data[i] - eps;
73 local_assembler.assemble(t, dt, _local_x_perturbed_data,
74 local_x_prev_data, _local_M_data,
75 _local_K_data, _local_b_data);
76
77 _local_x_perturbed_data[i] = local_x_data[i];
78
79 if (!local_M_data.empty())
80 {
81 auto const local_M_p =
82 MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
83 auto const local_M_m =
84 MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
85 local_Jac.col(i).noalias() +=
86 // dM/dxi * x_dot
87 (local_M_p - local_M_m) * local_xdot / (2.0 * eps);
88 local_M_data.clear();
89 _local_M_data.clear();
90 }
91 if (!local_K_data.empty())
92 {
93 auto const local_K_p =
94 MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
95 auto const local_K_m =
96 MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
97 local_Jac.col(i).noalias() +=
98 // dK/dxi * x
99 (local_K_p - local_K_m) * local_x / (2.0 * eps);
100 local_K_data.clear();
101 _local_K_data.clear();
102 }
103 if (!local_b_data.empty())
104 {
105 auto const local_b_p =
106 MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
107 auto const local_b_m =
108 MathLib::toVector<Eigen::VectorXd>(_local_b_data, num_r_c);
109 local_Jac.col(i).noalias() -=
110 // db/dxi
111 (local_b_p - local_b_m) / (2.0 * eps);
112 local_b_data.clear();
113 _local_b_data.clear();
114 }
115 }
116
117 // Assemble with unperturbed local x, i.e. compute M dxdot/dx + K dx/dx =
118 // M/dt + K
119 local_assembler.assemble(t, dt, local_x_data, local_x_prev_data,
120 local_M_data, local_K_data, local_b_data);
121
122 // Compute remaining terms of the Jacobian.
123 if (!local_M_data.empty())
124 {
125 auto local_M = MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
126 local_Jac.noalias() += local_M / dt;
127 }
128 if (!local_K_data.empty())
129 {
130 auto local_K = MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
131 local_Jac.noalias() += local_K;
132 }
133
134 // Move the M and K contributions to the residuum for evaluation of nodal
135 // forces, flow rates, and the like. Cleaning up the M's and K's storage so
136 // it is not accounted for twice.
137 auto b = [&]()
138 {
139 if (!local_b_data.empty())
140 {
141 return MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
142 }
143 return MathLib::createZeroedVector<Eigen::VectorXd>(local_b_data,
144 num_r_c);
145 }();
146
147 if (!local_M_data.empty())
148 {
149 auto M = MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
150 b -= M * local_xdot;
151 local_M_data.clear();
152 }
153 if (!local_K_data.empty())
154 {
155 auto K = MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
156
157 // Note: The deformation segment of \c b is already computed as
158 // int{B^T sigma}dA, which is identical to K_uu * u. Therefore the
159 // corresponding K block is set to zero.
160 if (nved != 0)
161 {
162 auto const dm_start_index = num_purterbated_colums;
163 auto const dm_size = nved;
164 K.block(dm_start_index, dm_start_index, dm_size, dm_size).setZero();
165 }
166
167 b -= K * local_x;
168 local_K_data.clear();
169 }
170}
MathLib::createZeroedVector
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Definition EigenMapTools.h:146
MathLib::toVector
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
Definition EigenMapTools.h:160
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:25
MathLib::toMatrix
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:65

References _local_b_data, _local_K_data, _local_M_data, _local_x_perturbed_data, ProcessLib::LocalAssemblerInterface::assemble(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), ProcessLib::LocalAssemblerInterface::getNumberOfVectorElementsForDeformation(), ProcessLib::NumericalJacobianAssembler::getVariableComponentEpsilonsView(), ProcessLib::NumericalJacobianAssembler::non_deformation_component_ids_, MathLib::toMatrix(), and MathLib::toVector().

◆ copy()

std::unique_ptr< AbstractJacobianAssembler > ProcessLib::CentralDifferencesJacobianAssembler::copy ( ) const
overridevirtual

Implements ProcessLib::AbstractJacobianAssembler.

Definition at line 173 of file CentralDifferencesJacobianAssembler.cpp.

174{
175 return std::make_unique<CentralDifferencesJacobianAssembler>(*this);
176}

Member Data Documentation

◆ _local_b_data

std::vector<double> ProcessLib::CentralDifferencesJacobianAssembler::_local_b_data
private

Definition at line 62 of file CentralDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _local_K_data

std::vector<double> ProcessLib::CentralDifferencesJacobianAssembler::_local_K_data
private

Definition at line 61 of file CentralDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _local_M_data

std::vector<double> ProcessLib::CentralDifferencesJacobianAssembler::_local_M_data
private

Definition at line 60 of file CentralDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _local_x_perturbed_data

std::vector<double> ProcessLib::CentralDifferencesJacobianAssembler::_local_x_perturbed_data
private

Definition at line 63 of file CentralDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

The documentation for this class was generated from the following files: