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