524 auto const& eps_m = std::get<MPL::SymmetricTensor<DisplacementDim>>(
525 variable_array[
static_cast<int>(MPL::Variable::mechanical_strain)]);
526 auto const& eps_m_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
527 variable_array_prev[
static_cast<int>(
528 MPL::Variable::mechanical_strain)]);
529 auto const& sigma_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
530 variable_array_prev[
static_cast<int>(MPL::Variable::stress)]);
532 assert(
dynamic_cast<StateVariables<DisplacementDim> const*
>(
533 &material_state_variables) !=
nullptr);
535 StateVariables<DisplacementDim> state =
536 static_cast<StateVariables<DisplacementDim> const&
>(
537 material_state_variables);
538 state.setInitialConditions();
543 double const eps_V = Invariants::trace(eps_m);
545 auto const& P_dev = Invariants::deviatoric_projection;
550 MaterialProperties
const mp(t, x,
_mp);
553 mp.G, mp.K, sigma_prev, eps_m, eps_m_prev, eps_V);
557 PhysicalStressWithInvariants<DisplacementDim> s{mp.G * sigma};
560 if ((sigma.squaredNorm() == 0 ||
564 state.eps_p.eff)) < 0))
566 tangentStiffness = elasticTangentStiffness<DisplacementDim>(
567 mp.K - 2. / 3 * mp.G, mp.G);
584 Eigen::Matrix<double, JacobianResidualSize, 1>;
594 solution << sigma, state.eps_p.D, state.eps_p.V, state.eps_p.eff, 0;
598 auto const& eps_p_D =
599 solution.template segment<KelvinVectorSize>(
602 (eps_p_D - state.eps_p_prev.D) / dt;
605 double const eps_p_V_dot = (eps_p_V - state.eps_p_prev.V) / dt;
608 double const eps_p_eff_dot =
609 (eps_p_eff - state.eps_p_prev.eff) / dt;
612 mp.kappa, mp.hardening_coefficient,
614 residual = calculatePlasticResidual<DisplacementDim>(
616 solution.template segment<KelvinVectorSize>(
625 jacobian = calculatePlasticJacobian<DisplacementDim>(
629 auto const update_solution =
632 solution += increment;
633 s = PhysicalStressWithInvariants<DisplacementDim>{
634 mp.G * solution.template segment<KelvinVectorSize>(0)};
640 decltype(update_residual), decltype(update_solution)>(
641 linear_solver, update_jacobian, update_residual,
644 auto const success_iterations = newton_solver.solve(jacobian);
646 if (!success_iterations)
655 if (*success_iterations == 0)
657 linear_solver.compute(jacobian);
660 std::tie(sigma, state.eps_p, std::ignore) =
661 splitSolutionVector<ResidualVectorType, KelvinVector>(solution);
669 Eigen::RowMajor>::Zero();
670 dresidual_deps.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
671 .noalias() = calculateDResidualDEps<DisplacementDim>(mp.K, mp.G);
676 elasticTangentStiffness<DisplacementDim>(mp.K, mp.G);
683 linear_solver.
solve(-dresidual_deps)
684 .template block<KelvinVectorSize, KelvinVectorSize>(0, 0);
687 tangentStiffness *= 1 - state.damage.value();
693 "Unimplemented tangent type behaviour for the tangent type "
701 return {std::make_tuple(
704 typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
705 new StateVariables<DisplacementDim>{
706 static_cast<StateVariables<DisplacementDim> const&
>(state)}},
static int const KelvinVectorSize
static int const JacobianResidualSize
Eigen::Matrix< double, JacobianResidualSize, JacobianResidualSize, Eigen::RowMajor > JacobianMatrix
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
Eigen::Matrix< double, JacobianResidualSize, 1 > ResidualVectorType
std::optional< int > solve(JacobianMatrix &jacobian) const
double calculateIsotropicHardening(double const kappa, double const hardening_coefficient, double const eps_p_eff)
double yieldFunction(MaterialProperties const &mp, PhysicalStressWithInvariants< DisplacementDim > const &s, double const k)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.