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 : #pragma once
11 :
12 : #include "SymmetricRankTwoTensor.h"
13 :
14 : // MOOSE includes
15 : #include "MooseEnum.h"
16 : #include "ColumnMajorMatrix.h"
17 : #include "MooseRandom.h"
18 : #include "SymmetricRankFourTensor.h"
19 : #include "Conversion.h"
20 : #include "MooseArray.h"
21 : #include "MathUtils.h"
22 :
23 : #include "libmesh/libmesh.h"
24 : #include "libmesh/vector_value.h"
25 : #include "libmesh/utility.h"
26 :
27 : // PETSc includes
28 : #include <petscblaslapack.h>
29 :
30 : // C++ includes/sqrt
31 : #include <iomanip>
32 : #include <ostream>
33 : #include <vector>
34 : #include <array>
35 :
36 : template <typename T>
37 : constexpr std::array<Real, SymmetricRankTwoTensorTempl<T>::N>
38 : SymmetricRankTwoTensorTempl<T>::identityCoords;
39 :
40 : namespace MathUtils
41 : {
42 : template <>
43 : void mooseSetToZero<SymmetricRankTwoTensor>(SymmetricRankTwoTensor & v);
44 : template <>
45 : void mooseSetToZero<ADSymmetricRankTwoTensor>(ADSymmetricRankTwoTensor & v);
46 : }
47 :
48 : template <typename T>
49 : MooseEnum
50 0 : SymmetricRankTwoTensorTempl<T>::fillMethodEnum()
51 : {
52 0 : return MooseEnum("autodetect=0 isotropic1=1 diagonal3=3 symmetric6=6", "autodetect");
53 : }
54 :
55 : template <typename T>
56 3211 : SymmetricRankTwoTensorTempl<T>::SymmetricRankTwoTensorTempl()
57 : {
58 3211 : std::fill(_vals.begin(), _vals.end(), 0.0);
59 3211 : }
60 :
61 : template <typename T>
62 42 : SymmetricRankTwoTensorTempl<T>::SymmetricRankTwoTensorTempl(const InitMethod init)
63 : {
64 42 : switch (init)
65 : {
66 40 : case initNone:
67 40 : break;
68 :
69 2 : case initIdentity:
70 2 : setToIdentity();
71 2 : break;
72 :
73 0 : default:
74 0 : mooseError("Unknown SymmetricRankTwoTensorTempl initialization pattern.");
75 : }
76 42 : }
77 :
78 : template <typename T>
79 452 : SymmetricRankTwoTensorTempl<T>::SymmetricRankTwoTensorTempl(
80 0 : const T & S11, const T & S22, const T & S33, const T & S23, const T & S13, const T & S12)
81 : {
82 452 : _vals[0] = S11;
83 452 : _vals[1] = S22;
84 452 : _vals[2] = S33;
85 452 : _vals[3] = mandelFactor(3) * S23;
86 452 : _vals[4] = mandelFactor(4) * S13;
87 452 : _vals[5] = mandelFactor(5) * S12;
88 452 : }
89 :
90 : template <typename T>
91 136 : SymmetricRankTwoTensorTempl<T>::operator RankTwoTensorTempl<T>()
92 : {
93 136 : return RankTwoTensorTempl<T>(_vals[0],
94 136 : _vals[5] / mandelFactor(5),
95 136 : _vals[4] / mandelFactor(4),
96 136 : _vals[5] / mandelFactor(5),
97 136 : _vals[1],
98 136 : _vals[3] / mandelFactor(3),
99 136 : _vals[4] / mandelFactor(4),
100 272 : _vals[3] / mandelFactor(3),
101 408 : _vals[2]);
102 : }
103 :
104 : template <typename T>
105 2 : SymmetricRankTwoTensorTempl<T>::SymmetricRankTwoTensorTempl(const T & S11,
106 : const T & S21,
107 : const T & S31,
108 : const T & S12,
109 : const T & S22,
110 : const T & S32,
111 : const T & S13,
112 : const T & S23,
113 0 : const T & S33)
114 : {
115 2 : _vals[0] = S11;
116 2 : _vals[1] = S22;
117 2 : _vals[2] = S33;
118 2 : _vals[3] = (S23 + S32) / mandelFactor(3);
119 2 : _vals[4] = (S13 + S31) / mandelFactor(4);
120 2 : _vals[5] = (S12 + S21) / mandelFactor(5);
121 2 : }
122 :
123 : template <typename T>
124 : SymmetricRankTwoTensorTempl<T>
125 6 : SymmetricRankTwoTensorTempl<T>::fromRawComponents(
126 : const T & S11, const T & S22, const T & S33, const T & S23, const T & S13, const T & S12)
127 : {
128 6 : SymmetricRankTwoTensorTempl<T> ret(SymmetricRankTwoTensorTempl<T>::initNone);
129 6 : ret._vals[0] = S11;
130 6 : ret._vals[1] = S22;
131 6 : ret._vals[2] = S33;
132 6 : ret._vals[3] = S23;
133 6 : ret._vals[4] = S13;
134 6 : ret._vals[5] = S12;
135 6 : return ret;
136 0 : }
137 :
138 : template <typename T>
139 0 : SymmetricRankTwoTensorTempl<T>::SymmetricRankTwoTensorTempl(const TensorValue<T> & a)
140 : {
141 0 : _vals[0] = a(0, 0);
142 0 : _vals[1] = a(1, 1);
143 0 : _vals[2] = a(2, 2);
144 0 : _vals[3] = (a(1, 2) + a(2, 1)) / mandelFactor(3);
145 0 : _vals[4] = (a(0, 2) + a(2, 0)) / mandelFactor(4);
146 0 : _vals[5] = (a(0, 1) + a(1, 0)) / mandelFactor(5);
147 0 : }
148 :
149 : template <typename T>
150 0 : SymmetricRankTwoTensorTempl<T>::SymmetricRankTwoTensorTempl(const TypeTensor<T> & a)
151 : {
152 0 : _vals[0] = a(0, 0);
153 0 : _vals[1] = a(1, 1);
154 0 : _vals[2] = a(2, 2);
155 0 : _vals[3] = (a(1, 2) + a(2, 1)) / mandelFactor(3);
156 0 : _vals[4] = (a(0, 2) + a(2, 0)) / mandelFactor(4);
157 0 : _vals[5] = (a(0, 1) + a(1, 0)) / mandelFactor(5);
158 0 : }
159 :
160 : /// named constructor for initializing symmetrically
161 : template <typename T>
162 : SymmetricRankTwoTensorTempl<T>
163 2 : SymmetricRankTwoTensorTempl<T>::initializeSymmetric(const TypeVector<T> & v0,
164 : const TypeVector<T> & v1,
165 : const TypeVector<T> & v2)
166 : {
167 : return SymmetricRankTwoTensorTempl<T>(
168 2 : v0(0), v1(0), v2(0), v0(1), v1(1), v2(1), v0(2), v1(2), v2(2));
169 : }
170 :
171 : template <typename T>
172 : void
173 92 : SymmetricRankTwoTensorTempl<T>::fillFromInputVector(const std::vector<T> & input,
174 : FillMethod fill_method)
175 : {
176 92 : if (fill_method != autodetect && fill_method != input.size())
177 0 : mooseError("Expected an input vector size of ",
178 : fill_method,
179 : " to fill the SymmetricRankTwoTensorTempl");
180 :
181 92 : switch (input.size())
182 : {
183 2 : case 1:
184 2 : _vals[0] = input[0];
185 2 : _vals[1] = input[0];
186 2 : _vals[2] = input[0];
187 2 : _vals[3] = 0.0;
188 2 : _vals[4] = 0.0;
189 2 : _vals[5] = 0.0;
190 2 : break;
191 :
192 2 : case 3:
193 2 : _vals[0] = input[0];
194 2 : _vals[1] = input[1];
195 2 : _vals[2] = input[2];
196 2 : _vals[3] = 0.0;
197 2 : _vals[4] = 0.0;
198 2 : _vals[5] = 0.0;
199 2 : break;
200 :
201 86 : case 6:
202 86 : _vals[0] = input[0];
203 86 : _vals[1] = input[1];
204 86 : _vals[2] = input[2];
205 86 : _vals[3] = mandelFactor(3) * input[3];
206 86 : _vals[4] = mandelFactor(4) * input[4];
207 86 : _vals[5] = mandelFactor(5) * input[5];
208 86 : break;
209 :
210 2 : default:
211 2 : mooseError("Please check the number of entries in the input vector for building "
212 : "a SymmetricRankTwoTensorTempl. It must be 1, 3, 6");
213 : }
214 90 : }
215 :
216 : template <typename T>
217 : void
218 8 : SymmetricRankTwoTensorTempl<T>::fillFromScalarVariable(const VariableValue & scalar_variable)
219 : {
220 8 : switch (scalar_variable.size())
221 : {
222 2 : case 1:
223 2 : _vals[0] = scalar_variable[0];
224 2 : _vals[1] = 0.0;
225 2 : _vals[2] = 0.0;
226 2 : _vals[3] = 0.0;
227 2 : _vals[4] = 0.0;
228 2 : _vals[5] = 0.0;
229 2 : break;
230 :
231 2 : case 3:
232 2 : _vals[0] = scalar_variable[0];
233 2 : _vals[1] = scalar_variable[1];
234 2 : _vals[2] = 0.0;
235 2 : _vals[3] = 0.0;
236 2 : _vals[4] = 0.0;
237 2 : _vals[5] = mandelFactor(5) * scalar_variable[2];
238 2 : break;
239 :
240 2 : case 6:
241 2 : _vals[0] = scalar_variable[0];
242 2 : _vals[1] = scalar_variable[1];
243 2 : _vals[2] = scalar_variable[2];
244 2 : _vals[3] = mandelFactor(3) * scalar_variable[3];
245 2 : _vals[4] = mandelFactor(4) * scalar_variable[4];
246 2 : _vals[5] = mandelFactor(5) * scalar_variable[5];
247 2 : break;
248 :
249 2 : default:
250 2 : mooseError("Only FIRST, THIRD, or SIXTH order scalar variable can be used to build "
251 : "a SymmetricRankTwoTensorTempl.");
252 : }
253 6 : }
254 :
255 : template <typename T>
256 : void
257 2 : SymmetricRankTwoTensorTempl<T>::rotate(const TypeTensor<T> & R)
258 : {
259 2 : auto M = SymmetricRankFourTensorTempl<T>::rotationMatrix(R);
260 2 : *this = M * (*this);
261 2 : }
262 :
263 : template <typename T>
264 : SymmetricRankTwoTensorTempl<T>
265 8 : SymmetricRankTwoTensorTempl<T>::transpose() const
266 : {
267 8 : return *this;
268 : }
269 :
270 : /// multiply vector v with row n of this tensor
271 : template <typename T>
272 : libMesh::VectorValue<T>
273 6 : SymmetricRankTwoTensorTempl<T>::row(const unsigned int n) const
274 : {
275 6 : switch (n)
276 : {
277 2 : case 0:
278 : return libMesh::VectorValue<T>(
279 2 : _vals[0], _vals[5] / mandelFactor(5), _vals[4] / mandelFactor(4));
280 2 : case 1:
281 : return libMesh::VectorValue<T>(
282 2 : _vals[5] / mandelFactor(5), _vals[1], _vals[3] / mandelFactor(3));
283 2 : case 2:
284 : return libMesh::VectorValue<T>(
285 2 : _vals[4] / mandelFactor(4), _vals[3] / mandelFactor(3), _vals[2]);
286 0 : default:
287 0 : mooseError("Invalid row");
288 : }
289 : }
290 :
291 : template <typename T>
292 : SymmetricRankTwoTensorTempl<T>
293 2 : SymmetricRankTwoTensorTempl<T>::timesTranspose(const RankTwoTensorTempl<T> & a)
294 : {
295 2 : return SymmetricRankTwoTensorTempl<T>(a(0, 0) * a(0, 0) + a(0, 1) * a(0, 1) + a(0, 2) * a(0, 2),
296 2 : a(1, 0) * a(1, 0) + a(1, 1) * a(1, 1) + a(1, 2) * a(1, 2),
297 2 : a(2, 0) * a(2, 0) + a(2, 1) * a(2, 1) + a(2, 2) * a(2, 2),
298 2 : a(1, 0) * a(2, 0) + a(1, 1) * a(2, 1) + a(1, 2) * a(2, 2),
299 2 : a(0, 0) * a(2, 0) + a(0, 1) * a(2, 1) + a(0, 2) * a(2, 2),
300 4 : a(0, 0) * a(1, 0) + a(0, 1) * a(1, 1) + a(0, 2) * a(1, 2));
301 : }
302 :
303 : template <typename T>
304 : SymmetricRankTwoTensorTempl<T>
305 6 : SymmetricRankTwoTensorTempl<T>::timesTranspose(const SymmetricRankTwoTensorTempl<T> & a)
306 : {
307 6 : return SymmetricRankTwoTensorTempl<T>::fromRawComponents(
308 6 : a._vals[0] * a._vals[0] + a._vals[4] * a._vals[4] / 2.0 + a._vals[5] * a._vals[5] / 2.0,
309 6 : a._vals[1] * a._vals[1] + a._vals[3] * a._vals[3] / 2.0 + a._vals[5] * a._vals[5] / 2.0,
310 6 : a._vals[2] * a._vals[2] + a._vals[3] * a._vals[3] / 2.0 + a._vals[4] * a._vals[4] / 2.0,
311 6 : a._vals[1] * a._vals[3] + a._vals[2] * a._vals[3] +
312 6 : a._vals[4] * a._vals[5] / MathUtils::sqrt2,
313 6 : a._vals[0] * a._vals[4] + a._vals[2] * a._vals[4] +
314 6 : a._vals[3] * a._vals[5] / MathUtils::sqrt2,
315 12 : a._vals[0] * a._vals[5] + a._vals[1] * a._vals[5] +
316 18 : a._vals[3] * a._vals[4] / MathUtils::sqrt2);
317 : }
318 :
319 : template <typename T>
320 : SymmetricRankTwoTensorTempl<T>
321 2 : SymmetricRankTwoTensorTempl<T>::transposeTimes(const RankTwoTensorTempl<T> & a)
322 : {
323 2 : return SymmetricRankTwoTensorTempl<T>(a(0, 0) * a(0, 0) + a(1, 0) * a(1, 0) + a(2, 0) * a(2, 0),
324 2 : a(0, 1) * a(0, 1) + a(1, 1) * a(1, 1) + a(2, 1) * a(2, 1),
325 2 : a(0, 2) * a(0, 2) + a(1, 2) * a(1, 2) + a(2, 2) * a(2, 2),
326 2 : a(0, 1) * a(0, 2) + a(1, 1) * a(1, 2) + a(2, 1) * a(2, 2),
327 2 : a(0, 0) * a(0, 2) + a(1, 0) * a(1, 2) + a(2, 0) * a(2, 2),
328 4 : a(0, 0) * a(0, 1) + a(1, 0) * a(1, 1) + a(2, 0) * a(2, 1));
329 : }
330 :
331 : template <typename T>
332 : SymmetricRankTwoTensorTempl<T>
333 2 : SymmetricRankTwoTensorTempl<T>::transposeTimes(const SymmetricRankTwoTensorTempl<T> & a)
334 : {
335 2 : return timesTranspose(a);
336 : }
337 :
338 : template <typename T>
339 : SymmetricRankTwoTensorTempl<T>
340 0 : SymmetricRankTwoTensorTempl<T>::plusTranspose(const RankTwoTensorTempl<T> & a)
341 : {
342 0 : return SymmetricRankTwoTensorTempl<T>(
343 0 : a(0, 0), a(0, 1), a(0, 2), a(1, 0), a(1, 1), a(1, 2), a(2, 0), a(2, 1), a(2, 2)) *
344 0 : 2.0;
345 : }
346 :
347 : template <typename T>
348 : SymmetricRankTwoTensorTempl<T>
349 20 : SymmetricRankTwoTensorTempl<T>::plusTranspose(const SymmetricRankTwoTensorTempl<T> & a)
350 : {
351 20 : return a * 2.0;
352 : }
353 :
354 : template <typename T>
355 : SymmetricRankTwoTensorTempl<T>
356 2 : SymmetricRankTwoTensorTempl<T>::square() const
357 : {
358 2 : return SymmetricRankTwoTensorTempl<T>::timesTranspose(*this);
359 : }
360 :
361 : template <typename T>
362 : void
363 2 : SymmetricRankTwoTensorTempl<T>::zero()
364 : {
365 14 : for (std::size_t i = 0; i < N; ++i)
366 12 : _vals[i] = 0.0;
367 2 : }
368 :
369 : template <typename T>
370 : SymmetricRankTwoTensorTempl<T> &
371 2 : SymmetricRankTwoTensorTempl<T>::operator+=(const SymmetricRankTwoTensorTempl<T> & a)
372 : {
373 14 : for (std::size_t i = 0; i < N; ++i)
374 12 : _vals[i] += a._vals[i];
375 2 : return *this;
376 : }
377 :
378 : template <typename T>
379 : SymmetricRankTwoTensorTempl<T> &
380 2 : SymmetricRankTwoTensorTempl<T>::operator-=(const SymmetricRankTwoTensorTempl<T> & a)
381 : {
382 14 : for (std::size_t i = 0; i < N; ++i)
383 12 : _vals[i] -= a._vals[i];
384 2 : return *this;
385 : }
386 :
387 : template <typename T>
388 : template <typename T2>
389 : SymmetricRankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
390 2 : SymmetricRankTwoTensorTempl<T>::operator+(const SymmetricRankTwoTensorTempl<T2> & a) const
391 : {
392 2 : SymmetricRankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype> result;
393 14 : for (std::size_t i = 0; i < N; ++i)
394 12 : result(i) = _vals[i] + a(i);
395 2 : return result;
396 0 : }
397 :
398 : template <typename T>
399 : template <typename T2>
400 : SymmetricRankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
401 34 : SymmetricRankTwoTensorTempl<T>::operator-(const SymmetricRankTwoTensorTempl<T2> & a) const
402 : {
403 34 : SymmetricRankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype> result(
404 : SymmetricRankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>::initNone);
405 238 : for (std::size_t i = 0; i < N; ++i)
406 204 : result(i) = _vals[i] - a(i);
407 34 : return result;
408 0 : }
409 :
410 : /// returns -_vals
411 : template <typename T>
412 : SymmetricRankTwoTensorTempl<T>
413 2 : SymmetricRankTwoTensorTempl<T>::operator-() const
414 : {
415 2 : return (*this) * -1.0;
416 : }
417 :
418 : template <typename T>
419 : SymmetricRankTwoTensorTempl<T> &
420 2 : SymmetricRankTwoTensorTempl<T>::operator*=(const T & a)
421 : {
422 14 : for (std::size_t i = 0; i < N; ++i)
423 12 : _vals[i] *= a;
424 2 : return *this;
425 : }
426 :
427 : template <typename T>
428 : SymmetricRankTwoTensorTempl<T> &
429 2 : SymmetricRankTwoTensorTempl<T>::operator/=(const T & a)
430 : {
431 14 : for (std::size_t i = 0; i < N; ++i)
432 12 : _vals[i] /= a;
433 2 : return *this;
434 : }
435 :
436 : template <typename T>
437 : T
438 2 : SymmetricRankTwoTensorTempl<T>::doubleContraction(const SymmetricRankTwoTensorTempl<T> & b) const
439 : {
440 2 : T sum = 0;
441 14 : for (unsigned int i = 0; i < N; ++i)
442 12 : sum += _vals[i] * b._vals[i];
443 2 : return sum;
444 0 : }
445 :
446 : template <typename T>
447 : SymmetricRankFourTensorTempl<T>
448 8 : SymmetricRankTwoTensorTempl<T>::outerProduct(const SymmetricRankTwoTensorTempl<T> & b) const
449 : {
450 8 : SymmetricRankFourTensorTempl<T> result;
451 8 : unsigned int index = 0;
452 56 : for (unsigned int i = 0; i < N; ++i)
453 : {
454 48 : const T & a = _vals[i];
455 336 : for (unsigned int j = 0; j < N; ++j)
456 288 : result._vals[index++] = a * b._vals[j];
457 : }
458 8 : return result;
459 0 : }
460 :
461 : template <typename T>
462 : SymmetricRankTwoTensorTempl<T>
463 22 : SymmetricRankTwoTensorTempl<T>::deviatoric() const
464 : {
465 22 : SymmetricRankTwoTensorTempl<T> deviatoric(*this);
466 22 : deviatoric.addIa(-1.0 / 3.0 * this->tr());
467 22 : return deviatoric;
468 0 : }
469 :
470 : template <typename T>
471 : T
472 4 : SymmetricRankTwoTensorTempl<T>::generalSecondInvariant() const
473 : {
474 4 : return _vals[0] * _vals[1] + _vals[0] * _vals[2] + _vals[1] * _vals[2] -
475 4 : (_vals[3] * _vals[3] + _vals[4] * _vals[4] + _vals[5] * _vals[5]) / 2.0;
476 : }
477 :
478 : template <typename T>
479 : T
480 18 : SymmetricRankTwoTensorTempl<T>::secondInvariant() const
481 : {
482 0 : T result;
483 18 : result = Utility::pow<2>(_vals[0] - _vals[1]) / 6.0;
484 18 : result += Utility::pow<2>(_vals[0] - _vals[2]) / 6.0;
485 18 : result += Utility::pow<2>(_vals[1] - _vals[2]) / 6.0;
486 18 : result += Utility::pow<2>(_vals[5]) / 2.0;
487 18 : result += Utility::pow<2>(_vals[4]) / 2.0;
488 18 : result += Utility::pow<2>(_vals[3]) / 2.0;
489 18 : return result;
490 0 : }
491 :
492 : template <typename T>
493 : SymmetricRankTwoTensorTempl<T>
494 4 : SymmetricRankTwoTensorTempl<T>::dsecondInvariant() const
495 : {
496 8 : return SymmetricRankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
497 : }
498 :
499 : template <typename T>
500 : SymmetricRankFourTensorTempl<T>
501 2 : SymmetricRankTwoTensorTempl<T>::d2secondInvariant() const
502 : {
503 2 : SymmetricRankFourTensorTempl<T> result;
504 :
505 14 : for (auto i : make_range(N))
506 84 : for (auto j : make_range(N))
507 72 : result(i, j) = 0.5 * (i == j) + 0.5 * (i == j) - (1.0 / 3.0) * (i < 3) * (j < 3);
508 :
509 2 : return result;
510 0 : }
511 :
512 : template <typename T>
513 : T
514 0 : SymmetricRankTwoTensorTempl<T>::trace() const
515 : {
516 0 : return tr();
517 : }
518 :
519 : template <typename T>
520 : SymmetricRankTwoTensorTempl<T>
521 4 : SymmetricRankTwoTensorTempl<T>::inverse() const
522 : {
523 4 : const auto d = det();
524 4 : if (d == 0.0)
525 2 : mooseException("Matrix not invertible");
526 2 : const SymmetricRankTwoTensorTempl<T> inv(
527 2 : _vals[2] * _vals[1] - _vals[3] * _vals[3] / 2.0,
528 2 : _vals[2] * _vals[0] - _vals[4] * _vals[4] / 2.0,
529 2 : _vals[0] * _vals[1] - _vals[5] * _vals[5] / 2.0,
530 2 : _vals[5] * _vals[4] / 2.0 - _vals[0] * _vals[3] / MathUtils::sqrt2,
531 2 : _vals[5] * _vals[3] / 2.0 - _vals[4] * _vals[1] / MathUtils::sqrt2,
532 2 : _vals[4] * _vals[3] / 2.0 - _vals[2] * _vals[5] / MathUtils::sqrt2);
533 2 : return inv * 1.0 / d;
534 0 : }
535 :
536 : template <typename T>
537 : SymmetricRankTwoTensorTempl<T>
538 4 : SymmetricRankTwoTensorTempl<T>::dtrace() const
539 : {
540 4 : return SymmetricRankTwoTensorTempl<T>(1, 1, 1, 0, 0, 0);
541 : }
542 :
543 : template <typename T>
544 : T
545 8 : SymmetricRankTwoTensorTempl<T>::thirdInvariant() const
546 : {
547 8 : const auto s = SymmetricRankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
548 8 : return s(0) * (s(1) * s(2) - s(3) / MathUtils::sqrt2 * s(3) / MathUtils::sqrt2) -
549 8 : s(5) / MathUtils::sqrt2 *
550 8 : (s(5) / MathUtils::sqrt2 * s(2) - s(3) / MathUtils::sqrt2 * s(4) / MathUtils::sqrt2) +
551 8 : s(4) / MathUtils::sqrt2 *
552 8 : (s(5) / MathUtils::sqrt2 * s(3) / MathUtils::sqrt2 - s(1) * s(4) / MathUtils::sqrt2);
553 0 : }
554 :
555 : template <typename T>
556 : SymmetricRankTwoTensorTempl<T>
557 4 : SymmetricRankTwoTensorTempl<T>::dthirdInvariant() const
558 : {
559 4 : const auto s = SymmetricRankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
560 4 : const T s3 = secondInvariant() / 3.0;
561 4 : return SymmetricRankTwoTensorTempl<T>(s(1) * s(2) - s(3) * s(3) / 2.0 + s3,
562 4 : s(0) * s(2) - s(4) * s(4) / 2.0 + s3,
563 4 : s(0) * s(1) - s(5) * s(5) / 2.0 + s3,
564 4 : s(5) * s(4) / 2.0 - s(0) * s(3) / MathUtils::sqrt2,
565 4 : s(5) * s(3) / 2.0 - s(1) * s(4) / MathUtils::sqrt2,
566 8 : s(4) * s(3) / 2.0 - s(5) * s(2) / MathUtils::sqrt2);
567 0 : }
568 :
569 : template <typename T>
570 : SymmetricRankFourTensorTempl<T>
571 2 : SymmetricRankTwoTensorTempl<T>::d2thirdInvariant() const
572 : {
573 2 : const auto s = SymmetricRankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
574 :
575 2 : SymmetricRankFourTensorTempl<T> d2;
576 14 : for (auto i : make_range(N))
577 84 : for (auto j : make_range(N))
578 144 : d2(i, j) = Real(i < 3) * s(j) / 3.0 / (j < 3 ? 1 : MathUtils::sqrt2) +
579 72 : Real(j < 3) * s(i) / 3.0 / (i < 3 ? 1 : MathUtils::sqrt2);
580 :
581 2 : d2(0, 1) += s(2);
582 2 : d2(0, 2) += s(1);
583 2 : d2(0, 3) -= s(3) / MathUtils::sqrt2;
584 :
585 2 : d2(1, 0) += s(2);
586 2 : d2(1, 2) += s(0);
587 2 : d2(1, 4) -= s(4) / MathUtils::sqrt2;
588 :
589 2 : d2(2, 0) += s(1);
590 2 : d2(2, 1) += s(0);
591 2 : d2(2, 5) -= s(5) / MathUtils::sqrt2;
592 :
593 2 : d2(3, 0) -= s(3) / MathUtils::sqrt2;
594 2 : d2(3, 3) -= s(0) / 2.0;
595 2 : d2(3, 4) += s(5) / 2.0 / MathUtils::sqrt2;
596 2 : d2(3, 5) += s(4) / 2.0 / MathUtils::sqrt2;
597 :
598 2 : d2(4, 1) -= s(4) / MathUtils::sqrt2;
599 2 : d2(4, 3) += s(5) / 2.0 / MathUtils::sqrt2;
600 2 : d2(4, 4) -= s(1) / 2.0;
601 2 : d2(4, 5) += s(3) / 2.0 / MathUtils::sqrt2;
602 :
603 2 : d2(5, 2) -= s(5) / MathUtils::sqrt2;
604 2 : d2(5, 3) += s(4) / 2.0 / MathUtils::sqrt2;
605 2 : d2(5, 4) += s(3) / 2.0 / MathUtils::sqrt2;
606 2 : d2(5, 5) -= s(2) / 2.0;
607 :
608 14 : for (auto i : make_range(N))
609 84 : for (auto j : make_range(N))
610 72 : d2(i, j) *= SymmetricRankFourTensorTempl<T>::mandelFactor(i, j);
611 :
612 4 : return d2;
613 0 : }
614 :
615 : template <typename T>
616 : T
617 12 : SymmetricRankTwoTensorTempl<T>::det() const
618 : {
619 12 : return _vals[0] * (_vals[2] * _vals[1] - _vals[3] * _vals[3] / 2.0) -
620 12 : _vals[5] / MathUtils::sqrt2 *
621 12 : (_vals[2] * _vals[5] / MathUtils::sqrt2 - _vals[3] * _vals[4] / 2.0) +
622 12 : _vals[4] / MathUtils::sqrt2 *
623 12 : (_vals[3] * _vals[5] / 2.0 - _vals[1] * _vals[4] / MathUtils::sqrt2);
624 : }
625 :
626 : template <typename T>
627 : SymmetricRankTwoTensorTempl<T>
628 2 : SymmetricRankTwoTensorTempl<T>::ddet() const
629 : {
630 : return SymmetricRankTwoTensorTempl<T>(
631 2 : _vals[2] * _vals[1] - _vals[3] * _vals[3] / 2.0,
632 2 : _vals[0] * _vals[2] - _vals[4] * _vals[4] / 2.0,
633 2 : _vals[0] * _vals[1] - _vals[5] * _vals[5] / 2.0,
634 2 : _vals[4] * _vals[5] / 2.0 - _vals[0] * _vals[3] / MathUtils::sqrt2,
635 2 : _vals[5] * _vals[3] / 2.0 - _vals[4] * _vals[1] / MathUtils::sqrt2,
636 4 : _vals[4] * _vals[3] / 2.0 - _vals[5] * _vals[2] / MathUtils::sqrt2);
637 : }
638 :
639 : template <typename T>
640 : void
641 2 : SymmetricRankTwoTensorTempl<T>::print(std::ostream & stm) const
642 : {
643 14 : for (std::size_t i = 0; i < N; ++i)
644 12 : stm << std::setw(15) << _vals[i] << std::endl;
645 2 : }
646 :
647 : template <typename T>
648 : void
649 2 : SymmetricRankTwoTensorTempl<T>::printReal(std::ostream & stm) const
650 : {
651 14 : for (std::size_t i = 0; i < N; ++i)
652 12 : stm << std::setw(15) << MetaPhysicL::raw_value(_vals[i]) << std::endl;
653 2 : }
654 :
655 : template <typename T>
656 : void
657 24 : SymmetricRankTwoTensorTempl<T>::addIa(const T & a)
658 : {
659 96 : for (unsigned int i = 0; i < 3; ++i)
660 72 : _vals[i] += a;
661 24 : }
662 :
663 : template <typename T>
664 : T
665 48 : SymmetricRankTwoTensorTempl<T>::L2norm() const
666 : {
667 48 : T norm = _vals[0] * _vals[0];
668 288 : for (unsigned int i = 1; i < N; ++i)
669 240 : norm += _vals[i] * _vals[i];
670 : using std::sqrt;
671 48 : return norm == 0.0 ? 0.0 : sqrt(norm);
672 0 : }
673 :
674 : template <typename T>
675 : void
676 0 : SymmetricRankTwoTensorTempl<T>::syev(const char *, std::vector<T> &, std::vector<T> &) const
677 : {
678 0 : mooseError("The syev method is only supported for Real valued tensors");
679 : }
680 :
681 : template <>
682 : void SymmetricRankTwoTensor::syev(const char * calculation_type,
683 : std::vector<Real> & eigvals,
684 : std::vector<Real> & a) const;
685 :
686 : template <typename T>
687 : void
688 0 : SymmetricRankTwoTensorTempl<T>::symmetricEigenvalues(std::vector<T> & eigvals) const
689 : {
690 0 : RankTwoTensorTempl<T> a;
691 0 : symmetricEigenvaluesEigenvectors(eigvals, a);
692 0 : }
693 :
694 : template <>
695 : void SymmetricRankTwoTensor::symmetricEigenvalues(std::vector<Real> & eigvals) const;
696 :
697 : template <typename T>
698 : void
699 2 : SymmetricRankTwoTensorTempl<T>::symmetricEigenvaluesEigenvectors(std::vector<T> &,
700 : RankTwoTensorTempl<T> &) const
701 : {
702 2 : mooseError(
703 : "symmetricEigenvaluesEigenvectors is only available for ordered tensor component types");
704 : }
705 :
706 : template <>
707 : void SymmetricRankTwoTensor::symmetricEigenvaluesEigenvectors(std::vector<Real> & eigvals,
708 : RankTwoTensor & eigvecs) const;
709 :
710 : template <typename T>
711 : void
712 2 : SymmetricRankTwoTensorTempl<T>::initRandom(unsigned int rand_seed)
713 : {
714 2 : MooseRandom::seed(rand_seed);
715 2 : }
716 :
717 : template <typename T>
718 : SymmetricRankTwoTensorTempl<T>
719 2 : SymmetricRankTwoTensorTempl<T>::genRandomSymmTensor(T scale, T offset)
720 : {
721 14 : auto r = [&]() { return (MooseRandom::rand() + offset) * scale; };
722 2 : return SymmetricRankTwoTensorTempl<T>(r(), r(), r(), r(), r(), r());
723 : }
724 :
725 : template <typename T>
726 : SymmetricRankTwoTensorTempl<T>
727 20 : SymmetricRankTwoTensorTempl<T>::selfOuterProduct(const TypeVector<T> & v)
728 : {
729 : return SymmetricRankTwoTensorTempl<T>(
730 20 : v(0) * v(0), v(1) * v(1), v(2) * v(2), v(1) * v(2), v(0) * v(2), v(0) * v(1));
731 : }
732 :
733 : template <typename T>
734 : SymmetricRankTwoTensorTempl<T>
735 2 : SymmetricRankTwoTensorTempl<T>::initialContraction(const SymmetricRankFourTensorTempl<T> & b) const
736 : {
737 2 : SymmetricRankTwoTensorTempl<T> result;
738 14 : for (auto i : make_range(N))
739 84 : for (auto j : make_range(N))
740 72 : result(j) += (*this)(i)*b(i, j);
741 2 : return result;
742 0 : }
743 :
744 : template <typename T>
745 : void
746 2 : SymmetricRankTwoTensorTempl<T>::setToIdentity()
747 : {
748 2 : std::copy(identityCoords.begin(), identityCoords.end(), _vals.begin());
749 2 : }
|
