www.mooseframework.org
Element.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 
10 #include "Element.h"
11 #include "Nonlinear3D.h"
12 #include "SolidModel.h"
13 
14 namespace SolidMechanics
15 {
16 
18  const std::string & /*name*/,
19  const InputParameters & /*parameters*/)
20  : Coupleable(&solid_model, false), _solid_model(solid_model)
21 {
22 }
23 
25 
27 
29 
30 Real
31 Element::detMatrix(const ColumnMajorMatrix & A)
32 {
33  mooseAssert(A.n() == 3 && A.m() == 3, "detMatrix requires 3x3 matrix");
34 
35  Real Axx = A(0, 0);
36  Real Axy = A(0, 1);
37  Real Axz = A(0, 2);
38  Real Ayx = A(1, 0);
39  Real Ayy = A(1, 1);
40  Real Ayz = A(1, 2);
41  Real Azx = A(2, 0);
42  Real Azy = A(2, 1);
43  Real Azz = A(2, 2);
44 
45  return Axx * Ayy * Azz + Axy * Ayz * Azx + Axz * Ayx * Azy - Azx * Ayy * Axz - Azy * Ayz * Axx -
46  Azz * Ayx * Axy;
47 }
48 
50 
51 void
52 Element::invertMatrix(const ColumnMajorMatrix & A, ColumnMajorMatrix & Ainv)
53 {
54  Real Axx = A(0, 0);
55  Real Axy = A(0, 1);
56  Real Axz = A(0, 2);
57  Real Ayx = A(1, 0);
58  Real Ayy = A(1, 1);
59  Real Ayz = A(1, 2);
60  Real Azx = A(2, 0);
61  Real Azy = A(2, 1);
62  Real Azz = A(2, 2);
63 
64  mooseAssert(detMatrix(A) > 0, "Matrix is not positive definite!");
65  Real detInv = 1 / detMatrix(A);
66 
67  Ainv(0, 0) = +(Ayy * Azz - Azy * Ayz) * detInv;
68  Ainv(0, 1) = -(Axy * Azz - Azy * Axz) * detInv;
69  Ainv(0, 2) = +(Axy * Ayz - Ayy * Axz) * detInv;
70  Ainv(1, 0) = -(Ayx * Azz - Azx * Ayz) * detInv;
71  Ainv(1, 1) = +(Axx * Azz - Azx * Axz) * detInv;
72  Ainv(1, 2) = -(Axx * Ayz - Ayx * Axz) * detInv;
73  Ainv(2, 0) = +(Ayx * Azy - Azx * Ayy) * detInv;
74  Ainv(2, 1) = -(Axx * Azy - Azx * Axy) * detInv;
75  Ainv(2, 2) = +(Axx * Ayy - Ayx * Axy) * detInv;
76 }
77 
79 
80 void
81 Element::rotateSymmetricTensor(const ColumnMajorMatrix & R,
82  const RealTensorValue & T,
83  RealTensorValue & result)
84 {
85 
86  // R T Rt
87  // 00 01 02 00 01 02 00 10 20
88  // 10 11 12 * 10 11 12 * 01 11 21
89  // 20 21 22 20 21 22 02 12 22
90  //
91  const Real T00 = R(0, 0) * T(0, 0) + R(0, 1) * T(1, 0) + R(0, 2) * T(2, 0);
92  const Real T01 = R(0, 0) * T(0, 1) + R(0, 1) * T(1, 1) + R(0, 2) * T(2, 1);
93  const Real T02 = R(0, 0) * T(0, 2) + R(0, 1) * T(1, 2) + R(0, 2) * T(2, 2);
94 
95  const Real T10 = R(1, 0) * T(0, 0) + R(1, 1) * T(1, 0) + R(1, 2) * T(2, 0);
96  const Real T11 = R(1, 0) * T(0, 1) + R(1, 1) * T(1, 1) + R(1, 2) * T(2, 1);
97  const Real T12 = R(1, 0) * T(0, 2) + R(1, 1) * T(1, 2) + R(1, 2) * T(2, 2);
98 
99  const Real T20 = R(2, 0) * T(0, 0) + R(2, 1) * T(1, 0) + R(2, 2) * T(2, 0);
100  const Real T21 = R(2, 0) * T(0, 1) + R(2, 1) * T(1, 1) + R(2, 2) * T(2, 1);
101  const Real T22 = R(2, 0) * T(0, 2) + R(2, 1) * T(1, 2) + R(2, 2) * T(2, 2);
102 
103  result = RealTensorValue(T00 * R(0, 0) + T01 * R(0, 1) + T02 * R(0, 2),
104  T00 * R(1, 0) + T01 * R(1, 1) + T02 * R(1, 2),
105  T00 * R(2, 0) + T01 * R(2, 1) + T02 * R(2, 2),
106 
107  T10 * R(0, 0) + T11 * R(0, 1) + T12 * R(0, 2),
108  T10 * R(1, 0) + T11 * R(1, 1) + T12 * R(1, 2),
109  T10 * R(2, 0) + T11 * R(2, 1) + T12 * R(2, 2),
110 
111  T20 * R(0, 0) + T21 * R(0, 1) + T22 * R(0, 2),
112  T20 * R(1, 0) + T21 * R(1, 1) + T22 * R(1, 2),
113  T20 * R(2, 0) + T21 * R(2, 1) + T22 * R(2, 2));
114 }
115 
117 
118 void
119 Element::rotateSymmetricTensor(const ColumnMajorMatrix & R,
120  const SymmTensor & T,
121  SymmTensor & result)
122 {
123 
124  // R T Rt
125  // 00 01 02 00 01 02 00 10 20
126  // 10 11 12 * 10 11 12 * 01 11 21
127  // 20 21 22 20 21 22 02 12 22
128  //
129  const Real T00 = R(0, 0) * T.xx() + R(0, 1) * T.xy() + R(0, 2) * T.zx();
130  const Real T01 = R(0, 0) * T.xy() + R(0, 1) * T.yy() + R(0, 2) * T.yz();
131  const Real T02 = R(0, 0) * T.zx() + R(0, 1) * T.yz() + R(0, 2) * T.zz();
132 
133  const Real T10 = R(1, 0) * T.xx() + R(1, 1) * T.xy() + R(1, 2) * T.zx();
134  const Real T11 = R(1, 0) * T.xy() + R(1, 1) * T.yy() + R(1, 2) * T.yz();
135  const Real T12 = R(1, 0) * T.zx() + R(1, 1) * T.yz() + R(1, 2) * T.zz();
136 
137  const Real T20 = R(2, 0) * T.xx() + R(2, 1) * T.xy() + R(2, 2) * T.zx();
138  const Real T21 = R(2, 0) * T.xy() + R(2, 1) * T.yy() + R(2, 2) * T.yz();
139  const Real T22 = R(2, 0) * T.zx() + R(2, 1) * T.yz() + R(2, 2) * T.zz();
140 
141  result.xx(T00 * R(0, 0) + T01 * R(0, 1) + T02 * R(0, 2));
142  result.yy(T10 * R(1, 0) + T11 * R(1, 1) + T12 * R(1, 2));
143  result.zz(T20 * R(2, 0) + T21 * R(2, 1) + T22 * R(2, 2));
144  result.xy(T00 * R(1, 0) + T01 * R(1, 1) + T02 * R(1, 2));
145  result.yz(T10 * R(2, 0) + T11 * R(2, 1) + T12 * R(2, 2));
146  result.zx(T00 * R(2, 0) + T01 * R(2, 1) + T02 * R(2, 2));
147 }
148 
150 
151 void
152 Element::unrotateSymmetricTensor(const ColumnMajorMatrix & R,
153  const SymmTensor & T,
154  SymmTensor & result)
155 {
156 
157  // Rt T R
158  // 00 10 20 00 01 02 00 01 02
159  // 01 11 21 * 10 11 12 * 10 11 12
160  // 02 12 22 20 21 22 20 21 22
161  //
162  const Real T00 = R(0, 0) * T.xx() + R(1, 0) * T.xy() + R(2, 0) * T.zx();
163  const Real T01 = R(0, 0) * T.xy() + R(1, 0) * T.yy() + R(2, 0) * T.yz();
164  const Real T02 = R(0, 0) * T.zx() + R(1, 0) * T.yz() + R(2, 0) * T.zz();
165 
166  const Real T10 = R(0, 1) * T.xx() + R(1, 1) * T.xy() + R(2, 1) * T.zx();
167  const Real T11 = R(0, 1) * T.xy() + R(1, 1) * T.yy() + R(2, 1) * T.yz();
168  const Real T12 = R(0, 1) * T.zx() + R(1, 1) * T.yz() + R(2, 1) * T.zz();
169 
170  const Real T20 = R(0, 2) * T.xx() + R(1, 2) * T.xy() + R(2, 2) * T.zx();
171  const Real T21 = R(0, 2) * T.xy() + R(1, 2) * T.yy() + R(2, 2) * T.yz();
172  const Real T22 = R(0, 2) * T.zx() + R(1, 2) * T.yz() + R(2, 2) * T.zz();
173 
174  result.xx(T00 * R(0, 0) + T01 * R(1, 0) + T02 * R(2, 0));
175  result.yy(T10 * R(0, 1) + T11 * R(1, 1) + T12 * R(2, 1));
176  result.zz(T20 * R(0, 2) + T21 * R(1, 2) + T22 * R(2, 2));
177  result.xy(T00 * R(0, 1) + T01 * R(1, 1) + T02 * R(2, 1));
178  result.yz(T10 * R(0, 2) + T11 * R(1, 2) + T12 * R(2, 2));
179  result.zx(T00 * R(0, 2) + T01 * R(1, 2) + T02 * R(2, 2));
180 }
181 
183 
184 void
185 Element::polarDecompositionEigen(const ColumnMajorMatrix & Fhat,
186  ColumnMajorMatrix & Rhat,
187  SymmTensor & strain_increment)
188 {
189  const int ND = 3;
190 
191  ColumnMajorMatrix eigen_value(ND, 1), eigen_vector(ND, ND);
192  ColumnMajorMatrix n1(ND, 1), n2(ND, 1), n3(ND, 1), N1(ND, 1), N2(ND, 1), N3(ND, 1);
193 
194  ColumnMajorMatrix Chat = Fhat.transpose() * Fhat;
195 
196  Chat.eigen(eigen_value, eigen_vector);
197 
198  for (int i = 0; i < ND; i++)
199  {
200  N1(i) = eigen_vector(i, 0);
201  N2(i) = eigen_vector(i, 1);
202  N3(i) = eigen_vector(i, 2);
203  }
204 
205  const Real lamda1 = std::sqrt(eigen_value(0));
206  const Real lamda2 = std::sqrt(eigen_value(1));
207  const Real lamda3 = std::sqrt(eigen_value(2));
208 
209  const Real log1 = std::log(lamda1);
210  const Real log2 = std::log(lamda2);
211  const Real log3 = std::log(lamda3);
212 
213  ColumnMajorMatrix Uhat =
214  N1 * N1.transpose() * lamda1 + N2 * N2.transpose() * lamda2 + N3 * N3.transpose() * lamda3;
215 
216  ColumnMajorMatrix invUhat(ND, ND);
217  invertMatrix(Uhat, invUhat);
218 
219  Rhat = Fhat * invUhat;
220 
221  strain_increment =
222  N1 * N1.transpose() * log1 + N2 * N2.transpose() * log2 + N3 * N3.transpose() * log3;
223 }
224 
226 
227 void
228 Element::fillMatrix(unsigned int qp,
229  const VariableGradient & grad_x,
230  const VariableGradient & grad_y,
231  const VariableGradient & grad_z,
232  ColumnMajorMatrix & A)
233 {
234  A(0, 0) = grad_x[qp](0);
235  A(0, 1) = grad_x[qp](1);
236  A(0, 2) = grad_x[qp](2);
237  A(1, 0) = grad_y[qp](0);
238  A(1, 1) = grad_y[qp](1);
239  A(1, 2) = grad_y[qp](2);
240  A(2, 0) = grad_z[qp](0);
241  A(2, 1) = grad_z[qp](1);
242  A(2, 2) = grad_z[qp](2);
243 }
244 
246 
247 } // namespace solid_mechanics
SymmTensor::xx
Real xx() const
Definition: SymmTensor.h:131
SymmTensor::zx
Real zx() const
Definition: SymmTensor.h:136
SolidMechanics::Element::unrotateSymmetricTensor
static void unrotateSymmetricTensor(const ColumnMajorMatrix &R, const SymmTensor &T, SymmTensor &result)
Definition: Element.C:152
SymmTensor::zz
Real zz() const
Definition: SymmTensor.h:133
SolidMechanics::Element::rotateSymmetricTensor
static void rotateSymmetricTensor(const ColumnMajorMatrix &R, const RealTensorValue &T, RealTensorValue &result)
Definition: Element.C:81
SolidMechanics::Element::polarDecompositionEigen
static void polarDecompositionEigen(const ColumnMajorMatrix &Fhat, ColumnMajorMatrix &Rhat, SymmTensor &strain_increment)
Definition: Element.C:185
SolidMechanics::Element::Element
Element(SolidModel &solid_model, const std::string &name, const InputParameters &parameters)
Definition: Element.C:17
Nonlinear3D.h
Element.h
SymmTensor::xy
Real xy() const
Definition: SymmTensor.h:134
SolidMechanics
Definition: AxisymmetricRZ.h:16
SolidMechanics::Element::fillMatrix
void fillMatrix(unsigned int qp, const VariableGradient &grad_x, const VariableGradient &grad_y, const VariableGradient &grad_z, ColumnMajorMatrix &A)
Definition: Element.C:228
SolidModel
SolidModel is the base class for all this module's solid mechanics material models.
Definition: SolidModel.h:33
SymmTensor::yz
Real yz() const
Definition: SymmTensor.h:135
SolidMechanics::Element::detMatrix
static Real detMatrix(const ColumnMajorMatrix &A)
Definition: Element.C:31
SymmTensor
Definition: SymmTensor.h:21
SymmTensor::yy
Real yy() const
Definition: SymmTensor.h:132
SolidMechanics::Element::~Element
virtual ~Element()
Definition: Element.C:26
SolidModel.h
SolidMechanics::Element::invertMatrix
static void invertMatrix(const ColumnMajorMatrix &A, ColumnMajorMatrix &Ainv)
Definition: Element.C:52