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 "TotalLagrangianStressDivergenceBase.h"
      11             : 
      12             : template <class G>
      13             : InputParameters
      14        4540 : TotalLagrangianStressDivergenceBase<G>::baseParams()
      15             : {
      16        4540 :   InputParameters params = LagrangianStressDivergenceBase::validParams();
      17             :   // This kernel requires use_displaced_mesh to be off
      18        4540 :   params.suppressParameter<bool>("use_displaced_mesh");
      19        4540 :   return params;
      20           0 : }
      21             : 
      22             : template <class G>
      23        2270 : TotalLagrangianStressDivergenceBase<G>::TotalLagrangianStressDivergenceBase(
      24             :     const InputParameters & parameters)
      25             :   : LagrangianStressDivergenceBase(parameters),
      26        2270 :     _pk1(getMaterialPropertyByName<RankTwoTensor>(_base_name + "pk1_stress")),
      27        6810 :     _dpk1(getMaterialPropertyByName<RankFourTensor>(_base_name + "pk1_jacobian"))
      28             : {
      29        2270 : }
      30             : 
      31             : template <class G>
      32             : RankTwoTensor
      33  1526343944 : TotalLagrangianStressDivergenceBase<G>::gradTest(unsigned int component)
      34             : {
      35             :   // F-bar doesn't modify the test function
      36  1526343944 :   return G::gradOp(component, _grad_test[_i][_qp], _test[_i][_qp], _q_point[_qp]);
      37             : }
      38             : 
      39             : template <class G>
      40             : RankTwoTensor
      41   583538080 : TotalLagrangianStressDivergenceBase<G>::gradTrial(unsigned int component)
      42             : {
      43   583538080 :   return _stabilize_strain ? gradTrialStabilized(component) : gradTrialUnstabilized(component);
      44             : }
      45             : 
      46             : template <class G>
      47             : RankTwoTensor
      48   542123296 : TotalLagrangianStressDivergenceBase<G>::gradTrialUnstabilized(unsigned int component)
      49             : {
      50             :   // Without F-bar stabilization, simply return the gradient of the trial functions
      51   542123296 :   return G::gradOp(component, _grad_phi[_j][_qp], _phi[_j][_qp], _q_point[_qp]);
      52             : }
      53             : 
      54             : template <class G>
      55             : RankTwoTensor
      56    41414784 : TotalLagrangianStressDivergenceBase<G>::gradTrialStabilized(unsigned int component)
      57             : {
      58             :   // The base unstabilized trial function gradient
      59    41414784 :   const auto Gb = G::gradOp(component, _grad_phi[_j][_qp], _phi[_j][_qp], _q_point[_qp]);
      60             :   // The average trial function gradient
      61    41414784 :   const auto Ga = _avg_grad_trial[component][_j];
      62             : 
      63             :   // The F-bar stabilization depends on kinematics
      64    41414784 :   if (_large_kinematics)
      65             :   {
      66             :     // Horrible thing, see the documentation for how we get here
      67     9324032 :     const Real dratio = std::pow(_F_avg[_qp].det() / _F_ust[_qp].det(), 1.0 / 3.0);
      68     9324032 :     const Real fact = (_F_avg[_qp].inverse().transpose().doubleContraction(Ga) -
      69     9324032 :                        _F_ust[_qp].inverse().transpose().doubleContraction(Gb)) /
      70             :                       3.0;
      71     9324032 :     return dratio * (Gb + fact * _F_ust[_qp]);
      72             :   }
      73             : 
      74             :   // The small kinematics modification is linear
      75    32090752 :   return Gb + (Ga.trace() - Gb.trace()) / 3.0 * RankTwoTensor::Identity();
      76             : }
      77             : 
      78             : template <class G>
      79             : void
      80       97496 : TotalLagrangianStressDivergenceBase<G>::precalculateJacobianDisplacement(unsigned int component)
      81             : {
      82             :   // For total Lagrangian, the averaging is taken on the reference frame regardless of geometric
      83             :   // nonlinearity. Convenient!
      84      781880 :   for (auto j : make_range(_phi.size()))
      85      684384 :     _avg_grad_trial[component][j] = StabilizationUtils::elementAverage(
      86     5785376 :         [this, component, j](unsigned int qp)
      87     5100992 :         { return G::gradOp(component, _grad_phi[j][qp], _phi[j][qp], _q_point[qp]); },
      88             :         _JxW,
      89             :         _coord);
      90       97496 : }
      91             : 
      92             : template <class G>
      93             : Real
      94   945789352 : TotalLagrangianStressDivergenceBase<G>::computeQpResidual()
      95             : {
      96   945789352 :   return gradTest(_alpha).doubleContraction(_pk1[_qp]);
      97             : }
      98             : 
      99             : template <class G>
     100             : Real
     101   536551776 : TotalLagrangianStressDivergenceBase<G>::computeQpJacobianDisplacement(unsigned int alpha,
     102             :                                                                       unsigned int beta)
     103             : {
     104             :   // J_{alpha beta} = phi^alpha_{i, J} T_{iJkL} G^beta_{kL}
     105   536551776 :   return gradTest(alpha).doubleContraction(_dpk1[_qp] * gradTrial(beta));
     106             : }
     107             : 
     108             : template <class G>
     109             : Real
     110     1171968 : TotalLagrangianStressDivergenceBase<G>::computeQpJacobianTemperature(unsigned int cvar)
     111             : {
     112             :   usingTensorIndices(i_, j_, k_, l_);
     113             :   // Multiple eigenstrains may depend on the same coupled var
     114     1171968 :   RankTwoTensor total_deigen;
     115     2343936 :   for (const auto deigen_darg : _deigenstrain_dargs[cvar])
     116     1171968 :     total_deigen += (*deigen_darg)[_qp];
     117             : 
     118     1171968 :   const auto A = _f_inv[_qp].inverse();
     119     1171968 :   const auto B = _F_inv[_qp].inverse();
     120     1171968 :   const auto U = 0.5 * (A.template times<i_, k_, l_, j_>(B) + A.template times<i_, l_, k_, j_>(B));
     121             : 
     122     1171968 :   return -(_dpk1[_qp] * U * total_deigen).doubleContraction(gradTest(_alpha)) *
     123     1171968 :          _temperature->phi()[_j][_qp];
     124             : }
     125             : 
     126             : template <class G>
     127             : Real
     128      168960 : TotalLagrangianStressDivergenceBase<G>::computeQpJacobianOutOfPlaneStrain()
     129             : {
     130      168960 :   return _dpk1[_qp].contractionKl(2, 2, gradTest(_alpha)) * _out_of_plane_strain->phi()[_j][_qp];
     131             : }
     132             : 
     133             : template class TotalLagrangianStressDivergenceBase<GradientOperatorCartesian>;
     134             : template class TotalLagrangianStressDivergenceBase<GradientOperatorAxisymmetricCylindrical>;
     135             : template class TotalLagrangianStressDivergenceBase<GradientOperatorCentrosymmetricSpherical>;