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         644 : ADComputeIsotropicElasticityTensorShell::validParams()
      24             : {
      25         644 :   InputParameters params = Material::validParams();
      26         644 :   params.addClassDescription("Compute a plane stress isotropic elasticity tensor.");
      27        1288 :   params.addRequiredRangeCheckedParam<Real>("poissons_ratio",
      28             :                                             "poissons_ratio >= -1.0 & poissons_ratio < 0.5",
      29             :                                             "Poisson's ratio for the material.");
      30        1288 :   params.addRequiredRangeCheckedParam<Real>(
      31             :       "youngs_modulus", "youngs_modulus > 0.0", "Young's modulus of the material.");
      32        1288 :   params.addRequiredParam<std::string>("through_thickness_order",
      33             :                                        "Quadrature order in out of plane direction");
      34         644 :   return params;
      35           0 : }
      36             : 
      37         483 : ADComputeIsotropicElasticityTensorShell::ADComputeIsotropicElasticityTensorShell(
      38         483 :     const InputParameters & parameters)
      39             :   : Material(parameters),
      40         483 :     _poissons_ratio(getParam<Real>("poissons_ratio")),
      41        1449 :     _youngs_modulus(getParam<Real>("youngs_modulus"))
      42             : {
      43         483 :   _Cijkl.fillSymmetricIsotropicEandNu(_youngs_modulus, _poissons_ratio);
      44             : 
      45             :   // correction for plane stress
      46         483 :   _Cijkl(0, 0, 0, 0) = _youngs_modulus / (1.0 - _poissons_ratio * _poissons_ratio);
      47         483 :   _Cijkl(1, 1, 1, 1) = _Cijkl(0, 0, 0, 0);
      48         483 :   _Cijkl(0, 0, 1, 1) = _Cijkl(0, 0, 0, 0) * _poissons_ratio;
      49         483 :   _Cijkl(1, 1, 0, 0) = _Cijkl(0, 0, 1, 1);
      50         483 :   _Cijkl(0, 0, 2, 2) = 0.0;
      51         483 :   _Cijkl(1, 1, 2, 2) = 0.0;
      52         483 :   _Cijkl(2, 2, 2, 2) = 0.0;
      53         483 :   _Cijkl(2, 2, 0, 0) = 0.0;
      54         483 :   _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         966 :       1, Utility::string_to_enum<Order>(getParam<std::string>("through_thickness_order")));
      59         483 :   _t_points = t_qrule->get_points();
      60         483 :   _elasticity_tensor.resize(_t_points.size());
      61         483 :   _ge.resize(_t_points.size());
      62        1449 :   for (unsigned int t = 0; t < _t_points.size(); ++t)
      63             :   {
      64         966 :     _elasticity_tensor[t] =
      65         966 :         &declareADProperty<RankFourTensor>("elasticity_tensor_t_points_" + std::to_string(t));
      66        2898 :     _ge[t] = &getADMaterialProperty<RankTwoTensor>("ge_t_points_" + std::to_string(t));
      67             :   }
      68         483 : }
      69             : 
      70             : void
      71      159480 : ADComputeIsotropicElasticityTensorShell::computeQpProperties()
      72             : {
      73      478440 :   for (unsigned int t = 0; t < _t_points.size(); ++t)
      74             :   {
      75      318960 :     (*_elasticity_tensor[t])[_qp].zero();
      76             :     // compute contravariant elasticity tensor
      77     1275840 :     for (unsigned int i = 0; i < 3; ++i)
      78     3827520 :       for (unsigned int j = 0; j < 3; ++j)
      79    11482560 :         for (unsigned int k = 0; k < 3; ++k)
      80    34447680 :           for (unsigned int l = 0; l < 3; ++l)
      81   103343040 :             for (unsigned int m = 0; m < 3; ++m)
      82   310029120 :               for (unsigned int n = 0; n < 3; ++n)
      83   930087360 :                 for (unsigned int o = 0; o < 3; ++o)
      84  2790262080 :                   for (unsigned int p = 0; p < 3; ++p)
      85  2092696560 :                     (*_elasticity_tensor[t])[_qp](i, j, k, l) +=
      86  2092696560 :                         (*_ge[t])[_qp](i, m) * (*_ge[t])[_qp](j, n) * (*_ge[t])[_qp](k, o) *
      87  4185393120 :                         (*_ge[t])[_qp](l, p) * _Cijkl(m, n, o, p);
      88             :   }
      89      159480 : }