OGS 6.1.0-1721-g6382411ad
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  ProcessLib::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, ProcessLib::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  return {};
633 
634  // If the Newton loop didn't run, the linear solver will not be
635  // initialized.
636  // This happens usually for the first iteration of the first
637  // timestep.
638  if (*success_iterations == 0)
639  linear_solver.compute(jacobian);
640 
641  std::tie(sigma, state.eps_p, std::ignore) =
642  splitSolutionVector<ResidualVectorType, KelvinVector>(solution);
643  }
644 
645  // Calculate residual derivative w.r.t. strain
646  Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
647  Eigen::RowMajor>
648  dresidual_deps =
649  Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
650  Eigen::RowMajor>::Zero();
651  dresidual_deps.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
652  .noalias() = calculateDResidualDEps<DisplacementDim>(mp.K, mp.G);
653 
655  {
656  tangentStiffness =
657  elasticTangentStiffness<DisplacementDim>(mp.K, mp.G);
658  }
659  else if (_tangent_type == TangentType::Plastic ||
661  {
662  tangentStiffness =
663  mp.G *
664  linear_solver.solve(-dresidual_deps)
665  .template block<KelvinVectorSize, KelvinVectorSize>(0, 0);
667  {
668  tangentStiffness *= 1 - state.damage.value();
669  }
670  }
671  else
672  {
673  OGS_FATAL(
674  "Unimplemented tangent type behaviour for the tangent type "
675  "'%d'.",
676  _tangent_type);
677  }
678  }
679 
680  KelvinVector sigma_final = mp.G * sigma;
681 
682  return {std::make_tuple(
683  sigma_final,
684  std::unique_ptr<
687  static_cast<StateVariables<DisplacementDim> const&>(state)}},
688  tangentStiffness)};
689 }
690 
691 template <int DisplacementDim>
692 std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
694 {
695  return {
696  {"damage.kappa_d", 1,
697  [](typename MechanicsBase<
698  DisplacementDim>::MaterialStateVariables const& state,
699  std::vector<double>& cache) -> std::vector<double> const& {
700  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
701  &state) != nullptr);
702  auto const& ehlers_state =
703  static_cast<StateVariables<DisplacementDim> const&>(state);
704 
705  cache.resize(1);
706  cache.front() = ehlers_state.damage.kappa_d();
707  return cache;
708  }},
709  {"damage.value", 1,
710  [](typename MechanicsBase<
711  DisplacementDim>::MaterialStateVariables const& state,
712  std::vector<double>& cache) -> std::vector<double> const& {
713  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
714  &state) != nullptr);
715  auto const& ehlers_state =
716  static_cast<StateVariables<DisplacementDim> const&>(state);
717 
718  cache.resize(1);
719  cache.front() = ehlers_state.damage.value();
720  return cache;
721  }},
722  {"eps_p.D", KelvinVector::RowsAtCompileTime,
723  [](typename MechanicsBase<
724  DisplacementDim>::MaterialStateVariables const& state,
725  std::vector<double>& cache) -> std::vector<double> const& {
726  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
727  &state) != nullptr);
728  auto const& ehlers_state =
729  static_cast<StateVariables<DisplacementDim> const&>(state);
730 
731  cache.resize(KelvinVector::RowsAtCompileTime);
732  MathLib::toVector<KelvinVector>(cache,
733  KelvinVector::RowsAtCompileTime) =
735  ehlers_state.eps_p.D);
736 
737  return cache;
738  }},
739  {"eps_p.V", 1,
740  [](typename MechanicsBase<
741  DisplacementDim>::MaterialStateVariables const& state,
742  std::vector<double>& cache) -> std::vector<double> const& {
743  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
744  &state) != nullptr);
745  auto const& ehlers_state =
746  static_cast<StateVariables<DisplacementDim> const&>(state);
747 
748  cache.resize(1);
749  cache.front() = ehlers_state.eps_p.V;
750  return cache;
751  }},
752  {"eps_p.eff", 1,
753  [](typename MechanicsBase<
754  DisplacementDim>::MaterialStateVariables const& state,
755  std::vector<double>& cache) -> std::vector<double> const& {
756  assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
757  &state) != nullptr);
758  auto const& ehlers_state =
759  static_cast<StateVariables<DisplacementDim> const&>(state);
760 
761  cache.resize(1);
762  cache.front() = ehlers_state.eps_p.eff;
763  return cache;
764  }}};
765 }
766 
767 template class SolidEhlers<2>;
768 template class SolidEhlers<3>;
769 
770 template <>
773 {
775 
776  result(0, 0) = v(0) * v(0);
777  result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
778  result(0, 2) = result(2, 0) = v(5) * v(5) / 2.;
779  result(0, 3) = result(3, 0) = v(0) * v(3);
780  result(0, 4) = result(4, 0) = v(3) * v(5) / std::sqrt(2.);
781  result(0, 5) = result(5, 0) = v(0) * v(5);
782 
783  result(1, 1) = v(1) * v(1);
784  result(1, 2) = result(2, 1) = v(4) * v(4) / 2.;
785  result(1, 3) = result(3, 1) = v(3) * v(1);
786  result(1, 4) = result(4, 1) = v(1) * v(4);
787  result(1, 5) = result(5, 1) = v(3) * v(4) / std::sqrt(2.);
788 
789  result(2, 2) = v(2) * v(2);
790  result(2, 3) = result(3, 2) = v(5) * v(4) / std::sqrt(2.);
791  result(2, 4) = result(4, 2) = v(4) * v(2);
792  result(2, 5) = result(5, 2) = v(5) * v(2);
793 
794  result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
795  result(3, 4) = result(4, 3) =
796  v(3) * v(4) / 2. + v(5) * v(1) / std::sqrt(2.);
797  result(3, 5) = result(5, 3) =
798  v(0) * v(4) / std::sqrt(2.) + v(3) * v(5) / 2.;
799 
800  result(4, 4) = v(1) * v(2) + v(4) * v(4) / 2.;
801  result(4, 5) = result(5, 4) =
802  v(3) * v(2) / std::sqrt(2.) + v(5) * v(4) / 2.;
803 
804  result(5, 5) = v(0) * v(2) + v(5) * v(5) / 2.;
805  return result;
806 }
807 
808 template <>
811 {
813 
814  result(0, 0) = v(0) * v(0);
815  result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
816  result(0, 2) = result(2, 0) = 0;
817  result(0, 3) = result(3, 0) = v(0) * v(3);
818 
819  result(1, 1) = v(1) * v(1);
820  result(1, 2) = result(2, 1) = 0;
821  result(1, 3) = result(3, 1) = v(3) * v(1);
822 
823  result(2, 2) = v(2) * v(2);
824  result(2, 3) = result(3, 2) = 0;
825 
826  result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
827 
828  return result;
829 }
830 
831 } // namespace Ehlers
832 } // namespace Solids
833 } // 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:771
Damage damage
damage part of the state.
Definition: Ehlers.h:243
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > calculateDResidualDEps(double const K, double const G)
Definition: Ehlers.cpp:407
NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters
Definition: Ehlers.h:343
Eigen::Matrix< double, JacobianResidualSize, 1 > ResidualVectorType
Definition: Ehlers.h:281
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:242
Eigen::Matrix< double, JacobianResidualSize, JacobianResidualSize, Eigen::RowMajor > JacobianMatrix
Definition: Ehlers.h:283
MathLib::KelvinVector::KelvinMatrixType< 2 > sOdotS< 2 >(MathLib::KelvinVector::KelvinVectorType< 2 > const &v)
Definition: Ehlers.cpp:809
std::vector< typename MechanicsBase< DisplacementDim >::InternalVariable > getInternalVariables() const override
Definition: Ehlers.cpp:693
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:297
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:345
PhysicalStressWithInvariants(KelvinVector const &stress)
Definition: Ehlers.cpp:66
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
double computeFreeEnergyDensity(double const, ProcessLib::SpatialPosition const &, double const, KelvinVector const &eps, KelvinVector const &sigma, typename MechanicsBase< DisplacementDim >::MaterialStateVariables const &material_state_variables) const override
Definition: Ehlers.cpp:485
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:71
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:295
PlasticStrain< KelvinVector > eps_p_prev
plastic part of the state.
Definition: Ehlers.h:246
boost::optional< std::tuple< KelvinVector, std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables >, KelvinMatrix > > integrateStress(double const t, ProcessLib::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
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