OGS 6.2.1-97-g73d1aeda3
Ehlers.cpp
Go to the documentation of this file.
1 
10 #include "Ehlers.h"
11 #include <boost/math/special_functions/pow.hpp>
13 
14 #include "LinearElasticIsotropic.h"
15 
40 namespace MaterialLib
41 {
42 namespace Solids
43 {
44 namespace Ehlers
45 {
46 
53 template <int DisplacementDim>
56 
57 template <int DisplacementDim>
59 {
60  static int const KelvinVectorSize =
63  using KelvinVector =
65 
67  : value{stress},
69  I_1{Invariants::trace(stress)},
72  {
73  }
74 
77  double I_1;
78  double J_2;
79  double J_3;
80 
82 };
83 
85 struct OnePlusGamma_pTheta final
86 {
87  OnePlusGamma_pTheta(double const gamma_p, double const theta,
88  double const m_p)
89  : value{1 + gamma_p * theta},
90  pow_m_p{std::pow(value, m_p)},
91  pow_m_p1{pow_m_p / value}
92  {
93  }
94 
95  double const value;
96  double const pow_m_p;
97  double const pow_m_p1;
98 };
99 
100 template <int DisplacementDim>
103  double const sqrtPhi, double const alpha_p, double const beta_p,
104  double const delta_p, double const epsilon_p)
105 {
106  return 3 * (alpha_p * s.I_1 +
107  4 * boost::math::pow<2>(delta_p) * boost::math::pow<3>(s.I_1)) /
108  (2 * sqrtPhi) +
109  3 * beta_p + 6 * epsilon_p * s.I_1;
110 }
111 
112 template <int DisplacementDim>
115  OnePlusGamma_pTheta const& one_gt, double const sqrtPhi,
116  typename SolidEhlers<DisplacementDim>::KelvinVector const& dtheta_dsigma,
117  double const gamma_p, double const m_p)
118 {
119  return (one_gt.pow_m_p *
120  (s.D + s.J_2 * m_p * gamma_p * dtheta_dsigma / one_gt.value)) /
121  (2 * sqrtPhi);
122 }
123 
124 template <int DisplacementDim>
127  double const k)
128 {
129  double const I_1_squared = boost::math::pow<2>(s.I_1);
130  assert(s.J_2 != 0);
131 
132  return std::sqrt(
133  s.J_2 *
134  std::pow(1 + mp.gamma * s.J_3 / (s.J_2 * std::sqrt(s.J_2)),
135  mp.m) +
136  mp.alpha / 2. * I_1_squared +
137  boost::math::pow<2>(mp.delta) *
138  boost::math::pow<2>(I_1_squared)) +
139  mp.beta * s.I_1 + mp.epsilon * I_1_squared - k;
140 }
141 
142 template <int DisplacementDim>
146  double const eps_V,
150  double const eps_p_V,
151  double const eps_p_V_dot,
152  double const eps_p_eff_dot,
153  double const lambda,
154  double const k,
155  MaterialProperties const& mp)
156 {
157  static int const KelvinVectorSize =
160  using KelvinVector =
162 
163  auto const& P_dev = Invariants::deviatoric_projection;
164  auto const& identity2 = Invariants::identity2;
165 
166  double const theta = s.J_3 / (s.J_2 * std::sqrt(s.J_2));
167 
169  // calculate stress residual
170  residual.template segment<KelvinVectorSize>(0).noalias() =
171  s.value / mp.G - 2 * (eps_D - eps_p_D) -
172  mp.K / mp.G * (eps_V - eps_p_V) * identity2;
173 
174  // deviatoric plastic strain
175  KelvinVector const sigma_D_inverse_D =
177  KelvinVector const dtheta_dsigma =
178  theta * sigma_D_inverse_D - 3. / 2. * theta / s.J_2 * s.D;
179 
180  OnePlusGamma_pTheta const one_gt{mp.gamma_p, theta, mp.m_p};
181  double const sqrtPhi = std::sqrt(
182  s.J_2 * one_gt.pow_m_p + mp.alpha_p / 2. * boost::math::pow<2>(s.I_1) +
183  boost::math::pow<2>(mp.delta_p) * boost::math::pow<4>(s.I_1));
185  s, one_gt, sqrtPhi, dtheta_dsigma, mp.gamma_p, mp.m_p);
186  KelvinVector const lambda_flow_D = lambda * flow_D;
187 
188  residual.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
189  eps_p_D_dot - lambda_flow_D;
190 
191  // plastic volume strain
192  {
193  double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
194  s, sqrtPhi, mp.alpha_p, mp.beta_p, mp.delta_p, mp.epsilon_p);
195  residual(2 * KelvinVectorSize, 0) = eps_p_V_dot - lambda * flow_V;
196  }
197 
198  // evolution of plastic equivalent strain
199  residual(2 * KelvinVectorSize + 1) =
200  eps_p_eff_dot -
201  std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
202 
203  // yield function (for plastic multiplier)
204  residual(2 * KelvinVectorSize + 2) = yieldFunction(mp, s, k) / mp.G;
205  return residual;
206 }
207 
208 template <int DisplacementDim>
210  double const dt,
212  double const lambda,
213  MaterialProperties const& mp)
214 {
215  static int const KelvinVectorSize =
218  using KelvinVector =
220  using KelvinMatrix =
222 
223  auto const& P_dev = Invariants::deviatoric_projection;
224  auto const& identity2 = Invariants::identity2;
225 
226  double const theta = s.J_3 / (s.J_2 * std::sqrt(s.J_2));
227  OnePlusGamma_pTheta const one_gt{mp.gamma_p, theta, mp.m_p};
228 
229  // inverse of deviatoric stress tensor
230  if (Invariants::determinant(s.D) == 0)
231  {
232  OGS_FATAL("Determinant is zero. Matrix is non-invertable.");
233  }
234  // inverse of sigma_D
235  KelvinVector const sigma_D_inverse = MathLib::KelvinVector::inverse(s.D);
236  KelvinVector const sigma_D_inverse_D = P_dev * sigma_D_inverse;
237 
238  KelvinVector const dtheta_dsigma =
239  theta * sigma_D_inverse_D - 3. / 2. * theta / s.J_2 * s.D;
240 
241  // deviatoric flow
242  double const sqrtPhi = std::sqrt(
243  s.J_2 * one_gt.pow_m_p + mp.alpha_p / 2. * boost::math::pow<2>(s.I_1) +
244  boost::math::pow<2>(mp.delta_p) * boost::math::pow<4>(s.I_1));
246  s, one_gt, sqrtPhi, dtheta_dsigma, mp.gamma_p, mp.m_p);
247  KelvinVector const lambda_flow_D = lambda * flow_D;
248 
251 
252  // G_11
253  jacobian.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
254  .noalias() = KelvinMatrix::Identity();
255 
256  // G_12
257  jacobian
258  .template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
259  .noalias() = 2 * KelvinMatrix::Identity();
260 
261  // G_13
262  jacobian.template block<KelvinVectorSize, 1>(0, 2 * KelvinVectorSize)
263  .noalias() = mp.K / mp.G * identity2;
264 
265  // G_14 and G_15 are zero
266 
267  // G_21 -- derivative of deviatoric flow
268 
269  double const gm_p = mp.gamma_p * mp.m_p;
270  // intermediate variable for derivative of deviatoric flow
271  KelvinVector const M0 = s.J_2 / one_gt.value * dtheta_dsigma;
272  // derivative of Phi w.r.t. sigma
273  KelvinVector const dPhi_dsigma =
274  one_gt.pow_m_p * (s.D + gm_p * M0) +
275  (mp.alpha_p * s.I_1 +
276  4 * boost::math::pow<2>(mp.delta_p) * boost::math::pow<3>(s.I_1)) *
277  identity2;
278 
279  // intermediate variable for derivative of deviatoric flow
280  KelvinMatrix const M1 =
281  one_gt.pow_m_p *
282  (s.D * dPhi_dsigma.transpose() + gm_p * M0 * dPhi_dsigma.transpose());
283  // intermediate variable for derivative of deviatoric flow
284  KelvinMatrix const M2 =
285  one_gt.pow_m_p * (P_dev + s.D * gm_p * M0.transpose());
286  // second derivative of theta
287  KelvinMatrix const d2theta_dsigma2 =
288  theta * P_dev * sOdotS<DisplacementDim>(sigma_D_inverse) * P_dev +
289  sigma_D_inverse_D * dtheta_dsigma.transpose() -
290  3. / 2. * theta / s.J_2 * P_dev -
291  3. / 2. * dtheta_dsigma / s.J_2 * s.D.transpose() +
292  3. / 2. * theta / boost::math::pow<2>(s.J_2) * s.D * s.D.transpose();
293 
294  // intermediate variable for derivative of deviatoric flow
295  KelvinMatrix const M3 =
296  gm_p * one_gt.pow_m_p1 *
297  ((s.D + (gm_p - mp.gamma_p) * M0) * dtheta_dsigma.transpose() +
298  s.J_2 * d2theta_dsigma2);
299 
300  // derivative of flow_D w.r.t. sigma
301  KelvinMatrix const dflow_D_dsigma =
302  (-M1 / (4 * boost::math::pow<3>(sqrtPhi)) + (M2 + M3) / (2 * sqrtPhi)) *
303  mp.G;
304  jacobian
305  .template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
306  .noalias() = -lambda * dflow_D_dsigma;
307 
308  // G_22
309  jacobian
310  .template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
312  .noalias() = KelvinMatrix::Identity() / dt;
313 
314  // G_23 and G_24 are zero
315 
316  // G_25
317  jacobian
318  .template block<KelvinVectorSize, 1>(KelvinVectorSize,
319  2 * KelvinVectorSize + 2)
320  .noalias() = -flow_D;
321 
322  // G_31
323  {
324  // derivative of flow_V w.r.t. sigma
325  KelvinVector const dflow_V_dsigma =
326  3 * mp.G *
327  (-(mp.alpha_p * s.I_1 +
328  4 * boost::math::pow<2>(mp.delta_p) *
329  boost::math::pow<3>(s.I_1)) /
330  (4 * boost::math::pow<3>(sqrtPhi)) * dPhi_dsigma +
331  (mp.alpha_p * identity2 +
332  12 * boost::math::pow<2>(mp.delta_p * s.I_1) * identity2) /
333  (2 * sqrtPhi) +
334  2 * mp.epsilon_p * identity2);
335 
336  jacobian.template block<1, KelvinVectorSize>(2 * KelvinVectorSize, 0)
337  .noalias() = -lambda * dflow_V_dsigma.transpose();
338  }
339 
340  // G_32 is zero
341 
342  // G_33
343  jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize) = 1. / dt;
344 
345  // G_34 is zero
346 
347  // G_35
348  {
349  double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
350  s, sqrtPhi, mp.alpha_p, mp.beta_p, mp.delta_p, mp.epsilon_p);
351  jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize + 2) = -flow_V;
352  }
353 
354  // increment of effectiv plastic strain
355  double const eff_flow =
356  std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
357 
358  if (eff_flow > 0)
359  {
360  // intermediate variable for derivative of plastic jacobian
361  KelvinVector const eff_flow23_lambda_flow_D =
362  -2 / 3. / eff_flow * lambda_flow_D;
363  // G_41
364  jacobian
365  .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 1, 0)
366  .noalias() = lambda * dflow_D_dsigma * eff_flow23_lambda_flow_D;
367  // G_45
368  jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 2) =
369  eff_flow23_lambda_flow_D.transpose() * flow_D;
370  }
371 
372  // G_42 and G_43 are zero
373 
374  // G_44
375  jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 1) = 1. / dt;
376 
377  // G_51
378  {
379  double const one_gt_pow_m = std::pow(one_gt.value, mp.m);
380  double const gm = mp.gamma * mp.m;
381  // derivative of yield function w.r.t. sigma
382  KelvinVector const dF_dsigma =
383  mp.G * (one_gt_pow_m * (s.D + gm * M0) +
384  (mp.alpha * s.I_1 +
385  4 * boost::math::pow<2>(mp.delta) *
386  boost::math::pow<3>(s.I_1)) *
387  identity2) /
388  (2. * sqrtPhi) +
389  mp.G * (mp.beta + 2 * mp.epsilon_p * s.I_1) * identity2;
390 
391  jacobian
392  .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 2, 0)
393  .noalias() = dF_dsigma.transpose() / mp.G;
394  }
395 
396  // G_54
397  jacobian(2 * KelvinVectorSize + 2, 2 * KelvinVectorSize + 1) =
398  -mp.kappa * mp.hardening_coefficient / mp.G;
399 
400  // G_52, G_53, G_55 are zero
401  return jacobian;
402 }
403 
406 template <int DisplacementDim>
408  double const K, double const G)
409 {
410  static int const KelvinVectorSize =
413 
414  auto const& P_dev = Invariants::deviatoric_projection;
415  auto const& P_sph = Invariants::spherical_projection;
416  auto const& I =
418 
419  return -2. * I * P_dev - 3. * K / G * I * P_sph;
420 }
421 
422 inline double calculateIsotropicHardening(double const kappa,
423  double const hardening_coefficient,
424  double const eps_p_eff)
425 {
426  return kappa * (1. + eps_p_eff * hardening_coefficient);
427 }
428 
429 template <int DisplacementDim>
431  double const G, double const K,
432  typename SolidEhlers<DisplacementDim>::KelvinVector const& sigma_prev,
433  typename SolidEhlers<DisplacementDim>::KelvinVector const& eps,
434  typename SolidEhlers<DisplacementDim>::KelvinVector const& eps_prev,
435  double const eps_V)
436 {
437  static int const KelvinVectorSize =
440  auto const& P_dev = Invariants::deviatoric_projection;
441 
442  // dimensionless initial hydrostatic pressure
443  double const pressure_prev = Invariants::trace(sigma_prev) / (-3. * G);
444  // initial strain invariant
445  double const e_prev = Invariants::trace(eps_prev);
446  // dimensioness hydrostatic stress increment
447  double const pressure = pressure_prev - K / G * (eps_V - e_prev);
448  // dimensionless deviatoric initial stress
449  typename SolidEhlers<DisplacementDim>::KelvinVector const sigma_D_prev =
450  P_dev * sigma_prev / G;
451  // dimensionless deviatoric stress
452  typename SolidEhlers<DisplacementDim>::KelvinVector const sigma_D =
453  sigma_D_prev + 2 * P_dev * (eps - eps_prev);
454  return sigma_D - pressure * Invariants::identity2;
455 }
456 
459 template <typename ResidualVector, typename KelvinVector>
460 std::tuple<KelvinVector, PlasticStrain<KelvinVector>, double>
461 splitSolutionVector(ResidualVector const& solution)
462 {
463  static auto const size = KelvinVector::SizeAtCompileTime;
464  return std::forward_as_tuple(
465  solution.template segment<size>(size * 0),
466  PlasticStrain<KelvinVector>{solution.template segment<size>(size * 1),
467  solution[size * 2], solution[size * 2 + 1]},
468  solution[size * 2 + 2]);
469 }
470 
471 template <int DisplacementDim>
473  NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters,
474  MaterialPropertiesParameters material_properties,
475  std::unique_ptr<DamagePropertiesParameters>&& damage_properties,
476  TangentType tangent_type)
477  : _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
478  _mp(std::move(material_properties)),
479  _damage_properties(std::move(damage_properties)),
480  _tangent_type(tangent_type)
481 {
482 }
483 
484 template <int DisplacementDim>
486  double const /*t*/,
487  ParameterLib::SpatialPosition const& /*x*/,
488  double const /*dt*/,
489  KelvinVector const& eps,
490  KelvinVector const& sigma,
492  material_state_variables) const
493 {
494  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
495  &material_state_variables) != nullptr);
496 
497  auto const& eps_p = static_cast<StateVariables<DisplacementDim> const&>(
498  material_state_variables)
499  .eps_p;
501  auto const& identity2 = Invariants::identity2;
502  return (eps - eps_p.D - eps_p.V / 3 * identity2).dot(sigma) / 2;
503 }
504 
505 template <int DisplacementDim>
506 boost::optional<std::tuple<typename SolidEhlers<DisplacementDim>::KelvinVector,
507  std::unique_ptr<typename MechanicsBase<
508  DisplacementDim>::MaterialStateVariables>,
511  double const t, ParameterLib::SpatialPosition const& x, double const dt,
512  KelvinVector const& eps_prev, KelvinVector const& eps,
513  KelvinVector const& sigma_prev,
515  material_state_variables,
516  double const /*T*/) const
517 {
518  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
519  &material_state_variables) != nullptr);
520 
522  static_cast<StateVariables<DisplacementDim> const&>(
523  material_state_variables);
524  state.setInitialConditions();
525 
527 
528  // volumetric strain
529  double const eps_V = Invariants::trace(eps);
530 
531  auto const& P_dev = Invariants::deviatoric_projection;
532  // deviatoric strain
533  KelvinVector const eps_D = P_dev * eps;
534 
535  // do the evaluation once per function call.
536  MaterialProperties const mp(t, x, _mp);
537 
538  KelvinVector sigma = predict_sigma<DisplacementDim>(mp.G, mp.K, sigma_prev,
539  eps, eps_prev, eps_V);
540 
541  KelvinMatrix tangentStiffness;
542 
544  // Quit early if sigma is zero (nothing to do) or if we are still in elastic
545  // zone.
546  if ((sigma.squaredNorm() == 0 ||
548  mp, s,
550  state.eps_p.eff)) < 0))
551  {
552  tangentStiffness = elasticTangentStiffness<DisplacementDim>(
553  mp.K - 2. / 3 * mp.G, mp.G);
554  }
555  else
556  {
557  // Linear solver for the newton loop is required after the loop with the
558  // same matrix. This saves one decomposition.
559  Eigen::FullPivLU<Eigen::Matrix<double, JacobianResidualSize,
560  JacobianResidualSize, Eigen::RowMajor>>
561  linear_solver;
562 
563  {
564  static int const KelvinVectorSize =
566  DisplacementDim>::value;
567  using KelvinVector =
569  using ResidualVectorType =
570  Eigen::Matrix<double, JacobianResidualSize, 1>;
571  using JacobianMatrix =
572  Eigen::Matrix<double, JacobianResidualSize,
573  JacobianResidualSize, Eigen::RowMajor>;
574 
575  JacobianMatrix jacobian;
576 
577  // Agglomerated solution vector construction. It is later split
578  // into individual parts by splitSolutionVector().
579  ResidualVectorType solution;
580  solution << sigma, state.eps_p.D, state.eps_p.V, state.eps_p.eff, 0;
581 
582  auto const update_residual = [&](ResidualVectorType& residual) {
583 
584  auto const& eps_p_D =
585  solution.template segment<KelvinVectorSize>(
587  KelvinVector const eps_p_D_dot =
588  (eps_p_D - state.eps_p_prev.D) / dt;
589 
590  double const& eps_p_V = solution[KelvinVectorSize * 2];
591  double const eps_p_V_dot =
592  (eps_p_V - state.eps_p_prev.V) / dt;
593 
594  double const& eps_p_eff = solution[KelvinVectorSize * 2 + 1];
595  double const eps_p_eff_dot =
596  (eps_p_eff - state.eps_p_prev.eff) / dt;
597 
598  double const k_hardening = calculateIsotropicHardening(
600  solution[KelvinVectorSize * 2 + 1]);
601  residual = calculatePlasticResidual<DisplacementDim>(
602  eps_D, eps_V, s,
603  solution.template segment<KelvinVectorSize>(
605  eps_p_D_dot, solution[KelvinVectorSize * 2], eps_p_V_dot,
606  eps_p_eff_dot, solution[KelvinVectorSize * 2 + 2],
607  k_hardening, mp);
608  };
609 
610  auto const update_jacobian = [&](JacobianMatrix& jacobian) {
611  jacobian = calculatePlasticJacobian<DisplacementDim>(
612  dt, s, solution[KelvinVectorSize * 2 + 2], mp);
613  };
614 
615  auto const update_solution =
616  [&](ResidualVectorType const& increment) {
617  solution += increment;
619  mp.G * solution.template segment<KelvinVectorSize>(0)};
620  };
621 
622  auto newton_solver = NumLib::NewtonRaphson<
623  decltype(linear_solver), JacobianMatrix,
624  decltype(update_jacobian), ResidualVectorType,
625  decltype(update_residual), decltype(update_solution)>(
626  linear_solver, update_jacobian, update_residual,
627  update_solution, _nonlinear_solver_parameters);
628 
629  auto const success_iterations = newton_solver.solve(jacobian);
630 
631  if (!success_iterations)
632  {
633  return {};
634  }
635 
636  // If the Newton loop didn't run, the linear solver will not be
637  // initialized.
638  // This happens usually for the first iteration of the first
639  // timestep.
640  if (*success_iterations == 0)
641  {
642  linear_solver.compute(jacobian);
643  }
644 
645  std::tie(sigma, state.eps_p, std::ignore) =
646  splitSolutionVector<ResidualVectorType, KelvinVector>(solution);
647  }
648 
649  // Calculate residual derivative w.r.t. strain
650  Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
651  Eigen::RowMajor>
652  dresidual_deps =
653  Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
654  Eigen::RowMajor>::Zero();
655  dresidual_deps.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
656  .noalias() = calculateDResidualDEps<DisplacementDim>(mp.K, mp.G);
657 
659  {
660  tangentStiffness =
661  elasticTangentStiffness<DisplacementDim>(mp.K, mp.G);
662  }
663  else if (_tangent_type == TangentType::Plastic ||
665  {
666  tangentStiffness =
667  mp.G *
668  linear_solver.solve(-dresidual_deps)
669  .template block<KelvinVectorSize, KelvinVectorSize>(0, 0);
671  {
672  tangentStiffness *= 1 - state.damage.value();
673  }
674  }
675  else
676  {
677  OGS_FATAL(
678  "Unimplemented tangent type behaviour for the tangent type "
679  "'%d'.",
680  _tangent_type);
681  }
682  }
683 
684  KelvinVector sigma_final = mp.G * sigma;
685 
686  return {std::make_tuple(
687  sigma_final,
688  std::unique_ptr<
691  static_cast<StateVariables<DisplacementDim> const&>(state)}},
692  tangentStiffness)};
693 }
694 
695 template <int DisplacementDim>
696 std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
698 {
699  return {
700  {"damage.kappa_d", 1,
701  [](typename MechanicsBase<
702  DisplacementDim>::MaterialStateVariables const& state,
703  std::vector<double>& cache) -> std::vector<double> const& {
704  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
705  &state) != nullptr);
706  auto const& ehlers_state =
707  static_cast<StateVariables<DisplacementDim> const&>(state);
708 
709  cache.resize(1);
710  cache.front() = ehlers_state.damage.kappa_d();
711  return cache;
712  }},
713  {"damage.value", 1,
714  [](typename MechanicsBase<
715  DisplacementDim>::MaterialStateVariables const& state,
716  std::vector<double>& cache) -> std::vector<double> const& {
717  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
718  &state) != nullptr);
719  auto const& ehlers_state =
720  static_cast<StateVariables<DisplacementDim> const&>(state);
721 
722  cache.resize(1);
723  cache.front() = ehlers_state.damage.value();
724  return cache;
725  }},
726  {"eps_p.D", KelvinVector::RowsAtCompileTime,
727  [](typename MechanicsBase<
728  DisplacementDim>::MaterialStateVariables const& state,
729  std::vector<double>& cache) -> std::vector<double> const& {
730  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
731  &state) != nullptr);
732  auto const& ehlers_state =
733  static_cast<StateVariables<DisplacementDim> const&>(state);
734 
735  cache.resize(KelvinVector::RowsAtCompileTime);
736  MathLib::toVector<KelvinVector>(cache,
737  KelvinVector::RowsAtCompileTime) =
739  ehlers_state.eps_p.D);
740 
741  return cache;
742  }},
743  {"eps_p.V", 1,
744  [](typename MechanicsBase<
745  DisplacementDim>::MaterialStateVariables const& state,
746  std::vector<double>& cache) -> std::vector<double> const& {
747  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
748  &state) != nullptr);
749  auto const& ehlers_state =
750  static_cast<StateVariables<DisplacementDim> const&>(state);
751 
752  cache.resize(1);
753  cache.front() = ehlers_state.eps_p.V;
754  return cache;
755  }},
756  {"eps_p.eff", 1,
757  [](typename MechanicsBase<
758  DisplacementDim>::MaterialStateVariables const& state,
759  std::vector<double>& cache) -> std::vector<double> const& {
760  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
761  &state) != nullptr);
762  auto const& ehlers_state =
763  static_cast<StateVariables<DisplacementDim> const&>(state);
764 
765  cache.resize(1);
766  cache.front() = ehlers_state.eps_p.eff;
767  return cache;
768  }}};
769 }
770 
771 template class SolidEhlers<2>;
772 template class SolidEhlers<3>;
773 
774 template <>
777 {
779 
780  result(0, 0) = v(0) * v(0);
781  result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
782  result(0, 2) = result(2, 0) = v(5) * v(5) / 2.;
783  result(0, 3) = result(3, 0) = v(0) * v(3);
784  result(0, 4) = result(4, 0) = v(3) * v(5) / std::sqrt(2.);
785  result(0, 5) = result(5, 0) = v(0) * v(5);
786 
787  result(1, 1) = v(1) * v(1);
788  result(1, 2) = result(2, 1) = v(4) * v(4) / 2.;
789  result(1, 3) = result(3, 1) = v(3) * v(1);
790  result(1, 4) = result(4, 1) = v(1) * v(4);
791  result(1, 5) = result(5, 1) = v(3) * v(4) / std::sqrt(2.);
792 
793  result(2, 2) = v(2) * v(2);
794  result(2, 3) = result(3, 2) = v(5) * v(4) / std::sqrt(2.);
795  result(2, 4) = result(4, 2) = v(4) * v(2);
796  result(2, 5) = result(5, 2) = v(5) * v(2);
797 
798  result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
799  result(3, 4) = result(4, 3) =
800  v(3) * v(4) / 2. + v(5) * v(1) / std::sqrt(2.);
801  result(3, 5) = result(5, 3) =
802  v(0) * v(4) / std::sqrt(2.) + v(3) * v(5) / 2.;
803 
804  result(4, 4) = v(1) * v(2) + v(4) * v(4) / 2.;
805  result(4, 5) = result(5, 4) =
806  v(3) * v(2) / std::sqrt(2.) + v(5) * v(4) / 2.;
807 
808  result(5, 5) = v(0) * v(2) + v(5) * v(5) / 2.;
809  return result;
810 }
811 
812 template <>
815 {
817 
818  result(0, 0) = v(0) * v(0);
819  result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
820  result(0, 2) = result(2, 0) = 0;
821  result(0, 3) = result(3, 0) = v(0) * v(3);
822 
823  result(1, 1) = v(1) * v(1);
824  result(1, 2) = result(2, 1) = 0;
825  result(1, 3) = result(3, 1) = v(3) * v(1);
826 
827  result(2, 2) = v(2) * v(2);
828  result(2, 3) = result(3, 2) = 0;
829 
830  result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
831 
832  return result;
833 }
834 
835 } // namespace Ehlers
836 } // namespace Solids
837 } // namespace MaterialLib
double calculateIsotropicHardening(double const kappa, double const hardening_coefficient, double const eps_p_eff)
Definition: Ehlers.cpp:422
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition: KelvinVector.h:79
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
MathLib::KelvinVector::KelvinMatrixType< 3 > sOdotS< 3 >(MathLib::KelvinVector::KelvinVectorType< 3 > const &v)
Definition: Ehlers.cpp:775
Damage damage
damage part of the state.
Definition: Ehlers.h:240
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > calculateDResidualDEps(double const K, double const G)
Definition: Ehlers.cpp:407
NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters
Definition: Ehlers.h:346
Eigen::Matrix< double, JacobianResidualSize, 1 > ResidualVectorType
Definition: Ehlers.h:278
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:59
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
Definition: Ehlers.cpp:64
PlasticStrain< KelvinVector > eps_p
plastic part of the state.
Definition: Ehlers.h:239
Eigen::Matrix< double, JacobianResidualSize, JacobianResidualSize, Eigen::RowMajor > JacobianMatrix
Definition: Ehlers.h:280
double computeFreeEnergyDensity(double const, ParameterLib::SpatialPosition const &, double const, KelvinVector const &eps, KelvinVector const &sigma, typename MechanicsBase< DisplacementDim >::MaterialStateVariables const &material_state_variables) const override
Definition: Ehlers.cpp:485
MathLib::KelvinVector::KelvinMatrixType< 2 > sOdotS< 2 >(MathLib::KelvinVector::KelvinVectorType< 2 > const &v)
Definition: Ehlers.cpp:813
std::vector< typename MechanicsBase< DisplacementDim >::InternalVariable > getInternalVariables() const override
Definition: Ehlers.cpp:697
SolidEhlers(NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters, MaterialPropertiesParameters material_properties, std::unique_ptr< DamagePropertiesParameters > &&damage_properties, TangentType tangent_type)
Definition: Ehlers.cpp:472
SolidEhlers< DisplacementDim >::JacobianMatrix calculatePlasticJacobian(double const dt, PhysicalStressWithInvariants< DisplacementDim > const &s, double const lambda, MaterialProperties const &mp)
Definition: Ehlers.cpp:209
Holds powers of 1 + gamma_p*theta to base 0, m_p, and m_p-1.
Definition: Ehlers.cpp:85
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
Definition: Ehlers.h:294
static Eigen::Matrix< double, KelvinVectorSize, KelvinVectorSize > const deviatoric_projection
Definition: KelvinVector.h:69
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
MathLib::KelvinVector::Invariants< KelvinVectorSize > Invariants
Definition: Ehlers.cpp:62
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.
MaterialPropertiesParameters _mp
Definition: Ehlers.h:348
PhysicalStressWithInvariants(KelvinVector const &stress)
Definition: Ehlers.cpp:66
boost::optional< std::tuple< KelvinVector, std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables >, KelvinMatrix > > integrateStress(double const t, ParameterLib::SpatialPosition const &x, double const dt, KelvinVector const &eps_prev, KelvinVector const &eps, KelvinVector const &sigma_prev, typename MechanicsBase< DisplacementDim >::MaterialStateVariables const &material_state_variables, double const T) const override
Definition: Ehlers.cpp:510
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > sOdotS(MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &v)
SolidEhlers< DisplacementDim >::KelvinVector predict_sigma(double const G, double const K, typename SolidEhlers< DisplacementDim >::KelvinVector const &sigma_prev, typename SolidEhlers< DisplacementDim >::KelvinVector const &eps, typename SolidEhlers< DisplacementDim >::KelvinVector const &eps_prev, double const eps_V)
Definition: Ehlers.cpp:430
double yieldFunction(MaterialProperties const &mp, PhysicalStressWithInvariants< DisplacementDim > const &s, double const k)
Definition: Ehlers.cpp:125
std::tuple< KelvinVector, PlasticStrain< KelvinVector >, double > splitSolutionVector(ResidualVector const &solution)
Definition: Ehlers.cpp:461
OnePlusGamma_pTheta(double const gamma_p, double const theta, double const m_p)
Definition: Ehlers.cpp:87
static double J3(Eigen::Matrix< double, KelvinVectorSize, 1 > const &deviatoric_v)
static double J2(Eigen::Matrix< double, KelvinVectorSize, 1 > const &deviatoric_v)
SolidEhlers< DisplacementDim >::ResidualVectorType calculatePlasticResidual(MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &eps_D, double const eps_V, PhysicalStressWithInvariants< DisplacementDim > const &s, MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &eps_p_D, MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &eps_p_D_dot, double const eps_p_V, double const eps_p_V_dot, double const eps_p_eff_dot, double const lambda, double const k, MaterialProperties const &mp)
Definition: Ehlers.cpp:144
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
static double determinant(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Determinant of a matrix in Kelvin vector representation.
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
Definition: Ehlers.h:292
PlasticStrain< KelvinVector > eps_p_prev
plastic part of the state.
Definition: Ehlers.h:243
static Eigen::Matrix< double, KelvinVectorSize, KelvinVectorSize > const spherical_projection
Definition: KelvinVector.h:77
Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > inverse(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
SolidEhlers< DisplacementDim >::KelvinVector plasticFlowDeviatoricPart(PhysicalStressWithInvariants< DisplacementDim > const &s, OnePlusGamma_pTheta const &one_gt, double const sqrtPhi, typename SolidEhlers< DisplacementDim >::KelvinVector const &dtheta_dsigma, double const gamma_p, double const m_p)
Definition: Ehlers.cpp:113
double plasticFlowVolumetricPart(PhysicalStressWithInvariants< DisplacementDim > const &s, double const sqrtPhi, double const alpha_p, double const beta_p, double const delta_p, double const epsilon_p)
Definition: Ehlers.cpp:101
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49