81{
82 auto const& eps_m = std::get<MPL::SymmetricTensor<DisplacementDim>>(
83 variable_array.mechanical_strain);
84 auto const& eps_m_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
85 variable_array_prev.mechanical_strain);
86 auto const& sigma_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
87 variable_array_prev.stress);
88
90
91 assert(dynamic_cast<MaterialStateVariables const*>(
92 &material_state_variables) != nullptr);
93 MaterialStateVariables state(
94 static_cast<MaterialStateVariables const&>(material_state_variables));
95 state.setInitialConditions();
96
97 auto local_lubby2_properties =
98 detail::LocalLubby2Properties<DisplacementDim>{t, x,
_mp};
99
100
101 auto const& P_dev = Invariants::deviatoric_projection;
104
105
106 KelvinVector sigd_j = 2.0 * (eps_m_d_i - state.eps_M_t - state.eps_K_t);
107
108 KelvinVector sigd_t = P_dev * sigma_prev / local_lubby2_properties.GM0;
109
110
111 double sig_eff = Invariants::equivalentStress(sigd_j);
112 local_lubby2_properties.update(sig_eff);
113
114 using LocalJacobianMatrix =
116 Eigen::RowMajor>;
117
118
119
120 Eigen::FullPivLU<LocalJacobianMatrix> linear_solver;
121
122
123
124
125
126
127
128
129
130
131
132
133 LocalJacobianMatrix K_loc;
134 {
135 using LocalResidualVector =
136 Eigen::Matrix<double, KelvinVectorSize * 3, 1>;
137
138 auto const update_residual = [&](LocalResidualVector& residual)
139 {
141 state.eps_K_j, state.eps_K_t,
142 state.eps_M_j, state.eps_M_t, residual,
143 local_lubby2_properties);
144 };
145
146 auto const update_jacobian = [&](LocalJacobianMatrix& jacobian)
147 {
149 t, x, dt, jacobian, sig_eff, sigd_j, state.eps_K_j,
150 local_lubby2_properties);
151 };
152
153 auto const update_solution = [&](LocalResidualVector const& increment)
154 {
155
156 sigd_j.noalias() += increment.template segment<KelvinVectorSize>(
158 state.eps_K_j.noalias() +=
160 1);
161 state.eps_M_j.noalias() +=
163 2);
164
165
168 local_lubby2_properties.update(sig_eff);
169 };
170
172 linear_solver, update_jacobian, update_residual, update_solution,
174
175 auto const success_iterations = newton_solver.solve(K_loc);
176
177 if (!success_iterations)
178 {
179 return {};
180 }
181
182
183
184
185 if (*success_iterations == 0)
186 {
187 linear_solver.compute(K_loc);
188 }
189 }
190
193 local_lubby2_properties.KM0,
194 linear_solver);
195
196
197 double const delta_eps_m_trace = Invariants::trace(eps_m - eps_m_prev);
198 double const sigma_trace_prev = Invariants::trace(sigma_prev);
200 local_lubby2_properties.GM0 * sigd_j +
201 (local_lubby2_properties.KM0 * delta_eps_m_trace +
202 sigma_trace_prev / 3.) *
203 Invariants::identity2;
204 return {std::make_tuple(
205 sigma,
206 std::unique_ptr<
207 typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
208 new MaterialStateVariables{state}},
210}
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
void calculateJacobianBurgers(double const t, ParameterLib::SpatialPosition const &x, double const dt, JacobianMatrix &Jac, double s_eff, const KelvinVector &sig_i, const KelvinVector &eps_K_i, detail::LocalLubby2Properties< DisplacementDim > const &properties) const
Calculates the 18x18 Jacobian.
void calculateResidualBurgers(double const dt, const KelvinVector &strain_curr, const KelvinVector &strain_t, const KelvinVector &stress_curr, const KelvinVector &stress_t, const KelvinVector &strain_Kel_curr, const KelvinVector &strain_Kel_t, const KelvinVector &strain_Max_curr, const KelvinVector &strain_Max_t, ResidualVector &res, detail::LocalLubby2Properties< DisplacementDim > const &properties) const
Calculates the 18x1 residual vector.
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > tangentStiffnessA(double const GM0, double const KM0, LinearSolver const &linear_solver)