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 "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("SolidMechanicsApp", ADStressDivergenceShell);
      28             : 
      29             : InputParameters
      30        1400 : ADStressDivergenceShell::validParams()
      31             : {
      32        1400 :   InputParameters params = ADKernel::validParams();
      33        1400 :   params.addClassDescription("Quasi-static stress divergence kernel for Shell element");
      34        2800 :   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        2800 :   params.addRequiredParam<std::string>("through_thickness_order",
      41             :                                        "Quadrature order in out of plane direction");
      42        2800 :   params.addParam<bool>(
      43        2800 :       "large_strain", false, "Set to true to turn on finite strain calculations.");
      44        1400 :   params.set<bool>("use_displaced_mesh") = false;
      45        1400 :   return params;
      46           0 : }
      47             : 
      48         700 : ADStressDivergenceShell::ADStressDivergenceShell(const InputParameters & parameters)
      49             :   : ADKernel(parameters),
      50         700 :     _component(getParam<unsigned int>("component")),
      51        2100 :     _large_strain(getParam<bool>("large_strain"))
      52             : {
      53         700 :   _t_qrule = std::make_unique<libMesh::QGauss>(
      54        2100 :       1, Utility::string_to_enum<Order>(getParam<std::string>("through_thickness_order")));
      55         700 :   _t_weights = _t_qrule->get_weights();
      56             : 
      57         700 :   _stress.resize(_t_weights.size());
      58         700 :   _stress_old.resize(_t_weights.size());
      59         700 :   _B_mat.resize(_t_weights.size());
      60         700 :   if (_large_strain)
      61          60 :     _B_nl.resize(_t_weights.size());
      62         700 :   _J_map.resize(_t_weights.size());
      63             : 
      64        2100 :   for (unsigned int i = 0; i < _t_weights.size(); ++i)
      65             :   {
      66        2800 :     _stress[i] = &getADMaterialProperty<RankTwoTensor>("stress_t_points_" + std::to_string(i));
      67        1400 :     _stress_old[i] =
      68        1400 :         &getMaterialPropertyOldByName<RankTwoTensor>("stress_t_points_" + std::to_string(i));
      69        2800 :     _B_mat[i] = &getADMaterialProperty<DenseMatrix<Real>>("B_t_points_" + std::to_string(i));
      70        1400 :     if (_large_strain)
      71         360 :       _B_nl[i] = &getADMaterialProperty<DenseMatrix<Real>>("B_nl_t_points_" + std::to_string(i));
      72             : 
      73        4200 :     _J_map[i] = &getADMaterialProperty<Real>("J_mapping_t_points_" + std::to_string(i));
      74             :   }
      75         700 : }
      76             : 
      77             : ADReal
      78     1675840 : ADStressDivergenceShell::computeQpResidual()
      79             : {
      80     1675840 :   _q_weights = _qrule->get_weights();
      81     1675840 :   ADReal residual = 0.0;
      82     1675840 :   ADReal residual1 = 0.0;
      83     5027520 :   for (_qp_z = 0; _qp_z < _t_weights.size(); ++_qp_z)
      84             :   {
      85     3351680 :     residual1 = (*_stress[_qp_z])[_qp](0, 0) * (*_B_mat[_qp_z])[_qp](0, _i + _component * 4) +
      86     3351680 :                 (*_stress[_qp_z])[_qp](1, 1) * (*_B_mat[_qp_z])[_qp](1, _i + _component * 4) +
      87     3351680 :                 2.0 * (*_stress[_qp_z])[_qp](0, 1) * (*_B_mat[_qp_z])[_qp](2, _i + _component * 4) +
      88     3351680 :                 2.0 * (*_stress[_qp_z])[_qp](0, 2) * (*_B_mat[_qp_z])[_qp](3, _i + _component * 4) +
      89     6703360 :                 2.0 * (*_stress[_qp_z])[_qp](1, 2) * (*_B_mat[_qp_z])[_qp](4, _i + _component * 4);
      90             : 
      91     3351680 :     if (_large_strain)
      92             :       residual1 +=
      93      310400 :           ((*_stress_old[_qp_z])[_qp](0, 0) * (*_B_nl[_qp_z])[_qp](0, _i + _component * 4) +
      94      310400 :            (*_stress_old[_qp_z])[_qp](1, 1) * (*_B_nl[_qp_z])[_qp](1, _i + _component * 4) +
      95      310400 :            2.0 * (*_stress_old[_qp_z])[_qp](0, 1) * (*_B_nl[_qp_z])[_qp](2, _i + _component * 4) +
      96      310400 :            2.0 * (*_stress_old[_qp_z])[_qp](0, 2) * (*_B_nl[_qp_z])[_qp](3, _i + _component * 4) +
      97      620800 :            2.0 * (*_stress_old[_qp_z])[_qp](1, 2) * (*_B_nl[_qp_z])[_qp](4, _i + _component * 4));
      98             : 
      99     3351680 :     residual += residual1 * (*_J_map[_qp_z])[_qp] * _q_weights[_qp] * _t_weights[_qp_z] /
     100     6703360 :                 (_ad_JxW[_qp] * _ad_coord[_qp]);
     101             :   }
     102     1675840 :   return residual;
     103             : }