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