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 "ComputeReducedOrderEigenstrain.h"
      11             : 
      12             : #include "MooseMesh.h"
      13             : #include "libmesh/quadrature.h"
      14             : 
      15             : registerMooseObject("SolidMechanicsApp", ComputeReducedOrderEigenstrain);
      16             : 
      17             : InputParameters
      18          96 : ComputeReducedOrderEigenstrain::validParams()
      19             : {
      20          96 :   InputParameters params = ComputeEigenstrainBase::validParams();
      21          96 :   params.addClassDescription("accepts eigenstrains and computes a reduced order eigenstrain for "
      22             :                              "consistency in the order of strain and eigenstrains.");
      23         192 :   params.addRequiredParam<std::vector<MaterialPropertyName>>(
      24             :       "input_eigenstrain_names", "List of eigenstrains to be applied in this strain calculation");
      25          96 :   return params;
      26           0 : }
      27             : 
      28          72 : ComputeReducedOrderEigenstrain::ComputeReducedOrderEigenstrain(const InputParameters & parameters)
      29             :   : ComputeEigenstrainBase(parameters),
      30          72 :     _input_eigenstrain_names(
      31             :         getParam<std::vector<MaterialPropertyName>>("input_eigenstrain_names")),
      32          72 :     _eigenstrains(_input_eigenstrain_names.size()),
      33          72 :     _subproblem(*parameters.get<SubProblem *>("_subproblem")),
      34          72 :     _ncols(1 + _subproblem.mesh().dimension()),
      35          72 :     _second_order(_subproblem.mesh().hasSecondOrderElements()),
      36             :     _eigsum(),
      37          72 :     _A(),
      38          72 :     _b(6),
      39          72 :     _AT(),
      40          72 :     _ATb(_ncols),
      41          72 :     _x(6, DenseVector<Real>(_ncols)),
      42          72 :     _vals(6),
      43         144 :     _adjusted_eigenstrain()
      44             : {
      45         144 :   for (unsigned int i = 0; i < _eigenstrains.size(); ++i)
      46             :   {
      47          72 :     _input_eigenstrain_names[i] = _base_name + _input_eigenstrain_names[i];
      48          72 :     _eigenstrains[i] = &getMaterialProperty<RankTwoTensor>(_input_eigenstrain_names[i]);
      49             :   }
      50          72 : }
      51             : 
      52             : void
      53         576 : ComputeReducedOrderEigenstrain::initQpStatefulProperties()
      54             : {
      55         576 :   _eigenstrain[_qp].zero();
      56         576 : }
      57             : 
      58             : void
      59        1720 : ComputeReducedOrderEigenstrain::computeProperties()
      60             : {
      61        1720 :   sumEigenstrain();
      62             : 
      63        1720 :   prepareEigenstrain();
      64             : 
      65        1720 :   Material::computeProperties();
      66        1720 : }
      67             : 
      68             : void
      69       15480 : ComputeReducedOrderEigenstrain::computeQpEigenstrain()
      70             : {
      71       15480 :   if (_second_order)
      72             :   {
      73       13608 :     for (unsigned i = 0; i < 6; ++i)
      74             :     {
      75       11664 :       _vals[i] = _x[i](0);
      76       34992 :       for (unsigned j = 1; j < _ncols; ++j)
      77       23328 :         _vals[i] += _x[i](j) * _q_point[_qp](j - 1);
      78             :     }
      79        1944 :     _adjusted_eigenstrain.fillFromInputVector(_vals);
      80             :   }
      81       15480 :   _eigenstrain[_qp] = _adjusted_eigenstrain;
      82       15480 : }
      83             : 
      84             : void
      85        1720 : ComputeReducedOrderEigenstrain::sumEigenstrain()
      86             : {
      87        1720 :   _eigsum.resize(_qrule->n_points());
      88       17200 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      89             :   {
      90       15480 :     _eigsum[_qp].zero();
      91       30960 :     for (auto es : _eigenstrains)
      92       15480 :       _eigsum[_qp] += (*es)[_qp];
      93             :   }
      94        1720 : }
      95             : 
      96             : void
      97        1720 : ComputeReducedOrderEigenstrain::prepareEigenstrain()
      98             : {
      99             :   // The eigenstrains can either be constant in an element or linear in x, y, z
     100             :   // If constant, do volume averaging.
     101        1720 :   if (!_second_order || _qrule->n_points() == 1)
     102             :   {
     103             :     // Volume average
     104             :     _adjusted_eigenstrain.zero();
     105        1504 :     Real vol = 0.0;
     106       15040 :     for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
     107             :     {
     108       13536 :       _adjusted_eigenstrain += _eigsum[qp] * _JxW[qp] * _coord[qp];
     109       13536 :       vol += _JxW[qp] * _coord[qp];
     110             :     }
     111        1504 :     _adjusted_eigenstrain /= vol;
     112             :   }
     113             :   else
     114             :   {
     115             :     // 1 x y z
     116             : 
     117             :     // Six sets (one for each unique component of eigen) of qp data
     118         216 :     _A.resize(_qrule->n_points(), _ncols);
     119        1512 :     for (auto && b : _b)
     120        1296 :       b.resize(_qrule->n_points());
     121             : 
     122        2160 :     for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
     123             :     {
     124        1944 :       _A(qp, 0) = 1.0;
     125        5832 :       for (unsigned j = 1; j < _ncols; ++j)
     126        3888 :         _A(qp, j) = _q_point[qp](j - 1);
     127             : 
     128        1944 :       _b[0](qp) = _eigsum[qp](0, 0);
     129        1944 :       _b[1](qp) = _eigsum[qp](1, 1);
     130        1944 :       _b[2](qp) = _eigsum[qp](2, 2);
     131        1944 :       _b[3](qp) = _eigsum[qp](1, 2);
     132        1944 :       _b[4](qp) = _eigsum[qp](0, 2);
     133        1944 :       _b[5](qp) = _eigsum[qp](0, 1);
     134             :     }
     135             : 
     136         216 :     _A.get_transpose(_AT);
     137         216 :     _A.left_multiply(_AT);
     138        1512 :     for (unsigned i = 0; i < 6; ++i)
     139             :     {
     140        1296 :       _AT.vector_mult(_ATb, _b[i]);
     141             : 
     142        1296 :       _A.cholesky_solve(_ATb, _x[i]);
     143             :     }
     144             :   }
     145        1720 : }