Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "ComputeLagrangianObjectiveStress.h"
      11             : 
      12             : InputParameters
      13        5268 : ComputeLagrangianObjectiveStress::validParams()
      14             : {
      15        5268 :   InputParameters params = ComputeLagrangianStressCauchy::validParams();
      16             : 
      17        5268 :   params.addClassDescription("Stress update based on the small (engineering) stress");
      18             : 
      19       10536 :   MooseEnum objectiveRate("truesdell jaumann green_naghdi rashid", "truesdell");
      20       10536 :   params.addParam<MooseEnum>(
      21             :       "objective_rate", objectiveRate, "Which type of objective integration to use");
      22       10536 :   params.addParam<bool>(
      23             :       "rotate_old_stress",
      24       10536 :       false,
      25             :       "If true, the rate runs in passthrough mode: `cauchy_stress = small_stress` (no outer "
      26             :       "rotation applied). Use when wrapping a material with `perform_finite_strain_rotations = "
      27             :       "true` whose own internal rotation already produces the correct cumulative Cauchy stress "
      28             :       "(needs the strain calculator to publish a non-identity `rotation_increment` via "
      29             :       "`publish_rotation_increment = true`). `cauchy_jacobian` is still computed via the rate's "
      30             :       "chain rule. Default false preserves the standard objective-rate pipeline.");
      31             : 
      32        5268 :   return params;
      33        5268 : }
      34             : 
      35        3951 : ComputeLagrangianObjectiveStress::ComputeLagrangianObjectiveStress(
      36        3951 :     const InputParameters & parameters)
      37             :   : ComputeLagrangianStressCauchy(parameters),
      38        3951 :     _small_stress(declareProperty<RankTwoTensor>(_base_name + "small_stress")),
      39        7902 :     _small_stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "small_stress")),
      40        3951 :     _small_jacobian(declareProperty<RankFourTensor>(_base_name + "small_jacobian")),
      41        3951 :     _dcauchy_stress_d_eigenstrain(
      42        3951 :         declareProperty<RankFourTensor>(_base_name + "dcauchy_stress_d_eigenstrain")),
      43        7902 :     _cauchy_stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "cauchy_stress")),
      44        7902 :     _mechanical_strain(getMaterialPropertyByName<RankTwoTensor>(_base_name + "mechanical_strain")),
      45        7902 :     _strain_increment(getMaterialPropertyByName<RankTwoTensor>(_base_name + "strain_increment")),
      46        3951 :     _vorticity_increment(
      47        3951 :         getMaterialPropertyByName<RankTwoTensor>(_base_name + "vorticity_increment")),
      48        3951 :     _deformation_gradient_increment(getMaterialPropertyByName<RankTwoTensor>(
      49        3951 :         _base_name + "spatial_deformation_gradient_increment")),
      50        3951 :     _d_deformation_gradient_increment_d_F(getMaterialPropertyByName<RankFourTensor>(
      51        3951 :         _base_name + "d_spatial_deformation_gradient_increment_d_deformation_gradient")),
      52        3951 :     _d_vorticity_increment_d_F(getMaterialPropertyByName<RankFourTensor>(
      53        3951 :         _base_name + "d_vorticity_increment_d_deformation_gradient")),
      54        3951 :     _def_grad(getMaterialPropertyByName<RankTwoTensor>(_base_name + "deformation_gradient")),
      55        7902 :     _def_grad_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "deformation_gradient")),
      56        7902 :     _objective_rate(getParam<MooseEnum>("objective_rate")),
      57       15804 :     _rotate_old_stress(getParam<bool>("rotate_old_stress"))
      58             : {
      59             :   // Only the Green-Naghdi rate consumes the polar-decomposition rotation. Fetch these
      60             :   // (which marks `rotation` active and triggers the strain calc's polar decomposition) only
      61             :   // for that rate -- Truesdell/Jaumann/Rashid leave them null and skip the eigensolve.
      62        3951 :   if (_objective_rate == "green_naghdi")
      63             :   {
      64         192 :     _rotation = &getMaterialPropertyByName<RankTwoTensor>(_base_name + "rotation");
      65         192 :     _rotation_old = &getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "rotation");
      66          96 :     _d_rotation_d_F = &getMaterialPropertyByName<RankFourTensor>(
      67         192 :         _base_name + "d_rotation_d_deformation_gradient");
      68             :   }
      69        3951 : }
      70             : 
      71             : void
      72      655730 : ComputeLagrangianObjectiveStress::initQpStatefulProperties()
      73             : {
      74      655730 :   ComputeLagrangianStressBase::initQpStatefulProperties();
      75             : 
      76      655730 :   _small_stress[_qp].zero();
      77      655730 :   _cauchy_stress[_qp].zero();
      78      655730 : }
      79             : 
      80             : void
      81    49631698 : ComputeLagrangianObjectiveStress::computeQpCauchyStress()
      82             : {
      83    49631698 :   computeQpSmallStress();
      84             : 
      85    49631698 :   const bool need_jacobian = _fe_problem.currentlyComputingJacobian() ||
      86    47999651 :                              _fe_problem.currentlyComputingResidualAndJacobian();
      87             : 
      88    49631698 :   if (!_large_kinematics)
      89             :   {
      90    16049094 :     _cauchy_stress[_qp] = _small_stress[_qp];
      91             :     // Small kinematics: no objective advection. dsigma/d_eigenstrain = -small_jacobian
      92             :     // (the negative because mechanical_strain = total_strain - eigenstrain).
      93             :     // Both Jacobian-side assignments are skipped on residual-only sweeps; the kernel
      94             :     // consumes neither during residual assembly.
      95    16049094 :     if (need_jacobian)
      96             :     {
      97      505199 :       _cauchy_jacobian[_qp] = _small_jacobian[_qp];
      98      505199 :       _dcauchy_stress_d_eigenstrain[_qp] = -_small_jacobian[_qp];
      99             :     }
     100             :   }
     101             :   else
     102             :   {
     103    33582604 :     const RankTwoTensor dS = _small_stress[_qp] - _small_stress_old[_qp];
     104             : 
     105             :     // Gather the per-qp quantities the rate needs, dispatch to the selected objective rate, and
     106             :     // scatter the outputs back into our properties.
     107    33582604 :     LagrangianObjectiveRates::Inputs in;
     108    33582604 :     in.dS = dS;
     109    33582604 :     in.cauchy_stress_old = _cauchy_stress_old[_qp];
     110    33582604 :     in.small_jacobian = _small_jacobian[_qp];
     111    33582604 :     in.dL = _deformation_gradient_increment[_qp];
     112    33582604 :     in.dW = _vorticity_increment[_qp];
     113    33582604 :     in.d_dL_d_F = _d_deformation_gradient_increment_d_F[_qp];
     114    33582604 :     in.d_dW_d_F = _d_vorticity_increment_d_F[_qp];
     115             :     // Green-Naghdi only: the rotation properties (and inv_df / F^{-1}) are fetched only for that
     116             :     // rate, so `_rotation` is non-null exactly when `greenNaghdi` will read them.
     117    33582604 :     if (_rotation)
     118             :     {
     119     3159348 :       in.rotation = (*_rotation)[_qp];
     120     3159348 :       in.rotation_old = (*_rotation_old)[_qp];
     121     3159348 :       in.d_rotation_d_F = (*_d_rotation_d_F)[_qp];
     122     3159348 :       in.inv_df = _inv_df[_qp];
     123     3159348 :       in.inv_def_grad = _inv_def_grad[_qp];
     124             :     }
     125             : 
     126    33582604 :     const auto out = LagrangianObjectiveRates::compute(_objective_rate, in, need_jacobian);
     127    33582604 :     _cauchy_stress[_qp] = out.cauchy_stress;
     128    33582604 :     if (need_jacobian)
     129             :     {
     130     1126848 :       _cauchy_jacobian[_qp] = out.cauchy_jacobian;
     131     1126848 :       _dcauchy_stress_d_eigenstrain[_qp] = out.dcauchy_stress_d_eigenstrain;
     132             :     }
     133             : 
     134    33582604 :     if (_rotate_old_stress)
     135             :     {
     136             :       // Passthrough mode: the wrapped stress material (e.g. `ComputeMultiPlasticityStress`
     137             :       // with `perform_finite_strain_rotations = true`, fed our published
     138             :       // `_rotation_increment`) has already produced the correctly-rotated cumulative Cauchy
     139             :       // stress. Use `_small_stress` directly and discard the rate's own outer rotation of
     140             :       // the stress. `_cauchy_jacobian` from the rate's chain rule is still correct because
     141             :       // it computes `d sigma_cauchy / d(dL)` from `_small_jacobian = dDeltasigma_const/deps` plus
     142             :       // the same rotation-derivative terms -- algebraically the same chain whether we view sigma as
     143             :       // produced by the rate's `R * (sigma_old + Deltasigma_const) * R^T` formula or by the wrapped
     144             :       // material's `R * (s_old + Deltasigma_const) * R^T = small_stress` storage.
     145      103168 :       _cauchy_stress[_qp] = _small_stress[_qp];
     146             :     }
     147             :   }
     148    49631698 : }