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 "ComputeLagrangianStressCauchy.h"
      11             : 
      12             : InputParameters
      13        5320 : ComputeLagrangianStressCauchy::validParams()
      14             : {
      15        5320 :   InputParameters params = ComputeLagrangianStressBase::validParams();
      16             : 
      17        5320 :   params.addClassDescription("Stress update based on the Cauchy stress");
      18             : 
      19        5320 :   return params;
      20           0 : }
      21             : 
      22        3990 : ComputeLagrangianStressCauchy::ComputeLagrangianStressCauchy(const InputParameters & parameters)
      23             :   : ComputeLagrangianStressBase(parameters),
      24        3990 :     _inv_df(getMaterialPropertyByName<RankTwoTensor>(_base_name +
      25             :                                                      "inverse_incremental_deformation_gradient")),
      26        3990 :     _inv_def_grad(
      27        3990 :         getMaterialPropertyByName<RankTwoTensor>(_base_name + "inverse_deformation_gradient")),
      28        7980 :     _F(getMaterialPropertyByName<RankTwoTensor>(_base_name + "deformation_gradient")),
      29        3990 :     _F_ust_inv(getMaterialPropertyByName<RankTwoTensor>(
      30        3990 :         _base_name + "inverse_unstabilized_deformation_gradient")),
      31        3990 :     _F_ust_det(
      32        7980 :         getMaterialPropertyByName<Real>(_base_name + "det_unstabilized_deformation_gradient"))
      33             : {
      34        3990 : }
      35             : 
      36             : void
      37    49655570 : ComputeLagrangianStressCauchy::computeQpStressUpdate()
      38             : {
      39    49655570 :   computeQpCauchyStress();
      40    49655570 :   computeQpPK1Stress(); // This could be "switched"
      41    49655570 : }
      42             : 
      43             : void
      44    49655570 : ComputeLagrangianStressCauchy::computeQpPK1Stress()
      45             : {
      46             :   // Wrap sigma -> PK1 using the *unstabilized* F (= F_actual at alpha = 1). F-bar enters
      47             :   // only through sigma via the strain calc's F-bar'd `_f_inv`, NOT through this kinematic
      48             :   // wrap -- this matches OLD's `StressDivergenceTensors`-on-displaced-mesh residual,
      49             :   // which integrates `gradTest_spatial * sigma dV` on the actual deformed mesh (so the
      50             :   // implicit "wrap F" is F_actual, NOT the F-bar'd F).
      51             :   //
      52             :   // The R4 algebra for `_pk1_jacobian` and `_pk1_jacobian_bypass_fbar` is Jacobian-only;
      53             :   // the wrap to `_pk1_stress` is not. Gate accordingly.
      54    49655570 :   const bool need_jacobian = _fe_problem.currentlyComputingJacobian() ||
      55    48010723 :                              _fe_problem.currentlyComputingResidualAndJacobian();
      56    49655570 :   if (_large_kinematics)
      57             :   {
      58             :     // F_ust^{-1} and det(F_ust) are published once per qp by the strain calculator; reuse them
      59             :     // instead of recomputing a 3x3 inverse + determinant here on every residual/Jacobian sweep.
      60    33606476 :     const RankTwoTensor & F_ust_inv = _F_ust_inv[_qp];
      61    33606476 :     const RankTwoTensor F_ust_invT = F_ust_inv.transpose();
      62    33606476 :     const Real det_F_ust = _F_ust_det[_qp];
      63    33606476 :     _pk1_stress[_qp] = det_F_ust * _cauchy_stress[_qp] * F_ust_invT;
      64             : 
      65    33606476 :     if (!need_jacobian)
      66             :       return;
      67             : 
      68             :     usingTensorIndices(i_, j_, k_, l_);
      69             :     // Geometric pieces of dPK1/d(F_ust) (independent of F-bar):
      70     1139648 :     const RankFourTensor geom = _pk1_stress[_qp].outerProduct(F_ust_invT) -
      71     2279296 :                                 _pk1_stress[_qp].times<i_, l_, j_, k_>(F_ust_inv);
      72             :     // sigma chain WITHOUT the F-bar `_d_F_stab_d_F_ust` factor -- for specialty kernels
      73             :     // whose coupled variable adds to F_ust AFTER F-bar (WPS strain_zz, homogenization
      74             :     // macro_grad), so the perturbation bypasses F-bar's chain.
      75             :     const RankFourTensor sigma_chain_no_fbar =
      76     1139648 :         det_F_ust * (_cauchy_jacobian[_qp] * _d_deformation_gradient_increment_d_F[_qp])
      77     1139648 :                         .singleProductJ(F_ust_inv);
      78     1139648 :     _pk1_jacobian_bypass_fbar[_qp] = geom + sigma_chain_no_fbar;
      79             :     // sigma chain WITH the F-bar chain -- the default `_pk1_jacobian` for disp Jacobian.
      80     1139648 :     _pk1_jacobian[_qp] =
      81     1139648 :         geom + det_F_ust * (_cauchy_jacobian[_qp] * _d_deformation_gradient_increment_d_F[_qp] *
      82     1139648 :                             _d_F_stab_d_F_ust[_qp])
      83     2279296 :                                .singleProductJ(F_ust_inv);
      84             :   }
      85             :   else
      86             :   {
      87             :     // Small kinematics: PK1 = sigma (no wrap). The F-bar chain enters via
      88             :     // `_d_F_stab_d_F_ust` for the disp Jacobian path; specialty paths bypass it.
      89    16049094 :     _pk1_stress[_qp] = _cauchy_stress[_qp];
      90    16049094 :     if (!need_jacobian)
      91             :       return;
      92      505199 :     _pk1_jacobian[_qp] = _cauchy_jacobian[_qp] * _d_F_stab_d_F_ust[_qp];
      93      505199 :     _pk1_jacobian_bypass_fbar[_qp] = _cauchy_jacobian[_qp];
      94             :   }
      95             : }