Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       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 "ADStressDivergenceShell.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Assembly.h"
      14             : #include "Material.h"
      15             : #include "MooseVariable.h"
      16             : #include "SystemBase.h"
      17             : #include "RankTwoTensor.h"
      18             : #include "NonlinearSystem.h"
      19             : #include "MooseMesh.h"
      20             : #include "ArbitraryQuadrature.h"
      21             : #include "DenseMatrix.h"
      22             : 
      23             : #include "libmesh/quadrature.h"
      24             : #include "libmesh/enum_quadrature_type.h"
      25             : #include "libmesh/string_to_enum.h"
      26             : 
      27             : registerMooseObject("TensorMechanicsApp", ADStressDivergenceShell);
      28             : 
      29             : InputParameters
      30         400 : ADStressDivergenceShell::validParams()
      31             : {
      32         400 :   InputParameters params = ADKernel::validParams();
      33         400 :   params.addClassDescription("Quasi-static stress divergence kernel for Shell element");
      34         800 :   params.addRequiredRangeCheckedParam<unsigned int>(
      35             :       "component",
      36             :       "component < 5",
      37             :       "An integer corresponding to the degree of freedom "
      38             :       "this kernel acts on. (0 for disp_x, "
      39             :       "1 for disp_y, 2 for disp_z, 3 for rot_x, 4 for rot_y)");
      40         800 :   params.addRequiredParam<std::string>("through_thickness_order",
      41             :                                        "Quadrature order in out of plane direction");
      42         800 :   params.addParam<bool>(
      43         800 :       "large_strain", false, "Set to true to turn on finite strain calculations.");
      44         400 :   params.set<bool>("use_displaced_mesh") = false;
      45         400 :   return params;
      46           0 : }
      47             : 
      48         200 : ADStressDivergenceShell::ADStressDivergenceShell(const InputParameters & parameters)
      49             :   : ADKernel(parameters),
      50         200 :     _component(getParam<unsigned int>("component")),
      51         600 :     _large_strain(getParam<bool>("large_strain"))
      52             : {
      53         200 :   _t_qrule = std::make_unique<QGauss>(
      54         600 :       1, Utility::string_to_enum<Order>(getParam<std::string>("through_thickness_order")));
      55         200 :   _t_weights = _t_qrule->get_weights();
      56             : 
      57         200 :   _stress.resize(_t_weights.size());
      58         200 :   _stress_old.resize(_t_weights.size());
      59         200 :   _B_mat.resize(_t_weights.size());
      60         200 :   if (_large_strain)
      61          30 :     _B_nl.resize(_t_weights.size());
      62         200 :   _J_map.resize(_t_weights.size());
      63             : 
      64         600 :   for (unsigned int i = 0; i < _t_weights.size(); ++i)
      65             :   {
      66         800 :     _stress[i] = &getADMaterialProperty<RankTwoTensor>("stress_t_points_" + std::to_string(i));
      67         400 :     _stress_old[i] =
      68         400 :         &getMaterialPropertyOldByName<RankTwoTensor>("stress_t_points_" + std::to_string(i));
      69         800 :     _B_mat[i] = &getADMaterialProperty<DenseMatrix<Real>>("B_t_points_" + std::to_string(i));
      70         400 :     if (_large_strain)
      71         180 :       _B_nl[i] = &getADMaterialProperty<DenseMatrix<Real>>("B_nl_t_points_" + std::to_string(i));
      72             : 
      73        1200 :     _J_map[i] = &getADMaterialProperty<Real>("J_mapping_t_points_" + std::to_string(i));
      74             :   }
      75         200 : }
      76             : 
      77             : ADReal
      78      323200 : ADStressDivergenceShell::computeQpResidual()
      79             : {
      80      323200 :   _q_weights = _qrule->get_weights();
      81      323200 :   ADReal residual = 0.0;
      82      323200 :   ADReal residual1 = 0.0;
      83      969600 :   for (_qp_z = 0; _qp_z < _t_weights.size(); ++_qp_z)
      84             :   {
      85      646400 :     residual1 = (*_stress[_qp_z])[_qp](0, 0) * (*_B_mat[_qp_z])[_qp](0, _i + _component * 4) +
      86      646400 :                 (*_stress[_qp_z])[_qp](1, 1) * (*_B_mat[_qp_z])[_qp](1, _i + _component * 4) +
      87      646400 :                 2.0 * (*_stress[_qp_z])[_qp](0, 1) * (*_B_mat[_qp_z])[_qp](2, _i + _component * 4) +
      88      646400 :                 2.0 * (*_stress[_qp_z])[_qp](0, 2) * (*_B_mat[_qp_z])[_qp](3, _i + _component * 4) +
      89     1292800 :                 2.0 * (*_stress[_qp_z])[_qp](1, 2) * (*_B_mat[_qp_z])[_qp](4, _i + _component * 4);
      90             : 
      91      646400 :     if (_large_strain)
      92             :       residual1 +=
      93      294400 :           ((*_stress_old[_qp_z])[_qp](0, 0) * (*_B_nl[_qp_z])[_qp](0, _i + _component * 4) +
      94      294400 :            (*_stress_old[_qp_z])[_qp](1, 1) * (*_B_nl[_qp_z])[_qp](1, _i + _component * 4) +
      95      294400 :            2.0 * (*_stress_old[_qp_z])[_qp](0, 1) * (*_B_nl[_qp_z])[_qp](2, _i + _component * 4) +
      96      294400 :            2.0 * (*_stress_old[_qp_z])[_qp](0, 2) * (*_B_nl[_qp_z])[_qp](3, _i + _component * 4) +
      97      588800 :            2.0 * (*_stress_old[_qp_z])[_qp](1, 2) * (*_B_nl[_qp_z])[_qp](4, _i + _component * 4));
      98             : 
      99      646400 :     residual += residual1 * (*_J_map[_qp_z])[_qp] * _q_weights[_qp] * _t_weights[_qp_z] /
     100     1292800 :                 (_ad_JxW[_qp] * _ad_coord[_qp]);
     101             :   }
     102      323200 :   return residual;
     103             : }