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 "ADComputeShellStress.h"
      11             : #include "RankTwoTensor.h"
      12             : #include "RankFourTensor.h"
      13             : #include "ADComputeIsotropicElasticityTensorShell.h"
      14             : 
      15             : #include "libmesh/quadrature.h"
      16             : #include "libmesh/utility.h"
      17             : #include "libmesh/enum_quadrature_type.h"
      18             : #include "libmesh/fe_type.h"
      19             : #include "libmesh/string_to_enum.h"
      20             : #include "libmesh/quadrature_gauss.h"
      21             : 
      22             : registerMooseObject("SolidMechanicsApp", ADComputeShellStress);
      23             : 
      24             : InputParameters
      25         560 : ADComputeShellStress::validParams()
      26             : {
      27         560 :   InputParameters params = Material::validParams();
      28         560 :   params.addClassDescription("Compute in-plane stress using elasticity for shell");
      29        1120 :   params.addRequiredParam<std::string>("through_thickness_order",
      30             :                                        "Quadrature order in out of plane direction");
      31         560 :   return params;
      32           0 : }
      33             : 
      34         420 : ADComputeShellStress::ADComputeShellStress(const InputParameters & parameters)
      35         420 :   : Material(parameters)
      36             : {
      37             :   // get number of quadrature points along thickness based on order
      38             :   std::unique_ptr<libMesh::QGauss> t_qrule = std::make_unique<libMesh::QGauss>(
      39         840 :       1, Utility::string_to_enum<Order>(getParam<std::string>("through_thickness_order")));
      40         420 :   _t_points = t_qrule->get_points();
      41         420 :   _elasticity_tensor.resize(_t_points.size());
      42         420 :   _stress.resize(_t_points.size());
      43         420 :   _stress_old.resize(_t_points.size());
      44         420 :   _strain_increment.resize(_t_points.size());
      45         420 :   _local_transformation_matrix.resize(_t_points.size());
      46         420 :   _covariant_transformation_matrix.resize(_t_points.size());
      47         420 :   _global_stress.resize(_t_points.size());
      48         420 :   _local_stress.resize(_t_points.size());
      49        1260 :   for (unsigned int t = 0; t < _t_points.size(); ++t)
      50             :   {
      51         840 :     _elasticity_tensor[t] =
      52         840 :         &getADMaterialProperty<RankFourTensor>("elasticity_tensor_t_points_" + std::to_string(t));
      53         840 :     _stress[t] = &declareADProperty<RankTwoTensor>("stress_t_points_" + std::to_string(t));
      54         840 :     _stress_old[t] =
      55         840 :         &getMaterialPropertyOldByName<RankTwoTensor>("stress_t_points_" + std::to_string(t));
      56         840 :     _strain_increment[t] =
      57         840 :         &getADMaterialProperty<RankTwoTensor>("strain_increment_t_points_" + std::to_string(t));
      58             :     // rotation matrix and stress for output purposes only
      59         840 :     _local_transformation_matrix[t] =
      60         840 :         &getMaterialProperty<RankTwoTensor>("local_transformation_t_points_" + std::to_string(t));
      61         840 :     _covariant_transformation_matrix[t] = &getMaterialProperty<RankTwoTensor>(
      62         840 :         "covariant_transformation_t_points_" + std::to_string(t));
      63         840 :     _global_stress[t] =
      64         840 :         &declareProperty<RankTwoTensor>("global_stress_t_points_" + std::to_string(t));
      65         840 :     _local_stress[t] =
      66        1680 :         &declareProperty<RankTwoTensor>("local_stress_t_points_" + std::to_string(t));
      67             :   }
      68         420 : }
      69             : 
      70             : void
      71       47176 : ADComputeShellStress::initQpStatefulProperties()
      72             : {
      73             :   // initialize stress tensor to zero
      74      141528 :   for (unsigned int i = 0; i < _t_points.size(); ++i)
      75       94352 :     (*_stress[i])[_qp].zero();
      76       47176 : }
      77             : 
      78             : void
      79      111520 : ADComputeShellStress::computeQpProperties()
      80             : {
      81      334560 :   for (unsigned int i = 0; i < _t_points.size(); ++i)
      82             :   {
      83      223040 :     (*_stress[i])[_qp] =
      84      223040 :         (*_stress_old[i])[_qp] + (*_elasticity_tensor[i])[_qp] * (*_strain_increment[i])[_qp];
      85             : 
      86      892160 :     for (unsigned int ii = 0; ii < 3; ++ii)
      87     2676480 :       for (unsigned int jj = 0; jj < 3; ++jj)
      88     2007360 :         _unrotated_stress(ii, jj) = MetaPhysicL::raw_value((*_stress[i])[_qp](ii, jj));
      89      223040 :     (*_global_stress[i])[_qp] = (*_covariant_transformation_matrix[i])[_qp].transpose() *
      90      223040 :                                 _unrotated_stress * (*_covariant_transformation_matrix[i])[_qp];
      91      223040 :     (*_local_stress[i])[_qp] = (*_local_transformation_matrix[i])[_qp] * (*_global_stress[i])[_qp] *
      92      223040 :                                (*_local_transformation_matrix[i])[_qp].transpose();
      93             :   }
      94      111520 : }