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 "ADComputeIsotropicElasticityTensorShell.h"
      11             : #include "RankFourTensor.h"
      12             : 
      13             : #include "libmesh/quadrature.h"
      14             : #include "libmesh/utility.h"
      15             : #include "libmesh/enum_quadrature_type.h"
      16             : #include "libmesh/fe_type.h"
      17             : #include "libmesh/string_to_enum.h"
      18             : #include "libmesh/quadrature_gauss.h"
      19             : 
      20             : registerMooseObject("SolidMechanicsApp", ADComputeIsotropicElasticityTensorShell);
      21             : 
      22             : InputParameters
      23         560 : ADComputeIsotropicElasticityTensorShell::validParams()
      24             : {
      25         560 :   InputParameters params = Material::validParams();
      26         560 :   params.addClassDescription("Compute a plane stress isotropic elasticity tensor.");
      27        1120 :   params.addRequiredRangeCheckedParam<Real>("poissons_ratio",
      28             :                                             "poissons_ratio >= -1.0 & poissons_ratio < 0.5",
      29             :                                             "Poisson's ratio for the material.");
      30        1120 :   params.addRequiredRangeCheckedParam<Real>(
      31             :       "youngs_modulus", "youngs_modulus > 0.0", "Young's modulus of the material.");
      32        1120 :   params.addRequiredParam<std::string>("through_thickness_order",
      33             :                                        "Quadrature order in out of plane direction");
      34         560 :   return params;
      35           0 : }
      36             : 
      37         420 : ADComputeIsotropicElasticityTensorShell::ADComputeIsotropicElasticityTensorShell(
      38         420 :     const InputParameters & parameters)
      39             :   : Material(parameters),
      40         420 :     _poissons_ratio(getParam<Real>("poissons_ratio")),
      41        1260 :     _youngs_modulus(getParam<Real>("youngs_modulus"))
      42             : {
      43         420 :   _Cijkl.fillSymmetricIsotropicEandNu(_youngs_modulus, _poissons_ratio);
      44             : 
      45             :   // correction for plane stress
      46         420 :   _Cijkl(0, 0, 0, 0) = _youngs_modulus / (1.0 - _poissons_ratio * _poissons_ratio);
      47         420 :   _Cijkl(1, 1, 1, 1) = _Cijkl(0, 0, 0, 0);
      48         420 :   _Cijkl(0, 0, 1, 1) = _Cijkl(0, 0, 0, 0) * _poissons_ratio;
      49         420 :   _Cijkl(1, 1, 0, 0) = _Cijkl(0, 0, 1, 1);
      50         420 :   _Cijkl(0, 0, 2, 2) = 0.0;
      51         420 :   _Cijkl(1, 1, 2, 2) = 0.0;
      52         420 :   _Cijkl(2, 2, 2, 2) = 0.0;
      53         420 :   _Cijkl(2, 2, 0, 0) = 0.0;
      54         420 :   _Cijkl(2, 2, 1, 1) = 0.0;
      55             : 
      56             :   // get number of quadrature points along thickness based on order
      57             :   std::unique_ptr<libMesh::QGauss> t_qrule = std::make_unique<libMesh::QGauss>(
      58         840 :       1, Utility::string_to_enum<Order>(getParam<std::string>("through_thickness_order")));
      59         420 :   _t_points = t_qrule->get_points();
      60         420 :   _elasticity_tensor.resize(_t_points.size());
      61         420 :   _ge.resize(_t_points.size());
      62        1260 :   for (unsigned int t = 0; t < _t_points.size(); ++t)
      63             :   {
      64         840 :     _elasticity_tensor[t] =
      65         840 :         &declareADProperty<RankFourTensor>("elasticity_tensor_t_points_" + std::to_string(t));
      66        2520 :     _ge[t] = &getADMaterialProperty<RankTwoTensor>("ge_t_points_" + std::to_string(t));
      67             :   }
      68         420 : }
      69             : 
      70             : void
      71      111520 : ADComputeIsotropicElasticityTensorShell::computeQpProperties()
      72             : {
      73      334560 :   for (unsigned int t = 0; t < _t_points.size(); ++t)
      74             :   {
      75      223040 :     (*_elasticity_tensor[t])[_qp].zero();
      76             :     // compute contravariant elasticity tensor
      77      892160 :     for (unsigned int i = 0; i < 3; ++i)
      78     2676480 :       for (unsigned int j = 0; j < 3; ++j)
      79     8029440 :         for (unsigned int k = 0; k < 3; ++k)
      80    24088320 :           for (unsigned int l = 0; l < 3; ++l)
      81    72264960 :             for (unsigned int m = 0; m < 3; ++m)
      82   216794880 :               for (unsigned int n = 0; n < 3; ++n)
      83   650384640 :                 for (unsigned int o = 0; o < 3; ++o)
      84  1951153920 :                   for (unsigned int p = 0; p < 3; ++p)
      85  1463365440 :                     (*_elasticity_tensor[t])[_qp](i, j, k, l) +=
      86  1463365440 :                         (*_ge[t])[_qp](i, m) * (*_ge[t])[_qp](j, n) * (*_ge[t])[_qp](k, o) *
      87  2926730880 :                         (*_ge[t])[_qp](l, p) * _Cijkl(m, n, o, p);
      88             :   }
      89      111520 : }