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 3029489 : RankTwoTensorTempl<T>::RankTwoTensorTempl()
73 : {
74 : mooseAssert(N == 3, "RankTwoTensorTempl is currently only tested for 3 dimensions.");
75 :
76 30294890 : for (unsigned int i = 0; i < N2; i++)
77 27265401 : _coords[i] = 0.0;
78 3029489 : }
79 :
80 : template <typename T>
81 608 : RankTwoTensorTempl<T>::RankTwoTensorTempl(const InitMethod init)
82 : {
83 608 : switch (init)
84 : {
85 608 : case initNone:
86 608 : 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 608 : }
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 2 : 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 2 : (v1(0) + v0(1)) / 2.0,
126 2 : (v2(0) + v0(2)) / 2.0,
127 2 : (v1(0) + v0(1)) / 2.0,
128 : v1(1),
129 2 : (v2(1) + v1(2)) / 2.0,
130 2 : (v2(0) + v0(2)) / 2.0,
131 4 : (v2(1) + v1(2)) / 2.0,
132 8 : v2(2));
133 : }
134 :
135 : template <typename T>
136 : RankTwoTensorTempl<T>
137 1 : 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 1 : 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 14 : RankTwoTensorTempl<T>::RankTwoTensorTempl(
157 14 : const T & S11, const T & S22, const T & S33, const T & S23, const T & S13, const T & S12)
158 : {
159 14 : (*this)(0, 0) = S11;
160 14 : (*this)(1, 1) = S22;
161 14 : (*this)(2, 2) = S33;
162 14 : (*this)(1, 2) = (*this)(2, 1) = S23;
163 14 : (*this)(0, 2) = (*this)(2, 0) = S13;
164 14 : (*this)(0, 1) = (*this)(1, 0) = S12;
165 14 : }
166 :
167 : template <typename T>
168 372 : 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 372 : const T & S33)
177 : {
178 372 : (*this)(0, 0) = S11;
179 372 : (*this)(1, 0) = S21;
180 372 : (*this)(2, 0) = S31;
181 372 : (*this)(0, 1) = S12;
182 372 : (*this)(1, 1) = S22;
183 372 : (*this)(2, 1) = S32;
184 372 : (*this)(0, 2) = S13;
185 372 : (*this)(1, 2) = S23;
186 372 : (*this)(2, 2) = S33;
187 372 : }
188 :
189 : template <typename T>
190 : void
191 217 : RankTwoTensorTempl<T>::fillFromInputVector(const std::vector<T> & input, FillMethod fill_method)
192 : {
193 217 : 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 217 : switch (input.size())
197 : {
198 1 : case 1:
199 1 : this->zero();
200 1 : (*this)(0, 0) = input[0]; // S11
201 1 : (*this)(1, 1) = input[0]; // S22
202 1 : (*this)(2, 2) = input[0]; // S33
203 1 : break;
204 :
205 7 : case 3:
206 7 : this->zero();
207 7 : (*this)(0, 0) = input[0]; // S11
208 7 : (*this)(1, 1) = input[1]; // S22
209 7 : (*this)(2, 2) = input[2]; // S33
210 7 : break;
211 :
212 40 : case 6:
213 40 : (*this)(0, 0) = input[0]; // S11
214 40 : (*this)(1, 1) = input[1]; // S22
215 40 : (*this)(2, 2) = input[2]; // S33
216 40 : (*this)(1, 2) = (*this)(2, 1) = input[3]; // S23
217 40 : (*this)(0, 2) = (*this)(2, 0) = input[4]; // S13
218 40 : (*this)(0, 1) = (*this)(1, 0) = input[5]; // S12
219 40 : break;
220 :
221 169 : case 9:
222 169 : (*this)(0, 0) = input[0]; // S11
223 169 : (*this)(1, 0) = input[1]; // S21
224 169 : (*this)(2, 0) = input[2]; // S31
225 169 : (*this)(0, 1) = input[3]; // S12
226 169 : (*this)(1, 1) = input[4]; // S22
227 169 : (*this)(2, 1) = input[5]; // S32
228 169 : (*this)(0, 2) = input[6]; // S13
229 169 : (*this)(1, 2) = input[7]; // S23
230 169 : (*this)(2, 2) = input[8]; // S33
231 169 : 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 217 : }
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 615 : RankTwoTensorTempl<T>::column(const unsigned int c) const
275 : {
276 615 : return libMesh::VectorValue<T>((*this)(0, c), (*this)(1, c), (*this)(2, c));
277 : }
278 :
279 : template <typename T>
280 : RankTwoTensorTempl<T>
281 2 : RankTwoTensorTempl<T>::timesTranspose(const RankTwoTensorTempl<T> & a)
282 : {
283 2 : return a * a.transpose();
284 : }
285 :
286 : template <typename T>
287 : RankTwoTensorTempl<T>
288 1 : RankTwoTensorTempl<T>::transposeTimes(const RankTwoTensorTempl<T> & a)
289 : {
290 1 : return a.transpose() * a;
291 : }
292 :
293 : template <typename T>
294 : RankTwoTensorTempl<T>
295 1320 : RankTwoTensorTempl<T>::plusTranspose(const RankTwoTensorTempl<T> & a)
296 : {
297 1320 : return a + a.transpose();
298 : }
299 :
300 : template <typename T>
301 : RankTwoTensorTempl<T>
302 2 : RankTwoTensorTempl<T>::square() const
303 : {
304 2 : return *this * *this;
305 : }
306 :
307 : template <typename T>
308 : RankTwoTensorTempl<T>
309 9 : RankTwoTensorTempl<T>::rotated(const RankTwoTensorTempl<T> & R) const
310 : {
311 9 : RankTwoTensorTempl<T> result(*this);
312 9 : result.rotate(R);
313 9 : return result;
314 0 : }
315 :
316 : template <typename T>
317 : void
318 17 : RankTwoTensorTempl<T>::rotate(const RankTwoTensorTempl<T> & R)
319 : {
320 17 : RankTwoTensorTempl<T> temp;
321 68 : for (const auto i : make_range(N))
322 : {
323 51 : const auto i1 = i * N;
324 204 : for (const auto j : make_range(N))
325 : {
326 : // tmp += R(i,k)*R(j,l)*(*this)(k,l);
327 : // clang-format off
328 153 : const auto j1 = j * N;
329 153 : T tmp = R._coords[i1 + 0] * R._coords[j1 + 0] * (*this)(0, 0) +
330 153 : R._coords[i1 + 0] * R._coords[j1 + 1] * (*this)(0, 1) +
331 153 : R._coords[i1 + 0] * R._coords[j1 + 2] * (*this)(0, 2) +
332 153 : R._coords[i1 + 1] * R._coords[j1 + 0] * (*this)(1, 0) +
333 153 : R._coords[i1 + 1] * R._coords[j1 + 1] * (*this)(1, 1) +
334 153 : R._coords[i1 + 1] * R._coords[j1 + 2] * (*this)(1, 2) +
335 153 : R._coords[i1 + 2] * R._coords[j1 + 0] * (*this)(2, 0) +
336 153 : R._coords[i1 + 2] * R._coords[j1 + 1] * (*this)(2, 1) +
337 153 : R._coords[i1 + 2] * R._coords[j1 + 2] * (*this)(2, 2);
338 : // clang-format on
339 153 : temp._coords[i1 + j] = tmp;
340 : }
341 : }
342 170 : for (unsigned int i = 0; i < N2; i++)
343 153 : _coords[i] = temp._coords[i];
344 17 : }
345 :
346 : template <typename T>
347 : RankTwoTensorTempl<T>
348 0 : RankTwoTensorTempl<T>::rotateXyPlane(T a)
349 : {
350 0 : T c = std::cos(a);
351 0 : T s = std::sin(a);
352 0 : T x = (*this)(0, 0) * c * c + (*this)(1, 1) * s * s + 2.0 * (*this)(0, 1) * c * s;
353 0 : T y = (*this)(0, 0) * s * s + (*this)(1, 1) * c * c - 2.0 * (*this)(0, 1) * c * s;
354 0 : T xy = ((*this)(1, 1) - (*this)(0, 0)) * c * s + (*this)(0, 1) * (c * c - s * s);
355 :
356 0 : RankTwoTensorTempl<T> b(*this);
357 :
358 0 : b(0, 0) = x;
359 0 : b(1, 1) = y;
360 0 : b(1, 0) = b(0, 1) = xy;
361 :
362 0 : return b;
363 0 : }
364 :
365 : template <typename T>
366 : RankTwoTensorTempl<T>
367 1406 : RankTwoTensorTempl<T>::transpose() const
368 : {
369 1406 : return libMesh::TensorValue<T>::transpose();
370 : }
371 :
372 : template <typename T>
373 : RankTwoTensorTempl<T> &
374 1 : RankTwoTensorTempl<T>::operator+=(const RankTwoTensorTempl<T> & a)
375 : {
376 1 : libMesh::TensorValue<T>::operator+=(a);
377 1 : return *this;
378 : }
379 :
380 : template <typename T>
381 : RankTwoTensorTempl<T> &
382 1 : RankTwoTensorTempl<T>::operator-=(const RankTwoTensorTempl<T> & a)
383 : {
384 1 : libMesh::TensorValue<T>::operator-=(a);
385 1 : return *this;
386 : }
387 :
388 : template <typename T>
389 : RankTwoTensorTempl<T>
390 1 : RankTwoTensorTempl<T>::operator-() const
391 : {
392 1 : return libMesh::TensorValue<T>::operator-();
393 : }
394 :
395 : template <typename T>
396 : RankTwoTensorTempl<T> &
397 2 : RankTwoTensorTempl<T>::operator*=(const T & a)
398 : {
399 2 : libMesh::TensorValue<T>::operator*=(a);
400 2 : return *this;
401 : }
402 :
403 : template <typename T>
404 : RankTwoTensorTempl<T> &
405 2 : RankTwoTensorTempl<T>::operator/=(const T & a)
406 : {
407 2 : libMesh::TensorValue<T>::operator/=(a);
408 2 : return *this;
409 : }
410 :
411 : template <typename T>
412 : RankTwoTensorTempl<T> &
413 0 : RankTwoTensorTempl<T>::operator*=(const libMesh::TypeTensor<T> & a)
414 : {
415 0 : *this = *this * a;
416 0 : return *this;
417 : }
418 :
419 : template <typename T>
420 : bool
421 3 : RankTwoTensorTempl<T>::operator==(const RankTwoTensorTempl<T> & a) const
422 : {
423 9 : for (const auto i : make_range(N))
424 26 : for (const auto j : make_range(N))
425 20 : if (!MooseUtils::absoluteFuzzyEqual((*this)(i, j), a(i, j)))
426 1 : return false;
427 :
428 2 : return true;
429 : }
430 :
431 : template <typename T>
432 : bool
433 23 : RankTwoTensorTempl<T>::isSymmetric() const
434 : {
435 23 : auto test = MetaPhysicL::raw_value(*this - transpose());
436 46 : return MooseUtils::absoluteFuzzyEqual(test.norm_sq(), 0);
437 23 : }
438 :
439 : template <typename T>
440 : RankTwoTensorTempl<T> &
441 0 : RankTwoTensorTempl<T>::operator=(const ColumnMajorMatrixTempl<T> & a)
442 : {
443 0 : if (a.n() != N || a.m() != N)
444 0 : mooseError("Dimensions of ColumnMajorMatrixTempl<T> are incompatible with RankTwoTensorTempl");
445 :
446 0 : const T * cmm_rawdata = a.rawData();
447 0 : for (const auto i : make_range(N))
448 0 : for (const auto j : make_range(N))
449 0 : _coords[i * N + j] = cmm_rawdata[i + j * N];
450 :
451 0 : return *this;
452 : }
453 :
454 : template <typename T>
455 : T
456 2 : RankTwoTensorTempl<T>::doubleContraction(const RankTwoTensorTempl<T> & b) const
457 : {
458 : // deprecate this!
459 2 : return libMesh::TensorValue<T>::contract(b);
460 : }
461 :
462 : template <typename T>
463 : RankThreeTensorTempl<T>
464 1 : RankTwoTensorTempl<T>::contraction(const RankThreeTensorTempl<T> & b) const
465 : {
466 1 : RankThreeTensorTempl<T> result;
467 :
468 4 : for (const auto i : make_range(N))
469 12 : for (const auto j : make_range(N))
470 36 : for (const auto k : make_range(N))
471 108 : for (const auto l : make_range(N))
472 81 : result(i, k, l) += (*this)(i, j) * b(j, k, l);
473 :
474 1 : return result;
475 0 : }
476 :
477 : template <typename T>
478 : RankThreeTensorTempl<T>
479 1 : RankTwoTensorTempl<T>::mixedProductJkI(const libMesh::VectorValue<T> & b) const
480 : {
481 1 : RankThreeTensorTempl<T> result;
482 :
483 4 : for (const auto i : make_range(N))
484 12 : for (const auto j : make_range(N))
485 36 : for (const auto k : make_range(N))
486 27 : result(i, j, k) += (*this)(j, k) * b(i);
487 :
488 1 : return result;
489 0 : }
490 :
491 : template <typename T>
492 : RankTwoTensorTempl<T>
493 1325 : RankTwoTensorTempl<T>::deviatoric() const
494 : {
495 1325 : RankTwoTensorTempl<T> deviatoric(*this);
496 : // actually construct deviatoric part
497 1325 : deviatoric.addIa(-1.0 / 3.0 * this->tr());
498 1325 : return deviatoric;
499 0 : }
500 :
501 : template <typename T>
502 : T
503 2 : RankTwoTensorTempl<T>::generalSecondInvariant() const
504 : {
505 2 : return (*this)(0, 0) * (*this)(1, 1) + (*this)(0, 0) * (*this)(2, 2) +
506 2 : (*this)(1, 1) * (*this)(2, 2) - (*this)(0, 1) * (*this)(1, 0) -
507 2 : (*this)(0, 2) * (*this)(2, 0) - (*this)(1, 2) * (*this)(2, 1);
508 : }
509 :
510 : template <typename T>
511 : T
512 786 : RankTwoTensorTempl<T>::secondInvariant() const
513 : {
514 0 : T result;
515 :
516 : // RankTwoTensorTempl<T> deviatoric(*this);
517 : // deviatoric.addIa(-1.0/3.0 * this->tr()); // actually construct deviatoric part
518 : // result = 0.5*(deviatoric + deviatoric.transpose()).doubleContraction(deviatoric +
519 : // deviatoric.transpose());
520 786 : result = Utility::pow<2>((*this)(0, 0) - (*this)(1, 1)) / 6.0;
521 786 : result += Utility::pow<2>((*this)(0, 0) - (*this)(2, 2)) / 6.0;
522 786 : result += Utility::pow<2>((*this)(1, 1) - (*this)(2, 2)) / 6.0;
523 786 : result += Utility::pow<2>((*this)(0, 1) + (*this)(1, 0)) / 4.0;
524 786 : result += Utility::pow<2>((*this)(0, 2) + (*this)(2, 0)) / 4.0;
525 786 : result += Utility::pow<2>((*this)(1, 2) + (*this)(2, 1)) / 4.0;
526 786 : return result;
527 0 : }
528 :
529 : template <typename T>
530 : RankTwoTensorTempl<T>
531 500 : RankTwoTensorTempl<T>::dsecondInvariant() const
532 : {
533 500 : return RankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
534 : }
535 :
536 : template <typename T>
537 : RankFourTensorTempl<T>
538 5 : RankTwoTensorTempl<T>::d2secondInvariant() const
539 : {
540 5 : RankFourTensorTempl<T> result;
541 20 : for (const auto i : make_range(N))
542 60 : for (const auto j : make_range(N))
543 180 : for (const auto k : make_range(N))
544 540 : for (const auto l : make_range(N))
545 810 : result(i, j, k, l) = 0.5 * (i == k) * (j == l) + 0.5 * (i == l) * (j == k) -
546 405 : (1.0 / 3.0) * (i == j) * (k == l);
547 5 : return result;
548 0 : }
549 :
550 : template <typename T>
551 : T
552 8640001 : RankTwoTensorTempl<T>::trace() const
553 : {
554 : // deprecate this!
555 8640001 : return this->tr();
556 : }
557 :
558 : template <typename T>
559 : RankTwoTensorTempl<T>
560 4 : RankTwoTensorTempl<T>::dtrace() const
561 : {
562 4 : return RankTwoTensorTempl<T>(1, 0, 0, 0, 1, 0, 0, 0, 1);
563 : }
564 :
565 : template <typename T>
566 : RankTwoTensorTempl<T>
567 6 : RankTwoTensorTempl<T>::inverse() const
568 : {
569 6 : return libMesh::TensorValue<T>::inverse();
570 : }
571 :
572 : template <typename T>
573 : T
574 315 : RankTwoTensorTempl<T>::thirdInvariant() const
575 : {
576 315 : const auto s = RankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
577 315 : return s(0, 0) * (s(1, 1) * s(2, 2) - s(2, 1) * s(1, 2)) -
578 315 : s(1, 0) * (s(0, 1) * s(2, 2) - s(2, 1) * s(0, 2)) +
579 630 : s(2, 0) * (s(0, 1) * s(1, 2) - s(1, 1) * s(0, 2));
580 315 : }
581 :
582 : template <typename T>
583 : RankTwoTensorTempl<T>
584 498 : RankTwoTensorTempl<T>::dthirdInvariant() const
585 : {
586 498 : const auto s = RankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
587 498 : const T s3 = secondInvariant() / 3.0;
588 :
589 498 : RankTwoTensorTempl<T> d;
590 498 : d(0, 0) = s(1, 1) * s(2, 2) - s(2, 1) * s(1, 2) + s3;
591 498 : d(0, 1) = s(2, 0) * s(1, 2) - s(1, 0) * s(2, 2);
592 498 : d(0, 2) = s(1, 0) * s(2, 1) - s(2, 0) * s(1, 1);
593 498 : d(1, 0) = s(2, 1) * s(0, 2) - s(0, 1) * s(2, 2);
594 498 : d(1, 1) = s(0, 0) * s(2, 2) - s(2, 0) * s(0, 2) + s3;
595 498 : d(1, 2) = s(2, 0) * s(0, 1) - s(0, 0) * s(2, 1);
596 498 : d(2, 0) = s(0, 1) * s(1, 2) - s(1, 1) * s(0, 2);
597 498 : d(2, 1) = s(1, 0) * s(0, 2) - s(0, 0) * s(1, 2);
598 498 : d(2, 2) = s(0, 0) * s(1, 1) - s(1, 0) * s(0, 1) + s3;
599 996 : return d;
600 498 : }
601 :
602 : template <typename T>
603 : RankFourTensorTempl<T>
604 5 : RankTwoTensorTempl<T>::d2thirdInvariant() const
605 : {
606 5 : const auto s = RankTwoTensorTempl<T>::plusTranspose(deviatoric()) * 0.5;
607 :
608 5 : RankFourTensorTempl<T> d2;
609 20 : for (const auto i : make_range(N))
610 60 : for (const auto j : make_range(N))
611 180 : for (const auto k : make_range(N))
612 540 : for (const auto l : make_range(N))
613 : {
614 405 : d2(i, j, k, l) = Real(i == j) * s(k, l) / 3.0 + Real(k == l) * s(i, j) / 3.0;
615 : // for (const auto a: make_range(N))
616 : // for (const auto b: make_range(N))
617 : // d2(i, j, k, l) += 0.5*(PermutationTensor::eps(i, k, a)*PermutationTensor::eps(j, l,
618 : // b) + PermutationTensor::eps(i, l, a)*PermutationTensor::eps(j, k, b))*s(a, b);
619 : }
620 :
621 : // I'm not sure which is more readable: the above
622 : // PermutationTensor stuff, or the stuff below.
623 : // Anyway, they yield the same result, and so i leave
624 : // both of them here to enlighten you!
625 :
626 5 : d2(0, 0, 1, 1) += s(2, 2);
627 5 : d2(0, 0, 1, 2) -= s(2, 1);
628 5 : d2(0, 0, 2, 1) -= s(1, 2);
629 5 : d2(0, 0, 2, 2) += s(1, 1);
630 :
631 5 : d2(0, 1, 0, 1) -= s(2, 2) / 2.0;
632 5 : d2(0, 1, 1, 0) -= s(2, 2) / 2.0;
633 5 : d2(0, 1, 0, 2) += s(1, 2) / 2.0;
634 5 : d2(0, 1, 2, 0) += s(1, 2) / 2.0;
635 5 : d2(0, 1, 1, 2) += s(2, 0) / 2.0;
636 5 : d2(0, 1, 2, 1) += s(2, 0) / 2.0;
637 5 : d2(0, 1, 2, 2) -= s(1, 0);
638 :
639 5 : d2(0, 2, 0, 1) += s(2, 1) / 2.0;
640 5 : d2(0, 2, 1, 0) += s(2, 1) / 2.0;
641 5 : d2(0, 2, 0, 2) -= s(1, 1) / 2.0;
642 5 : d2(0, 2, 2, 0) -= s(1, 1) / 2.0;
643 5 : d2(0, 2, 1, 1) -= s(2, 0);
644 5 : d2(0, 2, 1, 2) += s(1, 0) / 2.0;
645 5 : d2(0, 2, 2, 1) += s(1, 0) / 2.0;
646 :
647 5 : d2(1, 0, 0, 1) -= s(2, 2) / 2.0;
648 5 : d2(1, 0, 1, 0) -= s(2, 2) / 2.0;
649 5 : d2(1, 0, 0, 2) += s(1, 2) / 2.0;
650 5 : d2(1, 0, 2, 0) += s(1, 2) / 2.0;
651 5 : d2(1, 0, 1, 2) += s(2, 0) / 2.0;
652 5 : d2(1, 0, 2, 1) += s(2, 0) / 2.0;
653 5 : d2(1, 0, 2, 2) -= s(1, 0);
654 :
655 5 : d2(1, 1, 0, 0) += s(2, 2);
656 5 : d2(1, 1, 0, 2) -= s(2, 0);
657 5 : d2(1, 1, 2, 0) -= s(2, 0);
658 5 : d2(1, 1, 2, 2) += s(0, 0);
659 :
660 5 : d2(1, 2, 0, 0) -= s(2, 1);
661 5 : d2(1, 2, 0, 1) += s(2, 0) / 2.0;
662 5 : d2(1, 2, 1, 0) += s(2, 0) / 2.0;
663 5 : d2(1, 2, 0, 2) += s(0, 1) / 2.0;
664 5 : d2(1, 2, 2, 0) += s(0, 1) / 2.0;
665 5 : d2(1, 2, 1, 2) -= s(0, 0) / 2.0;
666 5 : d2(1, 2, 2, 1) -= s(0, 0) / 2.0;
667 :
668 5 : d2(2, 0, 0, 1) += s(2, 1) / 2.0;
669 5 : d2(2, 0, 1, 0) += s(2, 1) / 2.0;
670 5 : d2(2, 0, 0, 2) -= s(1, 1) / 2.0;
671 5 : d2(2, 0, 2, 0) -= s(1, 1) / 2.0;
672 5 : d2(2, 0, 1, 1) -= s(2, 0);
673 5 : d2(2, 0, 1, 2) += s(1, 0) / 2.0;
674 5 : d2(2, 0, 2, 1) += s(1, 0) / 2.0;
675 :
676 5 : d2(2, 1, 0, 0) -= s(2, 1);
677 5 : d2(2, 1, 0, 1) += s(2, 0) / 2.0;
678 5 : d2(2, 1, 1, 0) += s(2, 0) / 2.0;
679 5 : d2(2, 1, 0, 2) += s(0, 1) / 2.0;
680 5 : d2(2, 1, 2, 0) += s(0, 1) / 2.0;
681 5 : d2(2, 1, 1, 2) -= s(0, 0) / 2.0;
682 5 : d2(2, 1, 2, 1) -= s(0, 0) / 2.0;
683 :
684 5 : d2(2, 2, 0, 0) += s(1, 1);
685 5 : d2(2, 2, 0, 1) -= s(1, 0);
686 5 : d2(2, 2, 1, 0) -= s(1, 0);
687 5 : d2(2, 2, 1, 1) += s(0, 0);
688 :
689 10 : return d2;
690 5 : }
691 :
692 : template <typename T>
693 : RankTwoTensorTempl<T>
694 3 : RankTwoTensorTempl<T>::ddet() const
695 : {
696 3 : RankTwoTensorTempl<T> d;
697 :
698 3 : d(0, 0) = (*this)(1, 1) * (*this)(2, 2) - (*this)(2, 1) * (*this)(1, 2);
699 3 : d(0, 1) = (*this)(2, 0) * (*this)(1, 2) - (*this)(1, 0) * (*this)(2, 2);
700 3 : d(0, 2) = (*this)(1, 0) * (*this)(2, 1) - (*this)(2, 0) * (*this)(1, 1);
701 3 : d(1, 0) = (*this)(2, 1) * (*this)(0, 2) - (*this)(0, 1) * (*this)(2, 2);
702 3 : d(1, 1) = (*this)(0, 0) * (*this)(2, 2) - (*this)(2, 0) * (*this)(0, 2);
703 3 : d(1, 2) = (*this)(2, 0) * (*this)(0, 1) - (*this)(0, 0) * (*this)(2, 1);
704 3 : d(2, 0) = (*this)(0, 1) * (*this)(1, 2) - (*this)(1, 1) * (*this)(0, 2);
705 3 : d(2, 1) = (*this)(1, 0) * (*this)(0, 2) - (*this)(0, 0) * (*this)(1, 2);
706 3 : d(2, 2) = (*this)(0, 0) * (*this)(1, 1) - (*this)(1, 0) * (*this)(0, 1);
707 :
708 3 : return d;
709 0 : }
710 :
711 : template <typename T>
712 : void
713 0 : RankTwoTensorTempl<T>::print(std::ostream & stm) const
714 : {
715 0 : const RankTwoTensorTempl<T> & a = *this;
716 0 : for (const auto i : make_range(N))
717 : {
718 0 : for (const auto j : make_range(N))
719 0 : stm << std::setw(15) << a(i, j) << ' ';
720 0 : stm << std::endl;
721 : }
722 0 : }
723 :
724 : template <>
725 : void RankTwoTensor::printReal(std::ostream & stm) const;
726 :
727 : template <>
728 : void ADRankTwoTensor::printReal(std::ostream & stm) const;
729 :
730 : template <>
731 : void ADRankTwoTensor::printADReal(unsigned int nDual, std::ostream & stm) const;
732 :
733 : template <typename T>
734 : void
735 1327 : RankTwoTensorTempl<T>::addIa(const T & a)
736 : {
737 5308 : for (const auto i : make_range(N))
738 3981 : (*this)(i, i) += a;
739 1327 : }
740 :
741 : template <typename T>
742 : T
743 9432137 : RankTwoTensorTempl<T>::L2norm() const
744 : {
745 9432137 : T norm = 0.0;
746 94321370 : for (const auto i : make_range(N2))
747 : {
748 84889233 : T v = _coords[i];
749 84889233 : norm += v * v;
750 : }
751 9432138 : return norm == 0.0 ? 0.0 : std::sqrt(norm);
752 1 : }
753 :
754 : template <typename T>
755 : void
756 0 : RankTwoTensorTempl<T>::surfaceFillFromInputVector(const std::vector<T> & input)
757 : {
758 0 : if (input.size() == 4)
759 : {
760 : // initialize with zeros
761 0 : this->zero();
762 0 : (*this)(0, 0) = input[0];
763 0 : (*this)(0, 1) = input[1];
764 0 : (*this)(1, 0) = input[2];
765 0 : (*this)(1, 1) = input[3];
766 : }
767 : else
768 0 : mooseError("please provide correct number of values for surface RankTwoTensorTempl<T> "
769 : "initialization.");
770 0 : }
771 :
772 : template <typename T>
773 : void
774 0 : RankTwoTensorTempl<T>::syev(const char *, std::vector<T> &, std::vector<T> &) const
775 : {
776 0 : mooseError("The syev method is only supported for Real valued tensors");
777 : }
778 :
779 : template <>
780 : void RankTwoTensor::syev(const char * calculation_type,
781 : std::vector<Real> & eigvals,
782 : std::vector<Real> & a) const;
783 :
784 : template <typename T>
785 : void
786 0 : RankTwoTensorTempl<T>::symmetricEigenvalues(std::vector<T> & eigvals) const
787 : {
788 0 : RankTwoTensorTempl<T> a;
789 0 : symmetricEigenvaluesEigenvectors(eigvals, a);
790 0 : }
791 :
792 : template <>
793 : void RankTwoTensor::symmetricEigenvalues(std::vector<Real> & eigvals) const;
794 :
795 : template <typename T>
796 : void
797 : RankTwoTensorTempl<T>::symmetricEigenvaluesEigenvectors(std::vector<T> &,
798 : RankTwoTensorTempl<T> &) const
799 : {
800 : mooseError(
801 : "symmetricEigenvaluesEigenvectors is only available for ordered tensor component types");
802 : }
803 :
804 : template <>
805 : void RankTwoTensor::symmetricEigenvaluesEigenvectors(std::vector<Real> & eigvals,
806 : RankTwoTensor & eigvecs) const;
807 :
808 : template <>
809 : void ADRankTwoTensor::symmetricEigenvaluesEigenvectors(std::vector<ADReal> & eigvals,
810 : ADRankTwoTensor & eigvecs) const;
811 :
812 : template <typename T>
813 : void
814 494 : RankTwoTensorTempl<T>::dsymmetricEigenvalues(std::vector<T> & eigvals,
815 : std::vector<RankTwoTensorTempl<T>> & deigvals) const
816 : {
817 494 : deigvals.resize(N);
818 :
819 494 : std::vector<T> a;
820 494 : syev("V", eigvals, a);
821 :
822 : // now a contains the eigenvetors
823 : // extract these and place appropriately in deigvals
824 494 : std::vector<T> eig_vec;
825 494 : eig_vec.resize(N);
826 :
827 1976 : for (const auto i : make_range(N))
828 : {
829 5928 : for (const auto j : make_range(N))
830 4446 : eig_vec[j] = a[i * N + j];
831 5928 : for (const auto j : make_range(N))
832 17784 : for (const auto k : make_range(N))
833 13338 : deigvals[i](j, k) = eig_vec[j] * eig_vec[k];
834 : }
835 :
836 : // There are discontinuities in the derivative
837 : // for equal eigenvalues. The following is
838 : // an attempt to make a sensible choice for
839 : // the derivative. This agrees with a central-difference
840 : // approximation to the derivative.
841 494 : if (eigvals[0] == eigvals[1] && eigvals[0] == eigvals[2])
842 0 : deigvals[0] = deigvals[1] = deigvals[2] = (deigvals[0] + deigvals[1] + deigvals[2]) / 3.0;
843 494 : else if (eigvals[0] == eigvals[1])
844 2 : deigvals[0] = deigvals[1] = (deigvals[0] + deigvals[1]) / 2.0;
845 492 : else if (eigvals[0] == eigvals[2])
846 0 : deigvals[0] = deigvals[2] = (deigvals[0] + deigvals[2]) / 2.0;
847 492 : else if (eigvals[1] == eigvals[2])
848 2 : deigvals[1] = deigvals[2] = (deigvals[1] + deigvals[2]) / 2.0;
849 494 : }
850 :
851 : template <typename T>
852 : void
853 4 : RankTwoTensorTempl<T>::d2symmetricEigenvalues(std::vector<RankFourTensorTempl<T>> & deriv) const
854 : {
855 4 : std::vector<T> eigvec;
856 4 : std::vector<T> eigvals;
857 0 : T ev[N][N];
858 :
859 : // reset rank four tensor
860 4 : deriv.assign(N, RankFourTensorTempl<T>());
861 :
862 : // get eigen values and eigen vectors
863 4 : syev("V", eigvals, eigvec);
864 :
865 16 : for (const auto i : make_range(N))
866 48 : for (const auto j : make_range(N))
867 36 : ev[i][j] = eigvec[i * N + j];
868 :
869 16 : for (unsigned int alpha = 0; alpha < N; ++alpha)
870 48 : for (unsigned int beta = 0; beta < N; ++beta)
871 : {
872 36 : if (eigvals[alpha] == eigvals[beta])
873 12 : continue;
874 :
875 96 : for (const auto i : make_range(N))
876 288 : for (const auto j : make_range(N))
877 864 : for (const auto k : make_range(N))
878 2592 : for (const auto l : make_range(N))
879 : {
880 1944 : deriv[alpha](i, j, k, l) +=
881 1944 : 0.5 * (ev[beta][i] * ev[alpha][j] + ev[alpha][i] * ev[beta][j]) *
882 3888 : (ev[beta][k] * ev[alpha][l] + ev[beta][l] * ev[alpha][k]) /
883 1944 : (eigvals[alpha] - eigvals[beta]);
884 : }
885 : }
886 4 : }
887 :
888 : template <typename T>
889 : void
890 0 : RankTwoTensorTempl<T>::getRUDecompositionRotation(RankTwoTensorTempl<T> &) const
891 : {
892 0 : mooseError("getRUDecompositionRotation is only supported for Real valued tensors");
893 : }
894 :
895 : template <>
896 : void RankTwoTensor::getRUDecompositionRotation(RankTwoTensor & rot) const;
897 :
898 : template <typename T>
899 : void
900 1 : RankTwoTensorTempl<T>::initRandom(unsigned int rand_seed)
901 : {
902 1 : MooseRandom::seed(rand_seed);
903 1 : }
904 :
905 : template <typename T>
906 : RankTwoTensorTempl<T>
907 0 : RankTwoTensorTempl<T>::genRandomTensor(T stddev, T mean)
908 : {
909 0 : RankTwoTensorTempl<T> result;
910 0 : for (const auto i : make_range(N))
911 0 : for (const auto j : make_range(N))
912 0 : result(i, j) = (MooseRandom::rand() + mean) * stddev;
913 0 : return result;
914 0 : }
915 :
916 : template <typename T>
917 : RankTwoTensorTempl<T>
918 1 : RankTwoTensorTempl<T>::genRandomSymmTensor(T stddev, T mean)
919 : {
920 7 : auto r = [&]() { return (MooseRandom::rand() + mean) * stddev; };
921 2 : return RankTwoTensorTempl<T>(r(), r(), r(), r(), r(), r());
922 : }
923 :
924 : template <typename T>
925 : void
926 0 : RankTwoTensorTempl<T>::vectorOuterProduct(const libMesh::TypeVector<T> & v1,
927 : const libMesh::TypeVector<T> & v2)
928 : {
929 0 : RankTwoTensorTempl<T> & a = *this;
930 0 : for (const auto i : make_range(N))
931 0 : for (const auto j : make_range(N))
932 0 : a(i, j) = v1(i) * v2(j);
933 0 : }
934 :
935 : template <typename T>
936 : RankTwoTensorTempl<T>
937 3 : RankTwoTensorTempl<T>::outerProduct(const libMesh::TypeVector<T> & v1,
938 : const libMesh::TypeVector<T> & v2)
939 : {
940 3 : RankTwoTensorTempl<T> result;
941 12 : for (const auto i : make_range(N))
942 36 : for (const auto j : make_range(N))
943 27 : result(i, j) = v1(i) * v2(j);
944 3 : return result;
945 0 : }
946 :
947 : template <typename T>
948 : RankTwoTensorTempl<T>
949 608 : RankTwoTensorTempl<T>::selfOuterProduct(const libMesh::TypeVector<T> & v)
950 : {
951 608 : RankTwoTensorTempl<T> result(RankTwoTensorTempl<T>::initNone);
952 2432 : for (unsigned int i = 0; i < N; ++i)
953 7296 : for (unsigned int j = 0; j < N; ++j)
954 5472 : result(i, j) = v(i) * v(j);
955 608 : return result;
956 0 : }
957 :
958 : template <typename T>
959 : void
960 1 : RankTwoTensorTempl<T>::fillRealTensor(libMesh::TensorValue<T> & tensor)
961 : {
962 4 : for (const auto i : make_range(N))
963 12 : for (const auto j : make_range(N))
964 9 : tensor(i, j) = (*this)(i, j);
965 1 : }
966 :
967 : template <typename T>
968 : void
969 1 : RankTwoTensorTempl<T>::fillRow(unsigned int r, const libMesh::TypeVector<T> & v)
970 : {
971 4 : for (const auto i : make_range(N))
972 3 : (*this)(r, i) = v(i);
973 1 : }
974 :
975 : template <typename T>
976 : void
977 1 : RankTwoTensorTempl<T>::fillColumn(unsigned int c, const libMesh::TypeVector<T> & v)
978 : {
979 4 : for (const auto i : make_range(N))
980 3 : (*this)(i, c) = v(i);
981 1 : }
982 :
983 : template <typename T>
984 : RankTwoTensorTempl<T>
985 2 : RankTwoTensorTempl<T>::initialContraction(const RankFourTensorTempl<T> & b) const
986 : {
987 2 : RankTwoTensorTempl<T> result;
988 8 : for (const auto i : make_range(N))
989 24 : for (const auto j : make_range(N))
990 72 : for (const auto k : make_range(N))
991 216 : for (const auto l : make_range(N))
992 162 : result(k, l) += (*this)(i, j) * b(i, j, k, l);
993 2 : return result;
994 0 : }
995 :
996 : template <typename T>
997 : void
998 0 : RankTwoTensorTempl<T>::setToIdentity()
999 : {
1000 : mooseAssert(N2 == 9, "RankTwoTensorTempl is currently only tested for 3 dimensions.");
1001 0 : for (const auto i : make_range(N2))
1002 0 : _coords[i] = identityCoords[i];
1003 0 : }
|
