#include <StressDivergence.h>

Inheritance diagram for StressDivergence:

Public Member Functions
	StressDivergence (const InputParameters &parameters)

Protected Member Functions
virtual void	computeResidual () override

virtual void	computeJacobian () override

virtual void	computeOffDiagJacobian (MooseVariableFEBase &jvar) override

virtual Real	computeQpResidual () override

virtual Real	computeQpJacobian () override

virtual Real	computeQpOffDiagJacobian (unsigned int jvar) override

Protected Attributes
const MaterialProperty< SymmTensor > &	_stress_older

const MaterialProperty< SymmTensor > &	_stress_old

const MaterialProperty< SymmTensor > &	_stress

const MaterialProperty< SymmElasticityTensor > &	_Jacobian_mult

const MaterialProperty< SymmTensor > &	_d_stress_dT

Private Attributes
const unsigned int	_component

const bool	_xdisp_coupled

const bool	_ydisp_coupled

const bool	_zdisp_coupled

const bool	_temp_coupled

const unsigned int	_xdisp_var

const unsigned int	_ydisp_var

const unsigned int	_zdisp_var

const unsigned int	_temp_var

const Real	_zeta

const Real	_alpha

std::vector< std::vector< Real > >	_avg_grad_test

std::vector< std::vector< Real > >	_avg_grad_phi

bool	_volumetric_locking_correction

Detailed Description

Definition at line 22 of file StressDivergence.h.

Constructor & Destructor Documentation

◆ StressDivergence()

StressDivergence::StressDivergence ( const InputParameters & parameters )

Definition at line 50 of file StressDivergence.C.

   : Kernel(parameters),
     _stress_older(getMaterialPropertyOlder<SymmTensor>(
         "stress" + getParam<std::string>("appended_property_name"))),
     _stress_old(getMaterialPropertyOld<SymmTensor>(
         "stress" + getParam<std::string>("appended_property_name"))),
     _stress(getMaterialProperty<SymmTensor>("stress" +
                                             getParam<std::string>("appended_property_name"))),
     _Jacobian_mult(getMaterialProperty<SymmElasticityTensor>(
         "Jacobian_mult" + getParam<std::string>("appended_property_name"))),
     _d_stress_dT(getMaterialProperty<SymmTensor>("d_stress_dT" +
                                                  getParam<std::string>("appended_property_name"))),
     _component(getParam<unsigned int>("component")),
     _xdisp_coupled(isCoupled("disp_x")),
     _ydisp_coupled(isCoupled("disp_y")),
     _zdisp_coupled(isCoupled("disp_z")),
     _temp_coupled(isCoupled("temp")),
     _xdisp_var(_xdisp_coupled ? coupled("disp_x") : 0),
     _ydisp_var(_ydisp_coupled ? coupled("disp_y") : 0),
     _zdisp_var(_zdisp_coupled ? coupled("disp_z") : 0),
     _temp_var(_temp_coupled ? coupled("temp") : 0),
     _zeta(getParam<Real>("zeta")),
     _alpha(getParam<Real>("alpha")),
     _avg_grad_test(_test.size(), std::vector<Real>(3, 0.0)),
     _avg_grad_phi(_phi.size(), std::vector<Real>(3, 0.0)),
     _volumetric_locking_correction(getParam<bool>("volumetric_locking_correction"))
 {
   mooseDeprecated(name(), ": StressDivergence is deprecated.\n\
 The solid_mechanics module will be removed from MOOSE on July 31, 2020.\n\
 Please update your input files to utilize the tensor_mechanics equivalents of\n\
 models based on solid_mechanics. A detailed migration guide that was developed\n\
 for BISON, but which is generally applicable to any MOOSE model is available at:\n\
 https://mooseframework.org/bison/tutorials/mechanics_conversion/overview.html");
 }

Member Function Documentation

◆ computeJacobian()

void StressDivergence::computeJacobian ( )

overrideprotectedvirtual

Definition at line 152 of file StressDivergence.C.

 {
   prepareMatrixTag(_assembly, _var.number(), _var.number());
  
   if (_volumetric_locking_correction)
   {
     // calculate volume averaged value of shape function derivative
     _avg_grad_test.resize(_test.size());
     for (_i = 0; _i < _test.size(); _i++)
     {
       _avg_grad_test[_i].resize(3);
       _avg_grad_test[_i][_component] = 0.0;
       for (_qp = 0; _qp < _qrule->n_points(); _qp++)
         _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
  
       _avg_grad_test[_i][_component] /= _current_elem_volume;
     }
  
     _avg_grad_phi.resize(_phi.size());
     for (_i = 0; _i < _phi.size(); _i++)
     {
       _avg_grad_phi[_i].resize(3);
       for (unsigned int component = 0; component < _mesh.dimension(); component++)
       {
         _avg_grad_phi[_i][component] = 0.0;
         for (_qp = 0; _qp < _qrule->n_points(); _qp++)
           _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
  
         _avg_grad_phi[_i][component] /= _current_elem_volume;
       }
     }
   }
  
   for (_i = 0; _i < _test.size(); _i++)
     for (_j = 0; _j < _phi.size(); _j++)
       for (_qp = 0; _qp < _qrule->n_points(); _qp++)
         _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJacobian();
  
   accumulateTaggedLocalMatrix();
  
   if (_has_diag_save_in)
   {
     unsigned int rows = _local_ke.m();
     DenseVector<Number> diag(rows);
     for (unsigned int i = 0; i < rows; i++)
       diag(i) = _local_ke(i, i);
  
     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     for (const auto & var : _diag_save_in)
       var->sys().solution().add_vector(diag, var->dofIndices());
   }
 }

