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 "ComputeLagrangianStressCustomPK2.h"
      11             : 
      12             : registerMooseObject("SolidMechanicsApp", ComputeLagrangianStressCustomPK2);
      13             : 
      14             : InputParameters
      15          52 : ComputeLagrangianStressCustomPK2::validParams()
      16             : {
      17          52 :   InputParameters params = ComputeLagrangianStressPK1::validParams();
      18         104 :   params.addRequiredParam<MaterialPropertyName>("custom_pk2_stress", "The name of the PK2 stress.");
      19         104 :   params.addRequiredParam<MaterialPropertyName>(
      20             :       "custom_pk2_jacobian", "The name of the PK2 Jacobian (w.r.t. the deformation gradient).");
      21          52 :   return params;
      22           0 : }
      23             : 
      24          39 : ComputeLagrangianStressCustomPK2::ComputeLagrangianStressCustomPK2(
      25          39 :     const InputParameters & parameters)
      26             :   : ComputeLagrangianStressPK1(parameters),
      27          39 :     _pk2(getMaterialProperty<RankTwoTensor>("custom_pk2_stress")),
      28         117 :     _dpk2_dF(getMaterialProperty<RankFourTensor>("custom_pk2_jacobian"))
      29             : {
      30          39 : }
      31             : 
      32             : void
      33          39 : ComputeLagrangianStressCustomPK2::initialSetup()
      34             : {
      35             :   // Derives _large_kinematics from the strain calculator's LARGE_KINEMATICS guarantee.
      36          39 :   ComputeLagrangianStressBase::initialSetup();
      37             : 
      38          39 :   if (!_large_kinematics)
      39           0 :     mooseError("This material requires large kinematics; set large_kinematics = true on the "
      40             :                "ComputeLagrangianStrain calculator.");
      41          39 : }
      42             : 
      43             : void
      44       20416 : ComputeLagrangianStressCustomPK2::computeQpPK1Stress()
      45             : {
      46             :   // PK1 = F_ust * PK2(F_stab) -- F-bar enters only through PK2 via the strain calc's
      47             :   // F-bar'd `_F` (which is what NEML2 sees as `forces/F`). The wrap uses F_ust so the
      48             :   // residual matches the new convention.
      49       20416 :   _pk1_stress[_qp] = _F_ust[_qp] * _pk2[_qp];
      50             : 
      51             :   // The R4 algebra for `_pk1_jacobian` / `_pk1_jacobian_bypass_fbar` is read only on
      52             :   // Jacobian sweeps; skip it on residual-only sweeps.
      53       20416 :   if (!_fe_problem.currentlyComputingJacobian() &&
      54       11072 :       !_fe_problem.currentlyComputingResidualAndJacobian())
      55       11072 :     return;
      56             : 
      57             :   // dPK1/d(F_ust) = del_ik PK2_lj  +  F_ust * dPK2/d(F_stab) * d(F_stab)/d(F_ust).
      58             :   usingTensorIndices(i_, j_, k_, l_);
      59             :   const auto I = RankTwoTensor::Identity();
      60        9344 :   const RankFourTensor termA = I.template times<i_, k_, j_, l_>(_pk2[_qp].transpose());
      61             :   // Default path: full F-bar local correction. Equals `_dpk2_dF` when F-bar is off (since
      62             :   // `_d_F_stab_d_F_ust == IdentityFour` then).
      63        9344 :   const RankFourTensor dPK2_dFust = _dpk2_dF[_qp] * _d_F_stab_d_F_ust[_qp];
      64        9344 :   _pk1_jacobian[_qp] = termA + dPK2_dFust.singleProductI(_F_ust[_qp]);
      65             :   // Bypass-F-bar variant for specialty kernels (WPS, homogenization).
      66        9344 :   _pk1_jacobian_bypass_fbar[_qp] = termA + _dpk2_dF[_qp].singleProductI(_F_ust[_qp]);
      67             : }
      68             : 
      69             : void
      70       20416 : ComputeLagrangianStressCustomPK2::computeQpCauchyStress()
      71             : {
      72             :   // sigma = (1/J_ust) F_ust * PK2(F_stab) * F_ust^T. Computed directly because the PK1
      73             :   // base's sigma chain assumes `pk1_jacobian = constitutive dPK1/d(F_stab)`, but here PK1
      74             :   // depends on F_ust directly (Term A above) -- already folded into _pk1_jacobian.
      75       20416 :   const Real J_ust = _F_ust[_qp].det();
      76       20416 :   _cauchy_stress[_qp] = _F_ust[_qp] * _pk2[_qp] * _F_ust[_qp].transpose() / J_ust;
      77             : 
      78       20416 :   if (!_fe_problem.currentlyComputingJacobian() &&
      79       11072 :       !_fe_problem.currentlyComputingResidualAndJacobian())
      80       11072 :     return;
      81             : 
      82             :   // cauchy_jacobian = dsigma/d(dL). sigma depends on dL only through PK2(F_stab(dL))
      83             :   // (F_ust is constant w.r.t. dL since dL is computed from F_stab via the kinematic helper).
      84             :   //   dPK2/d(dL) = dPK2/d(F_stab) * d(F_stab)/d(dL) = _dpk2_dF *
      85             :   //   (_d_deformation_gradient_increment_d_F)^{-1}
      86             :   const RankFourTensor dPK2_d_dL =
      87        9344 :       _dpk2_dF[_qp] * _d_deformation_gradient_increment_d_F[_qp].inverse();
      88       18688 :   _cauchy_jacobian[_qp] = dPK2_d_dL.singleProductI(_F_ust[_qp]).singleProductJ(_F_ust[_qp]) / J_ust;
      89             : }