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 "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("TensorMechanicsApp", ADComputeIsotropicElasticityTensorShell);
      21             : 
      22             : InputParameters
      23         160 : ADComputeIsotropicElasticityTensorShell::validParams()
      24             : {
      25         160 :   InputParameters params = Material::validParams();
      26         160 :   params.addClassDescription("Compute a plane stress isotropic elasticity tensor.");
      27         320 :   params.addRequiredRangeCheckedParam<Real>("poissons_ratio",
      28             :                                             "poissons_ratio >= -1.0 & poissons_ratio < 0.5",
      29             :                                             "Poisson's ratio for the material.");
      30         320 :   params.addRequiredRangeCheckedParam<Real>(
      31             :       "youngs_modulus", "youngs_modulus > 0.0", "Young's modulus of the material.");
      32         320 :   params.addRequiredParam<std::string>("through_thickness_order",
      33             :                                        "Quadrature order in out of plane direction");
      34         160 :   return params;
      35           0 : }
      36             : 
      37         120 : ADComputeIsotropicElasticityTensorShell::ADComputeIsotropicElasticityTensorShell(
      38         120 :     const InputParameters & parameters)
      39             :   : Material(parameters),
      40         120 :     _poissons_ratio(getParam<Real>("poissons_ratio")),
      41         360 :     _youngs_modulus(getParam<Real>("youngs_modulus"))
      42             : {
      43         120 :   _Cijkl.fillSymmetricIsotropicEandNu(_youngs_modulus, _poissons_ratio);
      44             : 
      45             :   // correction for plane stress
      46         120 :   _Cijkl(0, 0, 0, 0) = _youngs_modulus / (1.0 - _poissons_ratio * _poissons_ratio);
      47         120 :   _Cijkl(1, 1, 1, 1) = _Cijkl(0, 0, 0, 0);
      48         120 :   _Cijkl(0, 0, 1, 1) = _Cijkl(0, 0, 0, 0) * _poissons_ratio;
      49         120 :   _Cijkl(1, 1, 0, 0) = _Cijkl(0, 0, 1, 1);
      50         120 :   _Cijkl(0, 0, 2, 2) = 0.0;
      51         120 :   _Cijkl(1, 1, 2, 2) = 0.0;
      52         120 :   _Cijkl(2, 2, 2, 2) = 0.0;
      53         120 :   _Cijkl(2, 2, 0, 0) = 0.0;
      54         120 :   _Cijkl(2, 2, 1, 1) = 0.0;
      55             : 
      56             :   // get number of quadrature points along thickness based on order
      57             :   std::unique_ptr<QGauss> t_qrule = std::make_unique<QGauss>(
      58         240 :       1, Utility::string_to_enum<Order>(getParam<std::string>("through_thickness_order")));
      59         120 :   _t_points = t_qrule->get_points();
      60         120 :   _elasticity_tensor.resize(_t_points.size());
      61         120 :   _ge.resize(_t_points.size());
      62         360 :   for (unsigned int t = 0; t < _t_points.size(); ++t)
      63             :   {
      64         240 :     _elasticity_tensor[t] =
      65         240 :         &declareADProperty<RankFourTensor>("elasticity_tensor_t_points_" + std::to_string(t));
      66         720 :     _ge[t] = &getADMaterialProperty<RankTwoTensor>("ge_t_points_" + std::to_string(t));
      67             :   }
      68         120 : }
      69             : 
      70             : void
      71       19880 : ADComputeIsotropicElasticityTensorShell::computeQpProperties()
      72             : {
      73       59640 :   for (unsigned int t = 0; t < _t_points.size(); ++t)
      74             :   {
      75       39760 :     (*_elasticity_tensor[t])[_qp].zero();
      76             :     // compute contravariant elasticity tensor
      77      159040 :     for (unsigned int i = 0; i < 3; ++i)
      78      477120 :       for (unsigned int j = 0; j < 3; ++j)
      79     1431360 :         for (unsigned int k = 0; k < 3; ++k)
      80     4294080 :           for (unsigned int l = 0; l < 3; ++l)
      81    12882240 :             for (unsigned int m = 0; m < 3; ++m)
      82    38646720 :               for (unsigned int n = 0; n < 3; ++n)
      83   115940160 :                 for (unsigned int o = 0; o < 3; ++o)
      84   347820480 :                   for (unsigned int p = 0; p < 3; ++p)
      85   260865360 :                     (*_elasticity_tensor[t])[_qp](i, j, k, l) +=
      86   260865360 :                         (*_ge[t])[_qp](i, m) * (*_ge[t])[_qp](j, n) * (*_ge[t])[_qp](k, o) *
      87   521730720 :                         (*_ge[t])[_qp](l, p) * _Cijkl(m, n, o, p);
      88             :   }
      89       19880 : }