Referenced by computeOffDiagJacobian().

◆ computeOffDiagJacobian()

void StressDivergence::computeOffDiagJacobian ( MooseVariableFEBase & jvar )

overrideprotectedvirtual

Definition at line 263 of file StressDivergence.C.

 {
   size_t jvar_num = jvar.number();
   if (jvar_num == _var.number())
     computeJacobian();
   else
   {
     // This (undisplaced) jvar could potentially yield the wrong phi size if this object is acting
     // on the displaced mesh
     auto phi_size = _sys.getVariable(_tid, jvar.number()).dofIndices().size();
  
     if (_volumetric_locking_correction)
     {
       // calculate volume averaged value of shape function derivative
       _avg_grad_test.resize(_test.size());
       for (_i = 0; _i < _test.size(); _i++)
       {
         _avg_grad_test[_i].resize(3);
         _avg_grad_test[_i][_component] = 0.0;
         for (_qp = 0; _qp < _qrule->n_points(); _qp++)
           _avg_grad_test[_i][_component] +=
               _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
  
         _avg_grad_test[_i][_component] /= _current_elem_volume;
       }
  
       _avg_grad_phi.resize(phi_size);
       for (_i = 0; _i < phi_size; _i++)
       {
         _avg_grad_phi[_i].resize(3);
         for (unsigned int component = 0; component < _mesh.dimension(); component++)
         {
           _avg_grad_phi[_i][component] = 0.0;
           for (_qp = 0; _qp < _qrule->n_points(); _qp++)
             _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
  
           _avg_grad_phi[_i][component] /= _current_elem_volume;
         }
       }
     }
  
     prepareMatrixTag(_assembly, _var.number(), jvar_num);
  
     for (_i = 0; _i < _test.size(); _i++)
       for (_j = 0; _j < phi_size; _j++)
         for (_qp = 0; _qp < _qrule->n_points(); _qp++)
           _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(jvar_num);
  
     accumulateTaggedLocalMatrix();
   }
 }

◆ computeQpJacobian()

Real StressDivergence::computeQpJacobian ( )

overrideprotectedvirtual

