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 "RankTwoTensor.h"
13 :
14 : // MOOSE includes
15 : #include "MooseEnum.h"
16 : #include "ColumnMajorMatrix.h"
17 : #include "MooseRandom.h"
18 : #include "RankFourTensor.h"
19 : #include "SymmetricRankTwoTensor.h"
20 : #include "Conversion.h"
21 : #include "MooseArray.h"
22 :
23 : #include "libmesh/libmesh.h"
24 : #include "libmesh/tensor_value.h"
25 : #include "libmesh/vector_value.h"
26 : #include "libmesh/utility.h"
27 :
28 : // PETSc includes
29 : #include <petscblaslapack.h>
30 :
31 : // C++ includes
32 : #include <iomanip>
33 : #include <ostream>
34 : #include <vector>
35 : #include <array>
36 :
37 : // Eigen includes
38 : #include <Eigen/Core>
39 : #include <Eigen/Eigenvalues>
40 :
41 : namespace Eigen
42 : {
43 : namespace internal
44 : {
45 : template <>
46 : struct cast_impl<ADReal, int>
47 : {
48 : static inline int run(const ADReal & x) { return static_cast<int>(MetaPhysicL::raw_value(x)); }
49 : };
50 : } // namespace internal
51 : } // namespace Eigen
52 :
53 : template <typename T>
54 : constexpr Real RankTwoTensorTempl<T>::identityCoords[];
55 :
56 : namespace MathUtils
57 : {
58 : template <>
59 : void mooseSetToZero<RankTwoTensor>(RankTwoTensor & v);
60 : template <>
61 : void mooseSetToZero<ADRankTwoTensor>(ADRankTwoTensor & v);
62 : }
63 :
64 : template <typename T>
65 : MooseEnum
66 0 : RankTwoTensorTempl<T>::fillMethodEnum()
67 : {
68 0 : return MooseEnum("autodetect=0 isotropic1=1 diagonal3=3 symmetric6=6 general=9", "autodetect");
69 : }
70 :
71 : template <typename T>
72 3389835 : RankTwoTensorTempl<T>::RankTwoTensorTempl()
73 : {
74 : mooseAssert(N == 3, "RankTwoTensorTempl is currently only tested for 3 dimensions.");
75 :
76 33898350 : for (unsigned int i = 0; i < N2; i++)
77 30508515 : _coords[i] = 0.0;
78 3389835 : }
79 :
80 : template <typename T>
81 1216 : RankTwoTensorTempl<T>::RankTwoTensorTempl(const InitMethod init)
82 : {
83 1216 : switch (init)
84 : {
85 1216 : case initNone:
86 1216 : break;
87 :
88 0 : case initIdentity:
89 0 : this->zero();
90 0 : for (const auto i : make_range(N))
91 0 : (*this)(i, i) = 1.0;
92 0 : break;
93 :
94 0 : default:
95 0 : mooseError("Unknown RankTwoTensorTempl initialization pattern.");
96 : }
97 1216 : }
98 :
99 : template <typename T>
100 0 : RankTwoTensorTempl<T>::RankTwoTensorTempl(const libMesh::TypeVector<T> & row1,
101 : const libMesh::TypeVector<T> & row2,
102 0 : const libMesh::TypeVector<T> & row3)
103 : {
104 0 : mooseDeprecated(
105 : "This constructor is deprecated in favor of RankTwoTensorTempl<T>::initializeFromRows");
106 :
107 : // Initialize the Tensor matrix from the passed in vectors
108 0 : for (const auto i : make_range(N))
109 0 : _coords[i] = row1(i);
110 :
111 0 : for (const auto i : make_range(N))
112 0 : _coords[N + i] = row2(i);
113 :
114 0 : for (const auto i : make_range(N))
115 0 : _coords[2 * N + i] = row3(i);
116 0 : }
117 :
118 : template <typename T>
119 : RankTwoTensorTempl<T>
120 4 : RankTwoTensorTempl<T>::initializeSymmetric(const libMesh::TypeVector<T> & v0,
121 : const libMesh::TypeVector<T> & v1,
122 : const libMesh::TypeVector<T> & v2)
123 : {
124 : return RankTwoTensorTempl<T>(v0(0),
125 4 : (v1(0) + v0(1)) / 2.0,
126 4 : (v2(0) + v0(2)) / 2.0,
127 4 : (v1(0) + v0(1)) / 2.0,
128 : v1(1),
129 4 : (v2(1) + v1(2)) / 2.0,
130 4 : (v2(0) + v0(2)) / 2.0,
131 8 : (v2(1) + v1(2)) / 2.0,
132 16 : v2(2));
133 : }
134 :
135 : template <typename T>
136 : RankTwoTensorTempl<T>
137 2 : RankTwoTensorTempl<T>::initializeFromRows(const libMesh::TypeVector<T> & row0,
138 : const libMesh::TypeVector<T> & row1,
139 : const libMesh::TypeVector<T> & row2)
140 : {
141 : return RankTwoTensorTempl<T>(
142 2 : row0(0), row1(0), row2(0), row0(1), row1(1), row2(1), row0(2), row1(2), row2(2));
143 : }
144 :
145 : template <typename T>
146 : RankTwoTensorTempl<T>
147 0 : RankTwoTensorTempl<T>::initializeFromColumns(const libMesh::TypeVector<T> & col0,
148 : const libMesh::TypeVector<T> & col1,
149 : const libMesh::TypeVector<T> & col2)
150 : {
151 : return RankTwoTensorTempl<T>(
152 0 : col0(0), col0(1), col0(2), col1(0), col1(1), col1(2), col2(0), col2(1), col2(2));
153 : }
154 :
155 : template <typename T>
156 28 : RankTwoTensorTempl<T>::RankTwoTensorTempl(
157 28 : const T & S11, const T & S22, const T & S33, const T & S23, const T & S13, const T & S12)
158 : {
159 28 : (*this)(0, 0) = S11;
160 28 : (*this)(1, 1) = S22;
161 28 : (*this)(2, 2) = S33;
162 28 : (*this)(1, 2) = (*this)(2, 1) = S23;
163 28 : (*this)(0, 2) = (*this)(2, 0) = S13;
164 28 : (*this)(0, 1) = (*this)(1, 0) = S12;
165 28 : }
166 :
167 : template <typename T>
168 746 : RankTwoTensorTempl<T>::RankTwoTensorTempl(const T & S11,
169 : const T & S21,
170 : const T & S31,
171 : const T & S12,
172 : const T & S22,
173 : const T & S32,
174 : const T & S13,
175 : const T & S23,
176 746 : const T & S33)
177 : {
178 746 : (*this)(0, 0) = S11;
179 746 : (*this)(1, 0) = S21;
180 746 : (*this)(2, 0) = S31;
181 746 : (*this)(0, 1) = S12;
182 746 : (*this)(1, 1) = S22;
183 746 : (*this)(2, 1) = S32;
184 746 : (*this)(0, 2) = S13;
185 746 : (*this)(1, 2) = S23;
186 746 : (*this)(2, 2) = S33;
187 746 : }
188 :
189 : template <typename T>
190 : void
191 248 : RankTwoTensorTempl<T>::fillFromInputVector(const std::vector<T> & input, FillMethod fill_method)
192 : {
193 248 : if (fill_method != autodetect && fill_method != input.size())
194 0 : mooseError("Expected an input vector size of ", fill_method, " to fill the RankTwoTensorTempl");
195 :
196 248 : switch (input.size())
197 : {
198 2 : case 1:
199 2 : this->zero();
200 2 : (*this)(0, 0) = input[0]; // S11
201 2 : (*this)(1, 1) = input[0]; // S22
202 2 : (*this)(2, 2) = input[0]; // S33
203 2 : break;
204 :
205 14 : case 3:
206 14 : this->zero();
207 14 : (*this)(0, 0) = input[0]; // S11
208 14 : (*this)(1, 1) = input[1]; // S22
209 14 : (*this)(2, 2) = input[2]; // S33
210 14 : break;
211 :
212 44 : case 6:
213 44 : (*this)(0, 0) = input[0]; // S11
214 44 : (*this)(1, 1) = input[1]; // S22
215 44 : (*this)(2, 2) = input[2]; // S33
216 44 : (*this)(1, 2) = (*this)(2, 1) = input[3]; // S23
217 44 : (*this)(0, 2) = (*this)(2, 0) = input[4]; // S13
218 44 : (*this)(0, 1) = (*this)(1, 0) = input[5]; // S12
219 44 : break;
220 :
221 188 : case 9:
222 188 : (*this)(0, 0) = input[0]; // S11
223 188 : (*this)(1, 0) = input[1]; // S21
224 188 : (*this)(2, 0) = input[2]; // S31
225 188 : (*this)(0, 1) = input[3]; // S12
226 188 : (*this)(1, 1) = input[4]; // S22
227 188 : (*this)(2, 1) = input[5]; // S32
228 188 : (*this)(0, 2) = input[6]; // S13
229 188 : (*this)(1, 2) = input[7]; // S23
230 188 : (*this)(2, 2) = input[8]; // S33
231 188 : break;
232 :
233 0 : default:
234 0 : mooseError("Please check the number of entries in the input vector for building "
235 : "a RankTwoTensorTempl. It must be 1, 3, 6, or 9");
236 : }
237 248 : }
238 :
239 : template <typename T>
240 : void
241 0 : RankTwoTensorTempl<T>::fillFromScalarVariable(const VariableValue & scalar_variable)
242 : {
243 0 : switch (scalar_variable.size())
244 : {
245 0 : case 1:
246 0 : this->zero();
247 0 : (*this)(0, 0) = scalar_variable[0]; // S11
248 0 : break;
249 :
250 0 : case 3:
251 0 : this->zero();
252 0 : (*this)(0, 0) = scalar_variable[0]; // S11
253 0 : (*this)(1, 1) = scalar_variable[1]; // S22
254 0 : (*this)(0, 1) = (*this)(1, 0) = scalar_variable[2]; // S12
255 0 : break;
256 :
257 0 : case 6:
258 0 : (*this)(0, 0) = scalar_variable[0]; // S11
259 0 : (*this)(1, 1) = scalar_variable[1]; // S22
260 0 : (*this)(2, 2) = scalar_variable[2]; // S33
261 0 : (*this)(1, 2) = (*this)(2, 1) = scalar_variable[3]; // S23
262 0 : (*this)(0, 2) = (*this)(2, 0) = scalar_variable[4]; // S13
263 0 : (*this)(0, 1) = (*this)(1, 0) = scalar_variable[5]; // S12
264 0 : break;
265 :
266 0 : default:
267 0 : mooseError("Only FIRST, THIRD, or SIXTH order scalar variable can be used to build "
268 : "a RankTwoTensorTempl.");
269 : }
270 0 : }
271 :
272 : template <typename T>
273 : libMesh::VectorValue<T>
274 1230 : RankTwoTensorTempl<T>::column(const unsigned int c) const
275 : {
276 1230 : return libMesh::VectorValue<T>((*this)(0, c), (*this)(1, c), (*this)(2, c));
277 : }
278 :
279 : template <typename T>
280 : RankTwoTensorTempl<T>
281 4 : RankTwoTensorTempl<T>::timesTranspose(const RankTwoTensorTempl<T> & a)
282 : {
283 4 : return a * a.transpose();
284 : }
285 :
286 : template <typename T>
287 : RankTwoTensorTempl<T>
288 2 : RankTwoTensorTempl<T>::transposeTimes(const RankTwoTensorTempl<T> & a)
289 : {
290 2 : return a.transpose() * a;
291 : }
292 :
293 : template <typename T>
294 : RankTwoTensorTempl<T>
295 2640 : RankTwoTensorTempl<T>::plusTranspose(const RankTwoTensorTempl<T> & a)
296 : {
297 2640 : return a + a.transpose();
298 : }
299 :
300 : template <typename T>
301 : RankTwoTensorTempl<T>
302 4 : RankTwoTensorTempl<T>::square() const
303 : {
304 4 : return *this * *this;
305 : }
306 :
307 : template <typename T>
308 : RankTwoTensorTempl<T>
309 18 : RankTwoTensorTempl<T>::rotated(const RankTwoTensorTempl<T> & R) const
310 : {
311 18 : RankTwoTensorTempl<T> result(*this);
312 18 : result.rotate(R);
313 18 : return result;
314 0 : }
315 :
316 : template <typename T>
317 : void
318 34 : RankTwoTensorTempl<T>::rotate(const RankTwoTensorTempl<T> & R)
319 : {
320 34 : RankTwoTensorTempl<T> temp;
321 136 : for (const auto i : make_range(N))
322 : {
323 102 : const auto i1 = i * N;
324 408 : for (const auto j : make_range(N))
325 : {
326 : // tmp += R(i,k)*R(j,l)*(*this)(k,l);
327 : // clang-format off
328 306 : const auto j1 = j * N;
329 306 : T tmp = R._coords[i1 + 0] * R._coords[j1 + 0] * (*this)(0, 0) +
330 306 : R._coords[i1 + 0] * R._coords[j1 + 1] * (*this)(0, 1) +
331 306 : R._coords[i1 + 0] * R._coords[j1 + 2] * (*this)(0, 2) +
332 306 : R._coords[i1 + 1] * R._coords[j1 + 0] * (*this)(1, 0) +
333 306 : R._coords[i1 + 1] * R._coords[j1 + 1] * (*this)(1, 1) +
334 306 : R._coords[i1 + 1] * R._coords[j1 + 2] * (*this)(1, 2) +
335 306 : R._coords[i1 + 2] * R._coords[j1 + 0] * (*this)(2, 0) +
336 306 : R._coords[i1 + 2] * R._coords[j1 + 1] * (*this)(2, 1) +
337 306 : R._coords[i1 + 2] * R._coords[j1 + 2] * (*this)(2, 2);
338 : // clang-format on
339 306 : temp._coords[i1 + j] = tmp;
340 : }
341 : }
342 340 : for (unsigned int i = 0; i < N2; i++)
343 306 : _coords[i] = temp._coords[i];
344 34 : }
345 :
346 : template <typename T>
347 : RankTwoTensorTempl<T>
348 0 : RankTwoTensorTempl<T>::rotateXyPlane(T a)
349 : {
350 : using std::cos, std::sin;
351 0 : T c = cos(a);
352 0 : T s = sin(a);
353 0 : T x = (*this)(0, 0) * c * c + (*this)(1, 1) * s * s + 2.0 * (*this)(0, 1) * c * s;
354 0 : T y = (*this)(0, 0) * s * s + (*this)(1, 1) * c * c - 2.0 * (*this)(0, 1) * c * s;
355 0 : T xy = ((*this)(1, 1) - (*this)(0, 0)) * c * s + (*this)(0, 1) * (c * c - s * s);
356 :
357 0 : RankTwoTensorTempl<T> b(*this);
358 :
359 0 : b(0, 0) = x;
360 0 : b(1, 1) = y;
361 0 : b(1, 0) = b(0, 1) = xy;
362 :
363 0 : return b;
364 0 : }
365 :
366 : template <typename T>
367 : RankTwoTensorTempl<T>
368 2812 : RankTwoTensorTempl<T>::transpose() const
369 : {
370 2812 : return libMesh::TensorValue<T>::transpose();
371 : }
372 :
373 : template <typename T>
374 : RankTwoTensorTempl<T> &
375 2 : RankTwoTensorTempl<T>::operator+=(const RankTwoTensorTempl<T> & a)
376 : {
377 2 : libMesh::TensorValue<T>::operator+=(a);
378 2 : return *this;
379 : }
380 :
381 : template <typename T>
382 : RankTwoTensorTempl<T> &
383 2 : RankTwoTensorTempl<T>::operator-=(const RankTwoTensorTempl<T> & a)
384 : {
385 2 : libMesh::TensorValue<T>::operator-=(a);
386 2 : return *this;
387 : }
388 :
389 : template <typename T>
390 : RankTwoTensorTempl<T>
391 2 : RankTwoTensorTempl<T>::operator-() const
392 : {
393 2 : return libMesh::TensorValue<T>::operator-();
394 : }
395 :
396 : template <typename T>
397 : RankTwoTensorTempl<T> &
398 4 : RankTwoTensorTempl<T>::operator*=(const T & a)
399 : {
400 4 : libMesh::TensorValue<T>::operator*=(a);
401 4 : return *this;
402 : }
403 :
404 : template <typename T>
405 : RankTwoTensorTempl<T> &
406 4 : RankTwoTensorTempl<T>::operator/=(const T & a)
407 : {
408 4 : libMesh::TensorValue<T>::operator/=(a);
409 4 : return *this;
410 : }
411 :
412 : template <typename T>
413 : RankTwoTensorTempl<T> &
414 0 : RankTwoTensorTempl<T>::operator*=(const libMesh::TypeTensor<T> & a)
415 : {
416 0 : *this = *this * a;
417 0 : return *this;
418 : }
419 :
420 : template <typename T>
421 : bool
422 6 : RankTwoTensorTempl<T>::operator==(const RankTwoTensorTempl<T> & a) const
423 : {
424 18 : for (const auto i : make_range(N))
425 52 : for (const auto j : make_range(N))
426 40 : if (!MooseUtils::absoluteFuzzyEqual((*this)(i, j), a(i, j)))
427 2 : return false;
428 :
429 4 : return true;
430 : }
431 :
432 : template <typename T>
433 : bool
434 46 : RankTwoTensorTempl<T>::isSymmetric() const
435 : {
436 46 : auto test = MetaPhysicL::raw_value(*this - transpose());
437 92 : return MooseUtils::absoluteFuzzyEqual(test.norm_sq(), 0);
438 46 : }
439 :
440 : template <typename T>
441 : RankTwoTensorTempl<T> &
442 0 : RankTwoTensorTempl<T>::operator=(const ColumnMajorMatrixTempl<T> & a)
443 : {
444 0 : if (a.n() != N || a.m() != N)
445 0 : mooseError("Dimensions of ColumnMajorMatrixTempl<T> are incompatible with RankTwoTensorTempl");
446 :
447 0 : const T * cmm_rawdata = a.rawData();
448 0 : for (const auto i : make_range(N))
449 0 : for (const auto j : make_range(N))
450 0 : _coords[i * N + j] = cmm_rawdata[i + j * N];
451 :
452 0 : return *this;
453 : }
454 :
455 : template <typename T>
456 : T
457 4 : RankTwoTensorTempl<T>::doubleContraction(const RankTwoTensorTempl<T> & b) const
458 : {
459 : // deprecate this!
460 4 : return libMesh::TensorValue<T>::contract(b);
461 : }
462 :
463 : template <typename T>
464 : RankThreeTensorTempl<T>
465 2 : RankTwoTensorTempl<T>::contraction(const RankThreeTensorTempl<T> & b) const
466 : {
467 2 : RankThreeTensorTempl<T> result;
468 :
469 8 : for (const auto i : make_range(N))
470 24 : for (const auto j : make_range(N))
471 72 : for (const auto k : make_range(N))
472 216 : for (const auto l : make_range(N))
473 162 : result(i, k, l) += (*this)(i, j) * b(j, k, l);
474 :
475 2 : return result;
476 0 : }
477 :
478 : template <typename T>
479 : RankThreeTensorTempl<T>
480 2 : RankTwoTensorTempl<T>::mixedProductJkI(const libMesh::VectorValue<T> & b) const
481 : {
482 2 : RankThreeTensorTempl<T> result;
483 :
484 8 : for (const auto i : make_range(N))
485 24 : for (const auto j : make_range(N))
486 72 : for (const auto k : make_range(N))
487 54 : result(i, j, k) += (*this)(j, k) * b(i);
488 :
489 2 : return result;
490 0 : }
491 :
492 : template <typename T>
493 : RankTwoTensorTempl<T>
494 2650 : RankTwoTensorTempl<T>::deviatoric() const
495 : {
496 2650 : RankTwoTensorTempl<T> deviatoric(*this);
497 : // actually construct deviatoric part
498 2650 : deviatoric.addIa(-1.0 / 3.0 * this->tr());
499 2650 : return deviatoric;
500 0 : }
501 :
502 : template <typename T>
503 : T
504 4 : RankTwoTensorTempl<T>::generalSecondInvariant() const
505 : {
506 4 : return (*this)(0, 0) * (*this)(1, 1) + (*this)(0, 0) * (*this)(2, 2) +
507 4 : (*this)(1, 1) * (*this)(2, 2) - (*this)(0, 1) * (*this)(1, 0) -
508 4 : (*this)(0, 2) * (*this)(2, 0) - (*this)(1, 2) * (*this)(2, 1);
509 : }
510 :
511 : template <typename T>
512 : T
513 1572 : RankTwoTensorTempl<T>::secondInvariant() const
514 : {
515 0 : T result;
516 :
517 : // RankTwoTensorTempl<T> deviatoric(*this);
518 : // deviatoric.addIa(-1.0/3.0 * this->tr()); // actually construct deviatoric part
519 : // result = 0.5*(deviatoric + deviatoric.transpose()).doubleContraction(deviatoric +
520 : // deviatoric.transpose());
521 1572 : result = Utility::pow<2>((*this)(0, 0) - (*this)(1, 1)) / 6.0;
522 1572 : result += Utility::pow<2>((*this)(0, 0) - (*this)(2, 2)) / 6.0;
523 1572 : result += Utility::pow<2>((*this)(1, 1) - (*this)(2, 2)) / 6.0;
524 1572 : result += Utility::pow<2>((*this)(0, 1) + (*this)(1, 0)) / 4.0;
525 1572 : result += Utility::pow<2>((*this)(0, 2) + (*this)(2, 0)) / 4.0;
526 1572 : result += Utility::pow<2>((*this)(1, 2) + (*this)(2, 1)) / 4.0;
527 1572 : return result;
528 0 : }
529 :
530 : template <typename T>
531 : RankTwoTensorTempl<T>
532 1000 : RankTwoTensorTempl<T>::dsecondInvariant() const
533 : {
534 1000 : return RankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
535 : }
536 :
537 : template <typename T>
538 : RankFourTensorTempl<T>
539 10 : RankTwoTensorTempl<T>::d2secondInvariant() const
540 : {
541 10 : RankFourTensorTempl<T> result;
542 40 : for (const auto i : make_range(N))
543 120 : for (const auto j : make_range(N))
544 360 : for (const auto k : make_range(N))
545 1080 : for (const auto l : make_range(N))
546 1620 : result(i, j, k, l) = 0.5 * (i == k) * (j == l) + 0.5 * (i == l) * (j == k) -
547 810 : (1.0 / 3.0) * (i == j) * (k == l);
548 10 : return result;
549 0 : }
550 :
551 : template <typename T>
552 : T
553 10080002 : RankTwoTensorTempl<T>::trace() const
554 : {
555 : // deprecate this!
556 10080002 : return this->tr();
557 : }
558 :
559 : template <typename T>
560 : RankTwoTensorTempl<T>
561 8 : RankTwoTensorTempl<T>::dtrace() const
562 : {
563 8 : return RankTwoTensorTempl<T>(1, 0, 0, 0, 1, 0, 0, 0, 1);
564 : }
565 :
566 : template <typename T>
567 : RankTwoTensorTempl<T>
568 12 : RankTwoTensorTempl<T>::inverse() const
569 : {
570 12 : return libMesh::TensorValue<T>::inverse();
571 : }
572 :
573 : template <typename T>
574 : T
575 630 : RankTwoTensorTempl<T>::thirdInvariant() const
576 : {
577 630 : const auto s = RankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
578 630 : return s(0, 0) * (s(1, 1) * s(2, 2) - s(2, 1) * s(1, 2)) -
579 630 : s(1, 0) * (s(0, 1) * s(2, 2) - s(2, 1) * s(0, 2)) +
580 1260 : s(2, 0) * (s(0, 1) * s(1, 2) - s(1, 1) * s(0, 2));
581 630 : }
582 :
583 : template <typename T>
584 : RankTwoTensorTempl<T>
585 996 : RankTwoTensorTempl<T>::dthirdInvariant() const
586 : {
587 996 : const auto s = RankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
588 996 : const T s3 = secondInvariant() / 3.0;
589 :
590 996 : RankTwoTensorTempl<T> d;
591 996 : d(0, 0) = s(1, 1) * s(2, 2) - s(2, 1) * s(1, 2) + s3;
592 996 : d(0, 1) = s(2, 0) * s(1, 2) - s(1, 0) * s(2, 2);
593 996 : d(0, 2) = s(1, 0) * s(2, 1) - s(2, 0) * s(1, 1);
594 996 : d(1, 0) = s(2, 1) * s(0, 2) - s(0, 1) * s(2, 2);
595 996 : d(1, 1) = s(0, 0) * s(2, 2) - s(2, 0) * s(0, 2) + s3;
596 996 : d(1, 2) = s(2, 0) * s(0, 1) - s(0, 0) * s(2, 1);
597 996 : d(2, 0) = s(0, 1) * s(1, 2) - s(1, 1) * s(0, 2);
598 996 : d(2, 1) = s(1, 0) * s(0, 2) - s(0, 0) * s(1, 2);
599 996 : d(2, 2) = s(0, 0) * s(1, 1) - s(1, 0) * s(0, 1) + s3;
600 1992 : return d;
601 996 : }
602 :
603 : template <typename T>
604 : RankFourTensorTempl<T>
605 10 : RankTwoTensorTempl<T>::d2thirdInvariant() const
606 : {
607 10 : const auto s = RankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
608 :
609 10 : RankFourTensorTempl<T> d2;
610 40 : for (const auto i : make_range(N))
611 120 : for (const auto j : make_range(N))
612 360 : for (const auto k : make_range(N))
613 1080 : for (const auto l : make_range(N))
614 : {
615 810 : d2(i, j, k, l) = Real(i == j) * s(k, l) / 3.0 + Real(k == l) * s(i, j) / 3.0;
616 : // for (const auto a: make_range(N))
617 : // for (const auto b: make_range(N))
618 : // d2(i, j, k, l) += 0.5*(PermutationTensor::eps(i, k, a)*PermutationTensor::eps(j, l,
619 : // b) + PermutationTensor::eps(i, l, a)*PermutationTensor::eps(j, k, b))*s(a, b);
620 : }
621 :
622 : // I'm not sure which is more readable: the above
623 : // PermutationTensor stuff, or the stuff below.
624 : // Anyway, they yield the same result, and so i leave
625 : // both of them here to enlighten you!
626 :
627 10 : d2(0, 0, 1, 1) += s(2, 2);
628 10 : d2(0, 0, 1, 2) -= s(2, 1);
629 10 : d2(0, 0, 2, 1) -= s(1, 2);
630 10 : d2(0, 0, 2, 2) += s(1, 1);
631 :
632 10 : d2(0, 1, 0, 1) -= s(2, 2) / 2.0;
633 10 : d2(0, 1, 1, 0) -= s(2, 2) / 2.0;
634 10 : d2(0, 1, 0, 2) += s(1, 2) / 2.0;
635 10 : d2(0, 1, 2, 0) += s(1, 2) / 2.0;
636 10 : d2(0, 1, 1, 2) += s(2, 0) / 2.0;
637 10 : d2(0, 1, 2, 1) += s(2, 0) / 2.0;
638 10 : d2(0, 1, 2, 2) -= s(1, 0);
639 :
640 10 : d2(0, 2, 0, 1) += s(2, 1) / 2.0;
641 10 : d2(0, 2, 1, 0) += s(2, 1) / 2.0;
642 10 : d2(0, 2, 0, 2) -= s(1, 1) / 2.0;
643 10 : d2(0, 2, 2, 0) -= s(1, 1) / 2.0;
644 10 : d2(0, 2, 1, 1) -= s(2, 0);
645 10 : d2(0, 2, 1, 2) += s(1, 0) / 2.0;
646 10 : d2(0, 2, 2, 1) += s(1, 0) / 2.0;
647 :
648 10 : d2(1, 0, 0, 1) -= s(2, 2) / 2.0;
649 10 : d2(1, 0, 1, 0) -= s(2, 2) / 2.0;
650 10 : d2(1, 0, 0, 2) += s(1, 2) / 2.0;
651 10 : d2(1, 0, 2, 0) += s(1, 2) / 2.0;
652 10 : d2(1, 0, 1, 2) += s(2, 0) / 2.0;
653 10 : d2(1, 0, 2, 1) += s(2, 0) / 2.0;
654 10 : d2(1, 0, 2, 2) -= s(1, 0);
655 :
656 10 : d2(1, 1, 0, 0) += s(2, 2);
657 10 : d2(1, 1, 0, 2) -= s(2, 0);
658 10 : d2(1, 1, 2, 0) -= s(2, 0);
659 10 : d2(1, 1, 2, 2) += s(0, 0);
660 :
661 10 : d2(1, 2, 0, 0) -= s(2, 1);
662 10 : d2(1, 2, 0, 1) += s(2, 0) / 2.0;
663 10 : d2(1, 2, 1, 0) += s(2, 0) / 2.0;
664 10 : d2(1, 2, 0, 2) += s(0, 1) / 2.0;
665 10 : d2(1, 2, 2, 0) += s(0, 1) / 2.0;
666 10 : d2(1, 2, 1, 2) -= s(0, 0) / 2.0;
667 10 : d2(1, 2, 2, 1) -= s(0, 0) / 2.0;
668 :
669 10 : d2(2, 0, 0, 1) += s(2, 1) / 2.0;
670 10 : d2(2, 0, 1, 0) += s(2, 1) / 2.0;
671 10 : d2(2, 0, 0, 2) -= s(1, 1) / 2.0;
672 10 : d2(2, 0, 2, 0) -= s(1, 1) / 2.0;
673 10 : d2(2, 0, 1, 1) -= s(2, 0);
674 10 : d2(2, 0, 1, 2) += s(1, 0) / 2.0;
675 10 : d2(2, 0, 2, 1) += s(1, 0) / 2.0;
676 :
677 10 : d2(2, 1, 0, 0) -= s(2, 1);
678 10 : d2(2, 1, 0, 1) += s(2, 0) / 2.0;
679 10 : d2(2, 1, 1, 0) += s(2, 0) / 2.0;
680 10 : d2(2, 1, 0, 2) += s(0, 1) / 2.0;
681 10 : d2(2, 1, 2, 0) += s(0, 1) / 2.0;
682 10 : d2(2, 1, 1, 2) -= s(0, 0) / 2.0;
683 10 : d2(2, 1, 2, 1) -= s(0, 0) / 2.0;
684 :
685 10 : d2(2, 2, 0, 0) += s(1, 1);
686 10 : d2(2, 2, 0, 1) -= s(1, 0);
687 10 : d2(2, 2, 1, 0) -= s(1, 0);
688 10 : d2(2, 2, 1, 1) += s(0, 0);
689 :
690 20 : return d2;
691 10 : }
692 :
693 : template <typename T>
694 : RankTwoTensorTempl<T>
695 6 : RankTwoTensorTempl<T>::ddet() const
696 : {
697 6 : RankTwoTensorTempl<T> d;
698 :
699 6 : d(0, 0) = (*this)(1, 1) * (*this)(2, 2) - (*this)(2, 1) * (*this)(1, 2);
700 6 : d(0, 1) = (*this)(2, 0) * (*this)(1, 2) - (*this)(1, 0) * (*this)(2, 2);
701 6 : d(0, 2) = (*this)(1, 0) * (*this)(2, 1) - (*this)(2, 0) * (*this)(1, 1);
702 6 : d(1, 0) = (*this)(2, 1) * (*this)(0, 2) - (*this)(0, 1) * (*this)(2, 2);
703 6 : d(1, 1) = (*this)(0, 0) * (*this)(2, 2) - (*this)(2, 0) * (*this)(0, 2);
704 6 : d(1, 2) = (*this)(2, 0) * (*this)(0, 1) - (*this)(0, 0) * (*this)(2, 1);
705 6 : d(2, 0) = (*this)(0, 1) * (*this)(1, 2) - (*this)(1, 1) * (*this)(0, 2);
706 6 : d(2, 1) = (*this)(1, 0) * (*this)(0, 2) - (*this)(0, 0) * (*this)(1, 2);
707 6 : d(2, 2) = (*this)(0, 0) * (*this)(1, 1) - (*this)(1, 0) * (*this)(0, 1);
708 :
709 6 : return d;
710 0 : }
711 :
712 : template <typename T>
713 : void
714 0 : RankTwoTensorTempl<T>::print(std::ostream & stm) const
715 : {
716 0 : const RankTwoTensorTempl<T> & a = *this;
717 0 : for (const auto i : make_range(N))
718 : {
719 0 : for (const auto j : make_range(N))
720 0 : stm << std::setw(15) << a(i, j) << ' ';
721 0 : stm << std::endl;
722 : }
723 0 : }
724 :
725 : template <>
726 : void RankTwoTensor::printReal(std::ostream & stm) const;
727 :
728 : template <>
729 : void ADRankTwoTensor::printReal(std::ostream & stm) const;
730 :
731 : template <>
732 : void ADRankTwoTensor::printADReal(unsigned int nDual, std::ostream & stm) const;
733 :
734 : template <typename T>
735 : void
736 2654 : RankTwoTensorTempl<T>::addIa(const T & a)
737 : {
738 10616 : for (const auto i : make_range(N))
739 7962 : (*this)(i, i) += a;
740 2654 : }
741 :
742 : template <typename T>
743 : T
744 10971794 : RankTwoTensorTempl<T>::L2norm() const
745 : {
746 10971794 : T norm = 0.0;
747 109717940 : for (const auto i : make_range(N2))
748 : {
749 98746146 : T v = _coords[i];
750 98746146 : norm += v * v;
751 : }
752 : using std::sqrt;
753 10971796 : return norm == 0.0 ? 0.0 : sqrt(norm);
754 2 : }
755 :
756 : template <typename T>
757 : void
758 0 : RankTwoTensorTempl<T>::surfaceFillFromInputVector(const std::vector<T> & input)
759 : {
760 0 : if (input.size() == 4)
761 : {
762 : // initialize with zeros
763 0 : this->zero();
764 0 : (*this)(0, 0) = input[0];
765 0 : (*this)(0, 1) = input[1];
766 0 : (*this)(1, 0) = input[2];
767 0 : (*this)(1, 1) = input[3];
768 : }
769 : else
770 0 : mooseError("please provide correct number of values for surface RankTwoTensorTempl<T> "
771 : "initialization.");
772 0 : }
773 :
774 : template <typename T>
775 : void
776 0 : RankTwoTensorTempl<T>::syev(const char *, std::vector<T> &, std::vector<T> &) const
777 : {
778 0 : mooseError("The syev method is only supported for Real valued tensors");
779 : }
780 :
781 : template <>
782 : void RankTwoTensor::syev(const char * calculation_type,
783 : std::vector<Real> & eigvals,
784 : std::vector<Real> & a) const;
785 :
786 : template <typename T>
787 : void
788 0 : RankTwoTensorTempl<T>::symmetricEigenvalues(std::vector<T> & eigvals) const
789 : {
790 0 : RankTwoTensorTempl<T> a;
791 0 : symmetricEigenvaluesEigenvectors(eigvals, a);
792 0 : }
793 :
794 : template <>
795 : void RankTwoTensor::symmetricEigenvalues(std::vector<Real> & eigvals) const;
796 :
797 : template <typename T>
798 : void
799 : RankTwoTensorTempl<T>::symmetricEigenvaluesEigenvectors(std::vector<T> &,
800 : RankTwoTensorTempl<T> &) const
801 : {
802 : mooseError(
803 : "symmetricEigenvaluesEigenvectors is only available for ordered tensor component types");
804 : }
805 :
806 : template <>
807 : void RankTwoTensor::symmetricEigenvaluesEigenvectors(std::vector<Real> & eigvals,
808 : RankTwoTensor & eigvecs) const;
809 :
810 : template <>
811 : void ADRankTwoTensor::symmetricEigenvaluesEigenvectors(std::vector<ADReal> & eigvals,
812 : ADRankTwoTensor & eigvecs) const;
813 :
814 : template <typename T>
815 : void
816 988 : RankTwoTensorTempl<T>::dsymmetricEigenvalues(std::vector<T> & eigvals,
817 : std::vector<RankTwoTensorTempl<T>> & deigvals) const
818 : {
819 988 : deigvals.resize(N);
820 :
821 988 : std::vector<T> a;
822 988 : syev("V", eigvals, a);
823 :
824 : // now a contains the eigenvetors
825 : // extract these and place appropriately in deigvals
826 988 : std::vector<T> eig_vec;
827 988 : eig_vec.resize(N);
828 :
829 3952 : for (const auto i : make_range(N))
830 : {
831 11856 : for (const auto j : make_range(N))
832 8892 : eig_vec[j] = a[i * N + j];
833 11856 : for (const auto j : make_range(N))
834 35568 : for (const auto k : make_range(N))
835 26676 : deigvals[i](j, k) = eig_vec[j] * eig_vec[k];
836 : }
837 :
838 : // There are discontinuities in the derivative
839 : // for equal eigenvalues. The following is
840 : // an attempt to make a sensible choice for
841 : // the derivative. This agrees with a central-difference
842 : // approximation to the derivative.
843 988 : if (eigvals[0] == eigvals[1] && eigvals[0] == eigvals[2])
844 0 : deigvals[0] = deigvals[1] = deigvals[2] = (deigvals[0] + deigvals[1] + deigvals[2]) / 3.0;
845 988 : else if (eigvals[0] == eigvals[1])
846 4 : deigvals[0] = deigvals[1] = (deigvals[0] + deigvals[1]) / 2.0;
847 984 : else if (eigvals[0] == eigvals[2])
848 0 : deigvals[0] = deigvals[2] = (deigvals[0] + deigvals[2]) / 2.0;
849 984 : else if (eigvals[1] == eigvals[2])
850 4 : deigvals[1] = deigvals[2] = (deigvals[1] + deigvals[2]) / 2.0;
851 988 : }
852 :
853 : template <typename T>
854 : void
855 8 : RankTwoTensorTempl<T>::d2symmetricEigenvalues(std::vector<RankFourTensorTempl<T>> & deriv) const
856 : {
857 8 : std::vector<T> eigvec;
858 8 : std::vector<T> eigvals;
859 0 : T ev[N][N];
860 :
861 : // reset rank four tensor
862 8 : deriv.assign(N, RankFourTensorTempl<T>());
863 :
864 : // get eigen values and eigen vectors
865 8 : syev("V", eigvals, eigvec);
866 :
867 32 : for (const auto i : make_range(N))
868 96 : for (const auto j : make_range(N))
869 72 : ev[i][j] = eigvec[i * N + j];
870 :
871 32 : for (unsigned int alpha = 0; alpha < N; ++alpha)
872 96 : for (unsigned int beta = 0; beta < N; ++beta)
873 : {
874 72 : if (eigvals[alpha] == eigvals[beta])
875 24 : continue;
876 :
877 192 : for (const auto i : make_range(N))
878 576 : for (const auto j : make_range(N))
879 1728 : for (const auto k : make_range(N))
880 5184 : for (const auto l : make_range(N))
881 : {
882 3888 : deriv[alpha](i, j, k, l) +=
883 3888 : 0.5 * (ev[beta][i] * ev[alpha][j] + ev[alpha][i] * ev[beta][j]) *
884 7776 : (ev[beta][k] * ev[alpha][l] + ev[beta][l] * ev[alpha][k]) /
885 3888 : (eigvals[alpha] - eigvals[beta]);
886 : }
887 : }
888 8 : }
889 :
890 : template <typename T>
891 : void
892 0 : RankTwoTensorTempl<T>::getRUDecompositionRotation(RankTwoTensorTempl<T> &) const
893 : {
894 0 : mooseError("getRUDecompositionRotation is only supported for Real valued tensors");
895 : }
896 :
897 : template <>
898 : void RankTwoTensor::getRUDecompositionRotation(RankTwoTensor & rot) const;
899 :
900 : template <typename T>
901 : void
902 2 : RankTwoTensorTempl<T>::initRandom(unsigned int rand_seed)
903 : {
904 2 : MooseRandom::seed(rand_seed);
905 2 : }
906 :
907 : template <typename T>
908 : RankTwoTensorTempl<T>
909 0 : RankTwoTensorTempl<T>::genRandomTensor(T stddev, T mean)
910 : {
911 0 : RankTwoTensorTempl<T> result;
912 0 : for (const auto i : make_range(N))
913 0 : for (const auto j : make_range(N))
914 0 : result(i, j) = (MooseRandom::rand() + mean) * stddev;
915 0 : return result;
916 0 : }
917 :
918 : template <typename T>
919 : RankTwoTensorTempl<T>
920 2 : RankTwoTensorTempl<T>::genRandomSymmTensor(T stddev, T mean)
921 : {
922 14 : auto r = [&]() { return (MooseRandom::rand() + mean) * stddev; };
923 2 : return RankTwoTensorTempl<T>(r(), r(), r(), r(), r(), r());
924 : }
925 :
926 : template <typename T>
927 : void
928 0 : RankTwoTensorTempl<T>::vectorOuterProduct(const libMesh::TypeVector<T> & v1,
929 : const libMesh::TypeVector<T> & v2)
930 : {
931 0 : RankTwoTensorTempl<T> & a = *this;
932 0 : for (const auto i : make_range(N))
933 0 : for (const auto j : make_range(N))
934 0 : a(i, j) = v1(i) * v2(j);
935 0 : }
936 :
937 : template <typename T>
938 : RankTwoTensorTempl<T>
939 6 : RankTwoTensorTempl<T>::outerProduct(const libMesh::TypeVector<T> & v1,
940 : const libMesh::TypeVector<T> & v2)
941 : {
942 6 : RankTwoTensorTempl<T> result;
943 24 : for (const auto i : make_range(N))
944 72 : for (const auto j : make_range(N))
945 54 : result(i, j) = v1(i) * v2(j);
946 6 : return result;
947 0 : }
948 :
949 : template <typename T>
950 : RankTwoTensorTempl<T>
951 1216 : RankTwoTensorTempl<T>::selfOuterProduct(const libMesh::TypeVector<T> & v)
952 : {
953 1216 : RankTwoTensorTempl<T> result(RankTwoTensorTempl<T>::initNone);
954 4864 : for (unsigned int i = 0; i < N; ++i)
955 14592 : for (unsigned int j = 0; j < N; ++j)
956 10944 : result(i, j) = v(i) * v(j);
957 1216 : return result;
958 0 : }
959 :
960 : template <typename T>
961 : void
962 2 : RankTwoTensorTempl<T>::fillRealTensor(libMesh::TensorValue<T> & tensor)
963 : {
964 8 : for (const auto i : make_range(N))
965 24 : for (const auto j : make_range(N))
966 18 : tensor(i, j) = (*this)(i, j);
967 2 : }
968 :
969 : template <typename T>
970 : void
971 2 : RankTwoTensorTempl<T>::fillRow(unsigned int r, const libMesh::TypeVector<T> & v)
972 : {
973 8 : for (const auto i : make_range(N))
974 6 : (*this)(r, i) = v(i);
975 2 : }
976 :
977 : template <typename T>
978 : void
979 2 : RankTwoTensorTempl<T>::fillColumn(unsigned int c, const libMesh::TypeVector<T> & v)
980 : {
981 8 : for (const auto i : make_range(N))
982 6 : (*this)(i, c) = v(i);
983 2 : }
984 :
985 : template <typename T>
986 : RankTwoTensorTempl<T>
987 4 : RankTwoTensorTempl<T>::initialContraction(const RankFourTensorTempl<T> & b) const
988 : {
989 4 : RankTwoTensorTempl<T> result;
990 16 : for (const auto i : make_range(N))
991 48 : for (const auto j : make_range(N))
992 144 : for (const auto k : make_range(N))
993 432 : for (const auto l : make_range(N))
994 324 : result(k, l) += (*this)(i, j) * b(i, j, k, l);
995 4 : return result;
996 0 : }
997 :
998 : template <typename T>
999 : void
1000 0 : RankTwoTensorTempl<T>::setToIdentity()
1001 : {
1002 : mooseAssert(N2 == 9, "RankTwoTensorTempl is currently only tested for 3 dimensions.");
1003 0 : for (const auto i : make_range(N2))
1004 0 : _coords[i] = identityCoords[i];
1005 0 : }
|
