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         108 : ComputeReducedOrderEigenstrain::validParams()
      19             : {
      20         108 :   InputParameters params = ComputeEigenstrainBase::validParams();
      21         108 :   params.addClassDescription("accepts eigenstrains and computes a reduced order eigenstrain for "
      22             :                              "consistency in the order of strain and eigenstrains.");
      23         216 :   params.addRequiredParam<std::vector<MaterialPropertyName>>(
      24             :       "input_eigenstrain_names", "List of eigenstrains to be applied in this strain calculation");
      25         108 :   return params;
      26           0 : }
      27             : 
      28          81 : ComputeReducedOrderEigenstrain::ComputeReducedOrderEigenstrain(const InputParameters & parameters)
      29             :   : ComputeEigenstrainBase(parameters),
      30          81 :     _input_eigenstrain_names(
      31             :         getParam<std::vector<MaterialPropertyName>>("input_eigenstrain_names")),
      32          81 :     _eigenstrains(_input_eigenstrain_names.size()),
      33          81 :     _subproblem(*parameters.get<SubProblem *>("_subproblem")),
      34          81 :     _ncols(1 + _subproblem.mesh().dimension()),
      35          81 :     _second_order(_subproblem.mesh().hasSecondOrderElements()),
      36             :     _eigsum(),
      37          81 :     _A(),
      38          81 :     _b(6),
      39          81 :     _AT(),
      40          81 :     _ATb(_ncols),
      41          81 :     _x(6, DenseVector<Real>(_ncols)),
      42          81 :     _vals(6),
      43         162 :     _adjusted_eigenstrain()
      44             : {
      45         162 :   for (unsigned int i = 0; i < _eigenstrains.size(); ++i)
      46             :   {
      47          81 :     _input_eigenstrain_names[i] = _base_name + _input_eigenstrain_names[i];
      48          81 :     _eigenstrains[i] = &getMaterialProperty<RankTwoTensor>(_input_eigenstrain_names[i]);
      49             :   }
      50          81 : }
      51             : 
      52             : void
      53         684 : ComputeReducedOrderEigenstrain::initQpStatefulProperties()
      54             : {
      55         684 :   _eigenstrain[_qp].zero();
      56         684 : }
      57             : 
      58             : void
      59        2270 : ComputeReducedOrderEigenstrain::computeProperties()
      60             : {
      61        2270 :   sumEigenstrain();
      62             : 
      63        2270 :   prepareEigenstrain();
      64             : 
      65        2270 :   Material::computeProperties();
      66        2270 : }
      67             : 
      68             : void
      69       20430 : ComputeReducedOrderEigenstrain::computeQpEigenstrain()
      70             : {
      71       20430 :   if (_second_order)
      72             :   {
      73       20412 :     for (unsigned i = 0; i < 6; ++i)
      74             :     {
      75       17496 :       _vals[i] = _x[i](0);
      76       52488 :       for (unsigned j = 1; j < _ncols; ++j)
      77       34992 :         _vals[i] += _x[i](j) * _q_point[_qp](j - 1);
      78             :     }
      79        2916 :     _adjusted_eigenstrain.fillFromInputVector(_vals);
      80             :   }
      81       20430 :   _eigenstrain[_qp] = _adjusted_eigenstrain;
      82       20430 : }
      83             : 
      84             : void
      85        2270 : ComputeReducedOrderEigenstrain::sumEigenstrain()
      86             : {
      87        2270 :   _eigsum.resize(_qrule->n_points());
      88       22700 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      89             :   {
      90       20430 :     _eigsum[_qp].zero();
      91       40860 :     for (auto es : _eigenstrains)
      92       20430 :       _eigsum[_qp] += (*es)[_qp];
      93             :   }
      94        2270 : }
      95             : 
      96             : void
      97        2270 : 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        2270 :   if (!_second_order || _qrule->n_points() == 1)
     102             :   {
     103             :     // Volume average
     104             :     _adjusted_eigenstrain.zero();
     105        1946 :     Real vol = 0.0;
     106       19460 :     for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
     107             :     {
     108       17514 :       _adjusted_eigenstrain += _eigsum[qp] * _JxW[qp] * _coord[qp];
     109       17514 :       vol += _JxW[qp] * _coord[qp];
     110             :     }
     111        1946 :     _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         324 :     _A.resize(_qrule->n_points(), _ncols);
     119        2268 :     for (auto && b : _b)
     120        1944 :       b.resize(_qrule->n_points());
     121             : 
     122        3240 :     for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
     123             :     {
     124        2916 :       _A(qp, 0) = 1.0;
     125        8748 :       for (unsigned j = 1; j < _ncols; ++j)
     126        5832 :         _A(qp, j) = _q_point[qp](j - 1);
     127             : 
     128        2916 :       _b[0](qp) = _eigsum[qp](0, 0);
     129        2916 :       _b[1](qp) = _eigsum[qp](1, 1);
     130        2916 :       _b[2](qp) = _eigsum[qp](2, 2);
     131        2916 :       _b[3](qp) = _eigsum[qp](1, 2);
     132        2916 :       _b[4](qp) = _eigsum[qp](0, 2);
     133        2916 :       _b[5](qp) = _eigsum[qp](0, 1);
     134             :     }
     135             : 
     136         324 :     _A.get_transpose(_AT);
     137         324 :     _A.left_multiply(_AT);
     138        2268 :     for (unsigned i = 0; i < 6; ++i)
     139             :     {
     140        1944 :       _AT.vector_mult(_ATb, _b[i]);
     141             : 
     142        1944 :       _A.cholesky_solve(_ATb, _x[i]);
     143             :     }
     144             :   }
     145        2270 : }