Definition at line 206 of file StressDivergence.C.

 {
   Real sum_C3x3 = _Jacobian_mult[_qp].sum_3x3();
   RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum_3x1();
  
   Real jacobian = 0.0;
   // B^T_i * C * B_j
   jacobian += _Jacobian_mult[_qp].stiffness(
       _component, _component, _grad_test[_i][_qp], _grad_phi[_j][_qp]); // B^T_i * C *B_j
  
   if (_volumetric_locking_correction)
   {
     // jacobian = Bbar^T_i * C * Bbar_j where Bbar = B + Bvol
     // jacobian = B^T_i * C * B_j + Bvol^T_i * C * Bvol_j +  Bvol^T_i * C * B_j + B^T_i * C * Bvol_j
  
     // Bvol^T_i * C * Bvol_j
     jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
                 (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 9.0;
  
     // B^T_i * C * Bvol_j
     jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
                 (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 3.0;
  
     // Bvol^T_i * C * B_j
     SymmTensor phi;
     if (_component == 0)
     {
       phi.xx() = _grad_phi[_j][_qp](0);
       phi.xy() = _grad_phi[_j][_qp](1);
       phi.xz() = _grad_phi[_j][_qp](2);
     }
     else if (_component == 1)
     {
       phi.yy() = _grad_phi[_j][_qp](1);
       phi.xy() = _grad_phi[_j][_qp](0);
       phi.yz() = _grad_phi[_j][_qp](2);
     }
     else if (_component == 2)
     {
       phi.zz() = _grad_phi[_j][_qp](2);
       phi.xz() = _grad_phi[_j][_qp](0);
       phi.yz() = _grad_phi[_j][_qp](1);
     }
  
     SymmTensor tmp(_Jacobian_mult[_qp] * phi);
  
     jacobian += (tmp.xx() + tmp.yy() + tmp.zz()) *
                 (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
   }
  
   if (_dt > 0)
     return jacobian * (1 + _alpha + _zeta / _dt);
   else
     return jacobian;
 }

Referenced by computeJacobian().

◆ computeQpOffDiagJacobian()

Real StressDivergence::computeQpOffDiagJacobian ( unsigned int jvar )

overrideprotectedvirtual

Definition at line 316 of file StressDivergence.C.

 {
   unsigned int coupled_component = 0;
  
   bool active = false;
  
   if (_xdisp_coupled && jvar == _xdisp_var)
   {
     coupled_component = 0;
     active = true;
   }
   else if (_ydisp_coupled && jvar == _ydisp_var)
   {
     coupled_component = 1;
     active = true;
   }
   else if (_zdisp_coupled && jvar == _zdisp_var)
   {
     coupled_component = 2;
     active = true;
   }
  
   if (active)
   {
     Real sum_C3x3 = _Jacobian_mult[_qp].sum_3x3();
     RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum_3x1();
     Real jacobian = 0.0;
     jacobian += _Jacobian_mult[_qp].stiffness(
         _component, coupled_component, _grad_test[_i][_qp], _grad_phi[_j][_qp]); // B^T_i * C *B_j
  
     if (_volumetric_locking_correction)
     {
       // jacobian = Bbar^T_i * C * Bbar_j where Bbar = B + Bvol
       // jacobian = B^T_i * C * B_j + Bvol^T_i * C * Bvol_j +  Bvol^T_i * C * B_j + B^T_i * C *
       // Bvol_j
  
       // Bvol^T_i * C * Bvol_j
       jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
                   (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
                   9.0;
  
       // B^T_i * C * Bvol_j
       jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
                   (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
                   3.0;
  
       // Bvol^T_i * C * B_i
       SymmTensor phi;
       if (coupled_component == 0)
       {
         phi.xx() = _grad_phi[_j][_qp](0);
         phi.xy() = _grad_phi[_j][_qp](1);
         phi.xz() = _grad_phi[_j][_qp](2);
       }
       else if (coupled_component == 1)
       {
         phi.yy() = _grad_phi[_j][_qp](1);
         phi.xy() = _grad_phi[_j][_qp](0);
         phi.yz() = _grad_phi[_j][_qp](2);
       }
       else if (coupled_component == 2)
       {
         phi.zz() = _grad_phi[_j][_qp](2);
         phi.xz() = _grad_phi[_j][_qp](0);
         phi.yz() = _grad_phi[_j][_qp](1);
       }
  
       SymmTensor tmp(_Jacobian_mult[_qp] * phi);
  
       jacobian += (tmp.xx() + tmp.yy() + tmp.zz()) *
                   (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
     }
  
     if (_dt > 0)
       return jacobian * (1 + _alpha + _zeta / _dt);
     else
       return jacobian;
   }
  
   if (_temp_coupled && jvar == _temp_var)
     return _d_stress_dT[_qp].rowDot(_component, _grad_test[_i][_qp]) * _phi[_j][_qp];
  
   return 0;
 }

Referenced by computeOffDiagJacobian().

◆ computeQpResidual()

Real StressDivergence::computeQpResidual ( )

overrideprotectedvirtual

Definition at line 121 of file StressDivergence.C.

 {
   Real residual = 0.0;
   if ((_dt > 0) && ((_zeta != 0) || (_alpha != 0)))
   {
     residual = _stress[_qp].rowDot(_component, _grad_test[_i][_qp]) *
                    (1 + _alpha + (1 + _alpha) * _zeta / _dt) -
                (_alpha + (1 + 2 * _alpha) * _zeta / _dt) *
                    _stress_old[_qp].rowDot(_component, _grad_test[_i][_qp]) +
                (_alpha * _zeta / _dt) * _stress_older[_qp].rowDot(_component, _grad_test[_i][_qp]);
  
     // volumetric locking correction
     if (_volumetric_locking_correction)
       residual += (_stress[_qp].trace() * (1 + _alpha + (1 + _alpha) * _zeta / _dt) -
                    (_alpha + (1 + 2 * _alpha) * _zeta / _dt) * _stress_old[_qp].trace() +
                    (_alpha * _zeta / _dt) * _stress_older[_qp].trace()) /
                   3.0 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
   }
   else
   {
     residual += _stress[_qp].rowDot(_component, _grad_test[_i][_qp]);
  
     // volumetric locking correction
     if (_volumetric_locking_correction)
       residual += _stress[_qp].trace() / 3.0 *
                   (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
   }
   return residual;
 }

Referenced by computeResidual().

◆ computeResidual()

void StressDivergence::computeResidual ( )

overrideprotectedvirtual

Definition at line 86 of file StressDivergence.C.

 {
   prepareVectorTag(_assembly, _var.number());
  
   if (_volumetric_locking_correction)
   {
     // calculate volume averaged value of shape function derivative
     _avg_grad_test.resize(_test.size());
     for (_i = 0; _i < _test.size(); _i++)
     {
       _avg_grad_test[_i].resize(3);
       _avg_grad_test[_i][_component] = 0.0;
       for (_qp = 0; _qp < _qrule->n_points(); _qp++)
         _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
  
       _avg_grad_test[_i][_component] /= _current_elem_volume;
     }
   }
  
   precalculateResidual();
   for (_i = 0; _i < _test.size(); _i++)
     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
       _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
  
   accumulateTaggedLocalResidual();
  
   if (_has_save_in)
   {
     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     for (const auto & var : _save_in)
       var->sys().solution().add_vector(_local_re, var->dofIndices());
   }
 }