www.mooseframework.org
ComputeReducedOrderEigenstrain.C
Go to the documentation of this file.
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 
11 
12 #include "MooseMesh.h"
13 #include "libmesh/quadrature.h"
14 
16 
19 {
21  params.addClassDescription("accepts eigenstrains and computes a reduced order eigenstrain for "
22  "consistency in the order of strain and eigenstrains.");
23  params.addRequiredParam<std::vector<MaterialPropertyName>>(
24  "input_eigenstrain_names", "List of eigenstrains to be applied in this strain calculation");
25  return params;
26 }
27 
29  : ComputeEigenstrainBase(parameters),
30  _input_eigenstrain_names(
31  getParam<std::vector<MaterialPropertyName>>("input_eigenstrain_names")),
32  _eigenstrains(_input_eigenstrain_names.size()),
33  _subproblem(*parameters.get<SubProblem *>("_subproblem")),
34  _ncols(1 + _subproblem.mesh().dimension()),
35  _second_order(_subproblem.mesh().hasSecondOrderElements()),
36  _eigsum(),
37  _A(),
38  _b(6),
39  _AT(),
40  _ATb(_ncols),
41  _x(6, DenseVector<Real>(_ncols)),
42  _vals(6),
43  _adjusted_eigenstrain()
44 {
45  for (unsigned int i = 0; i < _eigenstrains.size(); ++i)
46  {
48  _eigenstrains[i] = &getMaterialProperty<RankTwoTensor>(_input_eigenstrain_names[i]);
49  }
50 }
51 
52 void
54 {
55  _eigenstrain[_qp].zero();
56 }
57 
58 void
60 {
62 
64 
66 }
67 
68 void
70 {
71  if (_second_order)
72  {
73  for (unsigned i = 0; i < 6; ++i)
74  {
75  _vals[i] = _x[i](0);
76  for (unsigned j = 1; j < _ncols; ++j)
77  _vals[i] += _x[i](j) * _q_point[_qp](j - 1);
78  }
80  }
82 }
83 
84 void
86 {
87  _eigsum.resize(_qrule->n_points());
88  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
89  {
90  _eigsum[_qp].zero();
91  for (auto es : _eigenstrains)
92  _eigsum[_qp] += (*es)[_qp];
93  }
94 }
95 
96 void
98 {
99  // The eigenstrains can either be constant in an element or linear in x, y, z
100  // If constant, do volume averaging.
101  if (!_second_order || _qrule->n_points() == 1)
102  {
103  // Volume average
105  Real vol = 0.0;
106  for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
107  {
108  _adjusted_eigenstrain += _eigsum[qp] * _JxW[qp] * _coord[qp];
109  vol += _JxW[qp] * _coord[qp];
110  }
111  _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  _A.resize(_qrule->n_points(), _ncols);
119  for (auto && b : _b)
120  b.resize(_qrule->n_points());
121 
122  for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
123  {
124  _A(qp, 0) = 1.0;
125  for (unsigned j = 1; j < _ncols; ++j)
126  _A(qp, j) = _q_point[qp](j - 1);
127 
128  _b[0](qp) = _eigsum[qp](0, 0);
129  _b[1](qp) = _eigsum[qp](1, 1);
130  _b[2](qp) = _eigsum[qp](2, 2);
131  _b[3](qp) = _eigsum[qp](1, 2);
132  _b[4](qp) = _eigsum[qp](0, 2);
133  _b[5](qp) = _eigsum[qp](0, 1);
134  }
135 
138  for (unsigned i = 0; i < 6; ++i)
139  {
140  _AT.vector_mult(_ATb, _b[i]);
141 
142  _A.cholesky_solve(_ATb, _x[i]);
143  }
144  }
145 }
const MooseArray< Point > & _q_point
const QBase *const & _qrule
std::vector< RankTwoTensor > _eigsum
The sum of all eigenstrains at each integration point.
registerMooseObject("SolidMechanicsApp", ComputeReducedOrderEigenstrain)
void prepareEigenstrain()
Compute either the volume average or linear eigenstrain field in an element.
RankTwoTensor _adjusted_eigenstrain
Filled with _vals and subracted from strain.
virtual void computeProperties() override
MeshBase & mesh
const MooseArray< Real > & _JxW
const unsigned _ncols
Number of columns in A matrix (1 plus mesh dimension)
std::vector< DenseVector< Real > > _b
The b array holding the unique eigenstrain components for each integration point. ...
DenseMatrix< Real > _A
The (num q points x ncols) array for the least squares. Holds 1, xcoor, ycoor, zcoor.
void addRequiredParam(const std::string &name, const std::string &doc_string)
unsigned int _qp
void fillFromInputVector(const std::vector< Real > &input, FillMethod fill_method=autodetect)
static InputParameters validParams()
DenseVector< Real > _ATb
Transpose of A times b.
void get_transpose(DenseMatrix< Real > &dest) const
ComputeReducedOrderEigenstrain(const InputParameters &parameters)
ComputeEigenstrainBase is the base class for eigenstrain tensors.
const bool _second_order
Whether the mesh is made of second order elements.
std::vector< const MaterialProperty< RankTwoTensor > * > _eigenstrains
void sumEigenstrain()
Add contributions from every eigenstrain at each integration point.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
GenericMaterialProperty< RankTwoTensor, is_ad > & _eigenstrain
Stores the current total eigenstrain.
std::vector< DenseVector< Real > > _x
The solution vector for each unique component of the adjusted eigenstrain.
const std::string _base_name
Base name prepended to material property name.
void resize(const unsigned int new_m, const unsigned int new_n)
DenseMatrix< Real > _AT
Transpose of A.
void addClassDescription(const std::string &doc_string)
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
std::vector< Real > _vals
Vector to hold the adjusted strain as computed with _x.
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
const MooseArray< Real > & _coord
std::vector< MaterialPropertyName > _input_eigenstrain_names
void computeQpEigenstrain()
Compute the eigenstrain and store in _eigenstrain.
const Elem & get(const ElemType type_in)
virtual void left_multiply(const DenseMatrixBase< Real > &M2) override final
void vector_mult(DenseVector< Real > &dest, const DenseVector< Real > &arg) const