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 "LagrangianStressDivergenceBase.h"
      11             : #include "RankFourTensor.h"
      12             : 
      13             : InputParameters
      14        6774 : LagrangianStressDivergenceBase::validParams()
      15             : {
      16        6774 :   InputParameters params = KernelScalarBase::validParams();
      17             : 
      18       13548 :   params.addRequiredParam<unsigned int>("component", "Which direction this kernel acts in");
      19       13548 :   params.addRequiredCoupledVar("displacements", "The displacement components");
      20             : 
      21       20322 :   params.addDeprecatedParam<bool>(
      22             :       "large_kinematics",
      23       13548 :       false,
      24             :       "Use large displacement kinematics",
      25             :       "large_kinematics is no longer set on the stress-divergence kernel; it is derived from the "
      26             :       "ComputeLagrangianStrain calculator (the single source of truth) via the LARGE_KINEMATICS "
      27             :       "guarantee. Remove it here and set it only on the strain calculator.");
      28       13548 :   params.addParam<bool>("stabilize_strain", false, "Average the volumetric strains");
      29       13548 :   MooseEnum F_bar_mode("total incremental", "total");
      30       13548 :   params.addParam<MooseEnum>(
      31             :       "F_bar_mode",
      32             :       F_bar_mode,
      33             :       "Which F gets the F-bar volumetric correction (must match the strain calc's setting). "
      34             :       "'total' (default) reproduces existing behavior; 'incremental' makes the F-bar element-"
      35             :       "average operate on the incremental F so cumulative strain matches OLD's "
      36             :       "`ComputeFiniteStrain` + `volumetric_locking_correction = true`.");
      37             : 
      38       13548 :   params.addParam<std::string>("base_name", "Material property base name");
      39             : 
      40       13548 :   params.addCoupledVar("temperature",
      41             :                        "The name of the temperature variable used in the "
      42             :                        "ComputeThermalExpansionEigenstrain.  (Not required for "
      43             :                        "simulations without temperature coupling.)");
      44             : 
      45       13548 :   params.addParam<std::vector<MaterialPropertyName>>(
      46             :       "eigenstrain_names",
      47             :       {},
      48             :       "List of eigenstrains used in the strain calculation. Used for computing their derivatives "
      49             :       "for off-diagonal Jacobian terms.");
      50             : 
      51       13548 :   params.addCoupledVar("out_of_plane_strain",
      52             :                        "The out-of-plane strain variable for weak plane stress formulation.");
      53             : 
      54        6774 :   return params;
      55        6774 : }
      56             : 
      57        3387 : LagrangianStressDivergenceBase::LagrangianStressDivergenceBase(const InputParameters & parameters)
      58             :   : JvarMapKernelInterface<DerivativeMaterialInterface<KernelScalarBase>>(parameters),
      59             :     GuaranteeConsumer(this),
      60             :     // Derived from the strain calculator's guarantee in initialSetup(); the local parameter is
      61             :     // deprecated and only consulted for a consistency cross-check.
      62        3387 :     _large_kinematics(false),
      63        3387 :     _stabilize_strain(getParam<bool>("stabilize_strain")),
      64       10077 :     _F_bar_mode(getParam<MooseEnum>("F_bar_mode") == "incremental" ? FBarMode::Incremental
      65             :                                                                    : FBarMode::Total),
      66        6822 :     _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
      67        6774 :     _alpha(getParam<unsigned int>("component")),
      68        3387 :     _ndisp(coupledComponents("displacements")),
      69        3387 :     _disp_nums(_ndisp),
      70        3387 :     _avg_grad_trial(_ndisp),
      71        3387 :     _avg_grad_spatial_phi(_ndisp),
      72        3387 :     _avg_test_phi_cross(3, std::vector<std::vector<std::vector<Real>>>(3)),
      73        3387 :     _F_ust(
      74        3387 :         getMaterialPropertyByName<RankTwoTensor>(_base_name + "unstabilized_deformation_gradient")),
      75        6774 :     _F_ust_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name +
      76             :                                                            "unstabilized_deformation_gradient")),
      77        6774 :     _F_avg(getMaterialPropertyByName<RankTwoTensor>(_base_name + "average_deformation_gradient")),
      78        6774 :     _f_inv(getMaterialPropertyByName<RankTwoTensor>(_base_name +
      79             :                                                     "inverse_incremental_deformation_gradient")),
      80        6774 :     _F_inv(getMaterialPropertyByName<RankTwoTensor>(_base_name + "inverse_deformation_gradient")),
      81        3387 :     _F(getMaterialPropertyByName<RankTwoTensor>(_base_name + "deformation_gradient")),
      82        6774 :     _F_actual(getMaterialPropertyByName<RankTwoTensor>(_base_name + "actual_deformation_gradient")),
      83        3387 :     _d_deformation_gradient_increment_d_F(getMaterialPropertyByName<RankFourTensor>(
      84        3387 :         _base_name + "d_spatial_deformation_gradient_increment_d_deformation_gradient")),
      85        3387 :     _d_F_d_grad_u(getMaterialPropertyByName<RankFourTensor>(
      86        3387 :         _base_name + "d_deformation_gradient_d_grad_displacement")),
      87        3387 :     _d_F_stab_d_F_ust(
      88        3387 :         getMaterialPropertyByName<RankFourTensor>(_base_name + "d_F_stab_d_F_unstabilized")),
      89        3387 :     _d_F_stab_d_F_avg(
      90        3387 :         getMaterialPropertyByName<RankFourTensor>(_base_name + "d_F_stab_d_F_average")),
      91        3746 :     _temperature(isCoupled("temperature") ? getVar("temperature", 0) : nullptr),
      92        3419 :     _out_of_plane_strain(isCoupled("out_of_plane_strain") ? getVar("out_of_plane_strain", 0)
      93        6774 :                                                           : nullptr)
      94             : {
      95             :   // Do the vector coupling of the displacements
      96       12326 :   for (unsigned int i = 0; i < _ndisp; i++)
      97        8939 :     _disp_nums[i] = coupled("displacements", i);
      98             : 
      99             :   // We need to use identical discretizations for all displacement components
     100        3387 :   auto order_x = getVar("displacements", 0)->order();
     101        8939 :   for (unsigned int i = 1; i < _ndisp; i++)
     102             :   {
     103       11104 :     if (getVar("displacements", i)->order() != order_x)
     104           0 :       mooseError("The Lagrangian StressDivergence kernels require equal "
     105             :                  "order interpolation for all displacements.");
     106             :   }
     107             : 
     108             :   // fetch eigenstrain derivatives
     109             :   const auto nvar = _coupled_moose_vars.size();
     110        3387 :   _deigenstrain_dargs.resize(nvar);
     111       12717 :   for (std::size_t i = 0; i < nvar; ++i)
     112       20016 :     for (auto eigenstrain_name : getParam<std::vector<MaterialPropertyName>>("eigenstrain_names"))
     113        2712 :       _deigenstrain_dargs[i].push_back(&getMaterialPropertyDerivative<RankTwoTensor>(
     114        1356 :           eigenstrain_name, _coupled_moose_vars[i]->name()));
     115             : 
     116             :   // The direct-chain temperature off-diagonal Jacobian needs d_sigma/d_eigenstrain. Only
     117             :   // fetch it when eigenstrains are coupled -- otherwise the temperature Jacobian short-
     118             :   // circuits to zero and this property would be a needless dependency that other Cauchy-
     119             :   // providing materials would have to publish.
     120        6774 :   if (!getParam<std::vector<MaterialPropertyName>>("eigenstrain_names").empty())
     121         407 :     _dcauchy_stress_d_eigenstrain =
     122         814 :         &getMaterialPropertyByName<RankFourTensor>(_base_name + "dcauchy_stress_d_eigenstrain");
     123             : 
     124             :   // The F-bar spatial push-forward consumes F_ust^{-1}/det(F_ust) published by the strain calc.
     125             :   // Fetch them (and thereby mark them active for the strain material's isPropertyActive gate)
     126             :   // only when stabilization is on -- the only mode where the push-forward runs.
     127        3387 :   if (_stabilize_strain)
     128             :   {
     129        1204 :     _F_ust_inv = &getMaterialPropertyByName<RankTwoTensor>(
     130        1204 :         _base_name + "inverse_unstabilized_deformation_gradient");
     131        1204 :     _F_ust_det =
     132        1204 :         &getMaterialPropertyByName<Real>(_base_name + "det_unstabilized_deformation_gradient");
     133        3612 :     _d_nl_fbar = &getMaterialPropertyByName<RankFourTensor>(_base_name + "d_nl_fbar_operator");
     134             :   }
     135        3387 : }
     136             : 
     137             : void
     138        3381 : LagrangianStressDivergenceBase::initialSetup()
     139             : {
     140        3381 :   JvarMapKernelInterface<DerivativeMaterialInterface<KernelScalarBase>>::initialSetup();
     141             : 
     142             :   // Derive the kinematics regime from the strain calculator's LARGE_KINEMATICS guarantee -- the
     143             :   // single source of truth. hasGuaranteedMaterialProperty is block-restricted (per subdomain) and
     144             :   // keyed by the base_name-prefixed deformation_gradient, so a given base_name resolves against its
     145             :   // own strain calculator.
     146        6762 :   _large_kinematics = hasGuaranteedMaterialProperty(_base_name + "deformation_gradient",
     147             :                                                     Guarantee::LARGE_KINEMATICS);
     148             : 
     149             :   // The deprecated local parameter must not silently disagree with the strain calculator.
     150       10143 :   if (isParamSetByUser("large_kinematics") &&
     151       11190 :       getParam<bool>("large_kinematics") != _large_kinematics)
     152           0 :     paramError("large_kinematics",
     153             :                "large_kinematics disagrees with the ComputeLagrangianStrain calculator (which "
     154             :                "computes ",
     155           0 :                _large_kinematics ? "large" : "small",
     156             :                " kinematics). large_kinematics is deprecated here; set it only on the strain "
     157             :                "calculator.");
     158        3381 : }
     159             : 
     160             : void
     161      612029 : LagrangianStressDivergenceBase::precalculateJacobian()
     162             : {
     163             :   // Skip if we are not doing stabilization
     164      612029 :   if (!_stabilize_strain)
     165             :     return;
     166             : 
     167             :   // We need the gradients of shape functions in the reference frame
     168       89061 :   _fe_problem.prepareShapes(_var.number(), _tid);
     169       89061 :   _avg_grad_trial[_alpha].resize(_phi.size());
     170       89061 :   precalculateJacobianDisplacement(_alpha);
     171             : }
     172             : 
     173             : void
     174     1072896 : LagrangianStressDivergenceBase::precalculateOffDiagJacobian(unsigned int jvar)
     175             : {
     176             :   // Skip if we are not doing stabilization
     177     1072896 :   if (!_stabilize_strain)
     178             :     return;
     179             : 
     180      614348 :   for (auto beta : make_range(_ndisp))
     181      457586 :     if (jvar == _disp_nums[beta])
     182             :     {
     183             :       // We need the gradients of shape functions in the reference frame
     184      154240 :       _fe_problem.prepareShapes(jvar, _tid);
     185      154240 :       _avg_grad_trial[beta].resize(_phi.size());
     186      154240 :       precalculateJacobianDisplacement(beta);
     187             :     }
     188             : }
     189             : 
     190             : Real
     191   272699670 : LagrangianStressDivergenceBase::computeQpJacobian()
     192             : {
     193   272699670 :   return computeQpJacobianDisplacement(_alpha, _alpha);
     194             : }
     195             : 
     196             : RankTwoTensor
     197      376704 : LagrangianStressDivergenceBase::deltaPK1NonLocalFBar(const RankTwoTensor & delta_F_avg) const
     198             : {
     199      376704 :   if (!_stabilize_strain)
     200           0 :     return RankTwoTensor();
     201             :   // `_d_nl_fbar` (the composed non-local operator) and F_ust^{-1}/det come from the stress and
     202             :   // strain materials as per-qp properties -- no kernel-side cache to refresh.
     203      376704 :   const RankTwoTensor delta_sigma_nl = (*_d_nl_fbar)[_qp] * delta_F_avg;
     204      376704 :   if (_large_kinematics)
     205      145280 :     return (*_F_ust_det)[_qp] * delta_sigma_nl * (*_F_ust_inv)[_qp].transpose();
     206      231424 :   return delta_sigma_nl;
     207             : }
     208             : 
     209             : Real
     210   515841558 : LagrangianStressDivergenceBase::computeQpOffDiagJacobian(unsigned int jvar)
     211             : {
     212             :   // Bail if jvar not coupled
     213   515841558 :   if (getJvarMap()[jvar] < 0)
     214             :     return 0.0;
     215             : 
     216             :   // Off diagonal terms for other displacements
     217  1031590913 :   for (auto beta : make_range(_ndisp))
     218  1030416193 :     if (jvar == _disp_nums[beta])
     219   514666838 :       return computeQpJacobianDisplacement(_alpha, beta);
     220             : 
     221             :   // Off diagonal temperature term due to eigenstrain
     222     1174720 :   if (_temperature && jvar == _temperature->number())
     223     1088832 :     return computeQpJacobianTemperature(mapJvarToCvar(jvar));
     224             : 
     225             :   // Off diagonal term due to weak plane stress
     226       85888 :   if (_out_of_plane_strain && jvar == _out_of_plane_strain->number())
     227       85888 :     return computeQpJacobianOutOfPlaneStrain();
     228             : 
     229             :   return 0;
     230             : }