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 "Moose.h"
13 : #include "MooseTypes.h"
14 : #include "ADRankTwoTensorForward.h"
15 : #include "ADRankFourTensorForward.h"
16 : #include "ADRankThreeTensorForward.h"
17 : #include "ADSymmetricRankTwoTensorForward.h"
18 : #include "MooseUtils.h"
19 : #include "MathUtils.h"
20 :
21 : // Any requisite includes here
22 : #include "libmesh/libmesh.h"
23 : #include "libmesh/tensor_value.h"
24 : #include "libmesh/vector_value.h"
25 : #include "libmesh/int_range.h"
26 :
27 : #include "metaphysicl/raw_type.h"
28 :
29 : #include <petscsys.h>
30 : #include <vector>
31 :
32 : // Forward declarations
33 : class MooseEnum;
34 : template <typename T>
35 : class MooseArray;
36 : typedef MooseArray<Real> VariableValue;
37 : template <typename>
38 : class ColumnMajorMatrixTempl;
39 : namespace libMesh
40 : {
41 : template <typename>
42 : class TypeVector;
43 : template <typename>
44 : class TypeTensor;
45 : template <typename>
46 : class TensorValue;
47 : namespace TensorTools
48 : {
49 : template <>
50 : struct IncrementRank<RankTwoTensor>
51 : {
52 : typedef RankThreeTensor type;
53 : };
54 : }
55 : }
56 :
57 : namespace MathUtils
58 : {
59 : template <typename T>
60 : void mooseSetToZero(T & v);
61 :
62 : /**
63 : * Helper function template specialization to set an object to zero.
64 : * Needed by DerivativeMaterialInterface
65 : */
66 : template <>
67 : void mooseSetToZero<RankTwoTensor>(RankTwoTensor & v);
68 :
69 : /**
70 : * Helper function template specialization to set an object to zero.
71 : * Needed by DerivativeMaterialInterface
72 : */
73 : template <>
74 : void mooseSetToZero<ADRankTwoTensor>(ADRankTwoTensor & v);
75 : }
76 :
77 : /**
78 : * RankTwoTensorTempl is designed to handle the Stress or Strain Tensor for a fully anisotropic
79 : * material. It is designed to allow for maximum clarity of the mathematics and ease of use.
80 : * Original class authors: A. M. Jokisaari, O. Heinonen, M. R. Tonks
81 : *
82 : * RankTwoTensorTempl holds the 9 separate Sigma_ij or Epsilon_ij entries.
83 : * The entries are accessed by index, with i, j equal to 1, 2, or 3, or
84 : * internally i, j = 0, 1, 2.
85 : */
86 : template <typename T>
87 : class RankTwoTensorTempl : public libMesh::TensorValue<T>
88 : {
89 : public:
90 : /// @{
91 :
92 : /**
93 : * @brief Tensor dimension, i.e. number of rows/columns of the second order tensor
94 : */
95 : static constexpr unsigned int N = Moose::dim;
96 :
97 : /**
98 : * @brief The square of the tensor dimension.
99 : */
100 : static constexpr unsigned int N2 = N * N;
101 :
102 : /**
103 : * @brief The initialization method.
104 : * @see fillFromInputVector()
105 : * @see RankTwoTensorTempl(const std::vector<T> &)
106 : */
107 : enum InitMethod
108 : {
109 : initNone,
110 : initIdentity
111 : };
112 :
113 : /**
114 : * To fill up the 9 entries in the 2nd-order tensor, fillFromInputVector
115 : * is called with one of the following fill_methods.
116 : * See the fill*FromInputVector functions for more details
117 : */
118 : enum FillMethod
119 : {
120 : autodetect = 0,
121 : isotropic1 = 1,
122 : diagonal3 = 3,
123 : symmetric6 = 6,
124 : general = 9
125 : };
126 :
127 : /**
128 : * @brief Initialize the random seed based on an unsigned integer.
129 : * Deprecated in favor of MooseRandom::seed().
130 : */
131 : static void initRandom(unsigned int);
132 :
133 : /**
134 : * @brief Get the available `FillMethod` options.
135 : *
136 : * This method is useful in validParams().
137 : */
138 : [[nodiscard]] static MooseEnum fillMethodEnum();
139 :
140 : /**
141 : * @brief Fill a `libMesh::TensorValue<T>` from this second order tensor
142 : */
143 : void fillRealTensor(libMesh::TensorValue<T> &);
144 :
145 : /// Print the rank two tensor
146 : void print(std::ostream & stm = Moose::out) const;
147 :
148 : /// Print the Real part of the RankTwoTensorTempl<ADReal>
149 : void printReal(std::ostream & stm = Moose::out) const;
150 :
151 : /// Print the Real part of the RankTwoTensorTempl<ADReal> along with its first nDual dual numbers
152 : void printADReal(unsigned int nDual, std::ostream & stm = Moose::out) const;
153 :
154 : /// @}
155 :
156 : /// @{
157 :
158 : /**
159 : * @brief Empty constructor; fills to zero
160 : *
161 : * \f$ A_{ij} = 0 \f$
162 : *
163 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
164 : * RankTwoTensor A;
165 : * // A = [ 0 0 0
166 : * // 0 0 0
167 : * // 0 0 0 ]
168 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
169 : */
170 : RankTwoTensorTempl();
171 :
172 : /**
173 : * @brief Select specific initialization pattern.
174 : *
175 : * `initNone` initializes an empty second order tensor.
176 : *
177 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
178 : * RankTwoTensor A(RankTwoTensor::initNone);
179 : * // A = [ 0 0 0
180 : * // 0 0 0
181 : * // 0 0 0 ]
182 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
183 : *
184 : * `initIdentity` initializes a second order identity tensor.
185 : *
186 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
187 : * RankTwoTensor A(RankTwoTensor::initIdentity);
188 : * // A = [ 1 0 0
189 : * // 0 1 0
190 : * // 0 0 1 ]
191 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
192 : */
193 : RankTwoTensorTempl(const InitMethod);
194 :
195 : /**
196 : * @brief Initialize from row vectors.
197 : * @see initializeFromRows()
198 : *
199 : * Deprecated in favor of initializeFromRows()
200 : *
201 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
202 : * RealVectorValue row1(1, 2, 3);
203 : * RealVectorValue row2(4, 5, 6);
204 : * RealVectorValue row3(7, 8, 9);
205 : * RankTwoTensor A(row1, row2, row3);
206 : * // A = [ 1 2 3
207 : * // 4 5 6
208 : * // 7 8 9 ]
209 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
210 : */
211 : RankTwoTensorTempl(const libMesh::TypeVector<T> & row1,
212 : const libMesh::TypeVector<T> & row2,
213 : const libMesh::TypeVector<T> & row3);
214 :
215 : /**
216 : * @brief Constructor that proxies the fillFromInputVector() method
217 : * @see fillFromInputVector()
218 : */
219 9 : RankTwoTensorTempl(const std::vector<T> & input) { this->fillFromInputVector(input); };
220 :
221 : /**
222 : * @brief Initialize a symmetric second order tensor using the 6 arguments.
223 : *
224 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
225 : * RankTwoTensor A(1, 2, 3, 4, 5, 6);
226 : * // A = [ 1 6 5
227 : * // 6 2 4
228 : * // 5 4 3 ]
229 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
230 : */
231 : RankTwoTensorTempl(
232 : const T & S11, const T & S22, const T & S33, const T & S23, const T & S13, const T & S12);
233 :
234 : /**
235 : * @brief Initialize a second order tensor using the 9 arguments in a column-major fashion.
236 : *
237 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
238 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
239 : * // A = [ 1 4 7
240 : * // 2 5 8
241 : * // 3 6 9 ]
242 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
243 : */
244 : RankTwoTensorTempl(const T & S11,
245 : const T & S21,
246 : const T & S31,
247 : const T & S12,
248 : const T & S22,
249 : const T & S32,
250 : const T & S13,
251 : const T & S23,
252 : const T & S33);
253 :
254 : /**
255 : * @brief The copy constructor
256 : *
257 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
258 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
259 : * RankTwoTensor B = A;
260 : * // A = B = [ 1 4 7
261 : * // 2 5 8
262 : * // 3 6 9 ]
263 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
264 : */
265 22 : RankTwoTensorTempl(const RankTwoTensorTempl<T> & a) = default;
266 :
267 : /**
268 : * @brief The conversion operator from a `libMesh::TensorValue`
269 : */
270 264056 : RankTwoTensorTempl(const libMesh::TensorValue<T> & a) : libMesh::TensorValue<T>(a) {}
271 :
272 : /**
273 : * @brief The conversion operator from a `libMesh::TypeTensor`
274 : */
275 797213 : RankTwoTensorTempl(const libMesh::TypeTensor<T> & a) : libMesh::TensorValue<T>(a) {}
276 :
277 : /**
278 : * @brief The conversion operator from a `SymmetricRankTwoTensorTempl`
279 : */
280 : template <typename T2>
281 2 : RankTwoTensorTempl(const SymmetricRankTwoTensorTempl<T2> & a)
282 0 : : RankTwoTensorTempl<T>(a(0),
283 2 : a(1),
284 2 : a(2),
285 2 : a(3) / MathUtils::sqrt2,
286 2 : a(4) / MathUtils::sqrt2,
287 4 : a(5) / MathUtils::sqrt2)
288 : {
289 2 : }
290 :
291 : /**
292 : * @brief The conversion operator from `RankTwoTensorTempl<T2>` to `RankTwoTensorTempl<T>` where
293 : * `T2` is convertible to `T`.
294 : */
295 : template <typename T2>
296 36 : RankTwoTensorTempl(const RankTwoTensorTempl<T2> & a) : libMesh::TensorValue<T>(a)
297 : {
298 36 : }
299 : /// @} end Constructor
300 :
301 : /// @{
302 :
303 : /**
304 : * @brief Named constructor for initializing symmetrically.
305 : *
306 : * The supplied vectors are used as row and column vectors to construct two tensors respectively,
307 : * that are averaged to create a symmetric tensor.
308 : *
309 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
310 : * RealVectorValue row1(1, 2, 3);
311 : * RealVectorValue row2(4, 5, 6);
312 : * RealVectorValue row3(7, 8, 9);
313 : * auto A = RankTwoTensor::initializeSymmetric(row1, row2, row3);
314 : * // A = [ 1 3 5
315 : * // 3 5 7
316 : * // 5 7 9 ]
317 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
318 : */
319 : [[nodiscard]] static RankTwoTensorTempl initializeSymmetric(const libMesh::TypeVector<T> & v0,
320 : const libMesh::TypeVector<T> & v1,
321 : const libMesh::TypeVector<T> & v2);
322 :
323 : /**
324 : * @brief Named constructor for initializing from row vectors.
325 : *
326 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
327 : * RealVectorValue row1(1, 2, 3);
328 : * RealVectorValue row2(4, 5, 6);
329 : * RealVectorValue row3(7, 8, 9);
330 : * auto A = RankTwoTensor::initializeFromRows(row1, row2, row3);
331 : * // A = [ 1 2 3
332 : * // 4 5 6
333 : * // 7 8 9 ]
334 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
335 : */
336 : [[nodiscard]] static RankTwoTensorTempl initializeFromRows(const libMesh::TypeVector<T> & row0,
337 : const libMesh::TypeVector<T> & row1,
338 : const libMesh::TypeVector<T> & row2);
339 :
340 : /**
341 : * @brief Named constructor for initializing from row vectors.
342 : *
343 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
344 : * RealVectorValue col1(1, 2, 3);
345 : * RealVectorValue col2(4, 5, 6);
346 : * RealVectorValue col3(7, 8, 9);
347 : * auto A = RankTwoTensor::initializeFromColumns(col1, col2, col3);
348 : * // A = [ 1 4 7
349 : * // 2 5 8
350 : * // 3 6 9 ]
351 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
352 : */
353 : [[nodiscard]] static RankTwoTensorTempl
354 : initializeFromColumns(const libMesh::TypeVector<T> & col0,
355 : const libMesh::TypeVector<T> & col1,
356 : const libMesh::TypeVector<T> & col2);
357 :
358 : /**
359 : * @brief Initialize a second order identity tensor.
360 : *
361 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
362 : * auto A = RankTwoTensor::Identity();
363 : * // A = [ 1 0 0
364 : * // 0 1 0
365 : * // 0 0 1 ]
366 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
367 : */
368 0 : [[nodiscard]] static RankTwoTensorTempl Identity() { return RankTwoTensorTempl(initIdentity); }
369 :
370 : /**
371 : * @brief Initialize a second order tensor with expression \f$ B_{ij} = F_{ij} F_{ji} \f$.
372 : *
373 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
374 : * RankTwoTensor F(1, 2, 3, 4, 5, 6, 7, 8, 9);
375 : * // F = [ 1 4 7
376 : * // 2 5 8
377 : * // 3 6 9 ]
378 : * auto B = RankTwoTensor::timesTranspose(F);
379 : * // B = [ 66 78 90
380 : * // 78 93 108
381 : * // 90 108 126 ]
382 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
383 : */
384 : [[nodiscard]] static RankTwoTensorTempl<T> timesTranspose(const RankTwoTensorTempl<T> &);
385 :
386 : /**
387 : * @brief Initialize a second order tensor with expression \f$ C_{ij} = F_{ji} F_{ij} \f$.
388 : *
389 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
390 : * RankTwoTensor F(1, 2, 3, 4, 5, 6, 7, 8, 9);
391 : * // F = [ 1 4 7
392 : * // 2 5 8
393 : * // 3 6 9 ]
394 : * auto C = RankTwoTensor::transposeTimes(F);
395 : * // C = [ 14 32 50
396 : * // 32 77 122
397 : * // 50 122 194 ]
398 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
399 : */
400 : [[nodiscard]] static RankTwoTensorTempl<T> transposeTimes(const RankTwoTensorTempl<T> &);
401 :
402 : /**
403 : * @brief Initialize a second order tensor with expression \f$ E_{ij} = C_{ij} + C_{ji} \f$.
404 : *
405 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
406 : * RankTwoTensor C(1, 2, 3, 4, 5, 6, 7, 8, 9);
407 : * // C = [ 1 4 7
408 : * // 2 5 8
409 : * // 3 6 9 ]
410 : * auto E = RankTwoTensor::plusTranspose(C);
411 : * // E = [ 2 6 10
412 : * // 6 10 14
413 : * // 10 14 18 ]
414 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
415 : */
416 : [[nodiscard]] static RankTwoTensorTempl<T> plusTranspose(const RankTwoTensorTempl<T> &);
417 :
418 : /**
419 : * @brief Initialize a second order tensor as the outer product of two vectors, i.e. \f$ A_{ij} =
420 : * a_i b_j \f$.
421 : *
422 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
423 : * RealVectorValue a(1, 2, 3);
424 : * RealVectorValue b(4, 5, 6);
425 : * auto A = RankTwoTensor::outerProduct(a, b);
426 : * // A = [ 4 5 6
427 : * // 8 10 12
428 : * // 12 15 18 ]
429 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
430 : */
431 : [[nodiscard]] static RankTwoTensorTempl<T> outerProduct(const libMesh::TypeVector<T> &,
432 : const libMesh::TypeVector<T> &);
433 :
434 : /**
435 : * @brief Initialize a second order tensor as the outer product of a vector with itself, i.e. \f$
436 : * A_{ij} = a_i a_j \f$.
437 : *
438 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
439 : * RealVectorValue a(1, 2, 3);
440 : * auto A = RankTwoTensor::selfOuterProduct(a);
441 : * // A = [ 1 2 3
442 : * // 2 4 6
443 : * // 3 6 9 ]
444 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
445 : */
446 : [[nodiscard]] static RankTwoTensorTempl<T> selfOuterProduct(const libMesh::TypeVector<T> &);
447 :
448 : /**
449 : * @brief Generate a random second order tensor with all 9 components treated as independent
450 : * random variables following a Gaussian distribution.
451 : *
452 : * @param stddev The standard deviation of the Gaussian distribution
453 : * @param mean The mean of the Gaussian distribution
454 : */
455 : [[nodiscard]] static RankTwoTensorTempl<T> genRandomTensor(T stddev, T mean);
456 :
457 : /**
458 : * @brief Generate a random symmetric second order tensor with the 6 upper-triangular components
459 : * treated as independent random variables following a Gaussian distribution.
460 : *
461 : * @param stddev The standard deviation of the Gaussian distribution
462 : * @param mean The mean of the Gaussian distribution
463 : */
464 : [[nodiscard]] static RankTwoTensorTempl<T> genRandomSymmTensor(T stddev, T mean);
465 :
466 : /// @}
467 :
468 : /// @{
469 :
470 : /**
471 : * @brief Get the i-th column of the second order tensor.
472 : *
473 : * @param i The column number, i = 0, 1, 2
474 : * @return The i-th column of the second order tensor.
475 : *
476 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
477 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
478 : * RealVectorValue col = A.column(1);
479 : * // col = [ 4
480 : * // 5
481 : * // 6 ]
482 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
483 : */
484 : libMesh::VectorValue<T> column(const unsigned int i) const;
485 :
486 : /// @}
487 :
488 : /// @{
489 :
490 : /**
491 : * @brief Assignment operator
492 : *
493 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
494 : * RankTwoTensor A;
495 : * // A = [ 0 0 0
496 : * // 0 0 0
497 : * // 0 0 0 ]
498 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
499 : * // B = [ 9 6 3
500 : * // 8 5 2
501 : * // 7 4 1 ]
502 : * A = B;
503 : * // A = [ 9 6 3
504 : * // 8 5 2
505 : * // 7 4 1 ]
506 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
507 : */
508 96449 : RankTwoTensorTempl<T> & operator=(const RankTwoTensorTempl<T> & a) = default;
509 :
510 : /**
511 : * @brief Assignment operator (from a ColumnMajorMatrixTempl<T>)
512 : *
513 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
514 : * RankTwoTensor A;
515 : * // A = [ 0 0 0
516 : * // 0 0 0
517 : * // 0 0 0 ]
518 : * RealVectorValue col1(1, 2, 3);
519 : * RealVectorValue col2(4, 5, 6);
520 : * RealVectorValue col3(7, 8, 9);
521 : * ColumnMajorMatrix B(col1, col2, col3);
522 : * // B = [ 1 4 7
523 : * // 2 5 8
524 : * // 3 6 9 ]
525 : * A = B;
526 : * // A = [ 1 4 7
527 : * // 2 5 8
528 : * // 3 6 9 ]
529 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
530 : */
531 : RankTwoTensorTempl<T> & operator=(const ColumnMajorMatrixTempl<T> & a);
532 :
533 : /**
534 : * @brief Assignment-from-scalar operator. Used only to zero out the tensor.
535 : *
536 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
537 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
538 : * // A = [ 1 2 3
539 : * // 2 4 6
540 : * // 3 6 9 ]
541 : * A = 0;
542 : * // A = [ 0 0 0
543 : * // 0 0 0
544 : * // 0 0 0 ]
545 : * A = 1;
546 : * // This triggers an assertion failure.
547 : * A = 0.0;
548 : * // This triggers an assertion failure.
549 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
550 : */
551 : template <typename Scalar>
552 : typename libMesh::boostcopy::enable_if_c<libMesh::ScalarTraits<Scalar>::value,
553 : RankTwoTensorTempl &>::type
554 : operator=(const Scalar & libmesh_dbg_var(p))
555 : {
556 : libmesh_assert_equal_to(p, Scalar(0));
557 : this->zero();
558 : return *this;
559 : }
560 :
561 : /**
562 : * Add another second order tensor to this one
563 : *
564 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
565 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
566 : * // A = [ 1 2 3
567 : * // 2 4 6
568 : * // 3 6 9 ]
569 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
570 : * // B = [ 9 6 3
571 : * // 8 5 2
572 : * // 7 4 1 ]
573 : * A += B;
574 : * // A = [ 10 8 6
575 : * // 10 9 8
576 : * // 10 10 10 ]
577 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
578 : */
579 : RankTwoTensorTempl<T> & operator+=(const RankTwoTensorTempl<T> & a);
580 :
581 : /**
582 : * Subtract another second order tensor from this one
583 : *
584 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
585 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
586 : * // A = [ 1 2 3
587 : * // 2 4 6
588 : * // 3 6 9 ]
589 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
590 : * // B = [ 9 6 3
591 : * // 8 5 2
592 : * // 7 4 1 ]
593 : * A -= B;
594 : * // A = [ -8 -4 0
595 : * // -6 -1 4
596 : * // -4 2 8 ]
597 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
598 : */
599 : RankTwoTensorTempl<T> & operator-=(const RankTwoTensorTempl<T> & a);
600 :
601 : /**
602 : * Multiply this tensor by a scalar (component-wise)
603 : *
604 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
605 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
606 : * // A = [ 1 2 3
607 : * // 2 4 6
608 : * // 3 6 9 ]
609 : * A *= 2;
610 : * // A = [ 2 4 6
611 : * // 4 8 12
612 : * // 6 12 18 ]
613 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
614 : */
615 : RankTwoTensorTempl<T> & operator*=(const T & a);
616 :
617 : /**
618 : * Divide this tensor by a scalar (component-wise)
619 : *
620 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
621 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
622 : * // A = [ 1 2 3
623 : * // 2 4 6
624 : * // 3 6 9 ]
625 : * A /= 2;
626 : * // A = [ 0.5 1.0 1.5
627 : * // 1.0 2.0 3.0
628 : * // 1.5 3.0 4.5 ]
629 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
630 : */
631 : RankTwoTensorTempl<T> & operator/=(const T & a);
632 :
633 : /**
634 : * Multiplication with another second order tensor
635 : *
636 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
637 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
638 : * // A = [ 1 2 3
639 : * // 2 4 6
640 : * // 3 6 9 ]
641 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
642 : * // B = [ 9 6 3
643 : * // 8 5 2
644 : * // 7 4 1 ]
645 : * A *= B;
646 : * // A = [ 90 54 18
647 : * // 114 69 24
648 : * // 138 84 30 ]
649 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
650 : */
651 : RankTwoTensorTempl<T> & operator*=(const libMesh::TypeTensor<T> & a);
652 :
653 : /**
654 : * @brief The smart mutator that determines how to fill the second order tensor based on the size
655 : * of the input vector.
656 : *
657 : * @param input The input vector, can be of size 1, 3, 6, or 9
658 : * @param fill_method The fill method, default to autodetect.
659 : *
660 : * When `input.size() == 1`, the vector value is used to fill the diagonal components of the
661 : * second order tensor:
662 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
663 : * RankTwoTensor A;
664 : * A.fillFromInputVector({1.5});
665 : * // A = [ 1.5 0 0
666 : * // 0 1.5 0
667 : * // 0 0 1.5 ]
668 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
669 : *
670 : * When `input.size() == 3`, the vector values are used to fill the diagonal components of the
671 : * second order tensor:
672 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
673 : * RankTwoTensor A;
674 : * A.fillFromInputVector({1, 2, 3});
675 : * // A = [ 1 0 0
676 : * // 0 2 0
677 : * // 0 0 3 ]
678 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
679 : *
680 : * When `input.size() == 6`, the second order tensor is filled symmetrically using the Voigt
681 : * notation:
682 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
683 : * RankTwoTensor A;
684 : * A.fillFromInputVector({1, 2, 3, 4, 5, 6});
685 : * // A = [ 1 6 5
686 : * // 6 2 4
687 : * // 5 4 3 ]
688 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
689 : *
690 : * When `input.size() == 9`, all components of the second order tensor are filled in a
691 : * column-major fashion:
692 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
693 : * RankTwoTensor A;
694 : * A.fillFromInputVector({1, 2, 3, 4, 5, 6, 7, 8, 9});
695 : * // A = [ 1 4 7
696 : * // 2 5 8
697 : * // 3 6 9 ]
698 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
699 : */
700 : void fillFromInputVector(const std::vector<T> & input, FillMethod fill_method = autodetect);
701 :
702 : /**
703 : * @brief The smart mutator that determines how to fill the second order tensor based on the order
704 : * of the scalar_variable.
705 : *
706 : * @param scalar_variable The input scalar variable. Supported orders are FIRST, THIRD, and SIXTH.
707 : *
708 : * When `scalar_variable.size() == 1`, the scalar value is used to fill the very first component
709 : * of the second order tensor:
710 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
711 : & // Suppose v[0] = 1
712 : * RankTwoTensor A;
713 : * A.fillFromScalarVariable(v);
714 : * // A = [ 1 0 0
715 : * // 0 0 0
716 : * // 0 0 0 ]
717 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
718 : *
719 : * When `scalar_variable.size() == 3`, the scalar values are used to fill the in-plane components
720 : * of the second order tensor using the Voigt notation:
721 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
722 : * // Suppose v[0] = 1
723 : * // v[1] = 2
724 : * // v[2] = 3
725 : * RankTwoTensor A;
726 : * A.fillFromScalarVariable(v);
727 : * // A = [ 1 3 0
728 : * // 3 2 0
729 : * // 0 0 0 ]
730 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
731 : *
732 : * When `scalar_variable.size() == 6`, the second order tensor is filled symmetrically using the
733 : * Voigt notation:
734 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
735 : * // Suppose v[0] = 1
736 : * // v[1] = 2
737 : * // v[2] = 3
738 : * // v[3] = 4
739 : * // v[4] = 5
740 : * // v[5] = 6
741 : * RankTwoTensor A;
742 : * A.fillFromScalarVariable(v);
743 : * // A = [ 1 6 5
744 : * // 6 2 4
745 : * // 5 4 3 ]
746 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
747 : */
748 : void fillFromScalarVariable(const VariableValue & scalar_variable);
749 :
750 : /**
751 : * sets _coords[0][0], _coords[0][1], _coords[1][0], _coords[1][1] to input,
752 : * and the remainder to zero
753 : */
754 : void surfaceFillFromInputVector(const std::vector<T> & input);
755 :
756 : /**
757 : * @brief Set the values of the second order tensor to be the outer product of two vectors, i.e.
758 : * \f$ A_{ij} = a_i b_j \f$.
759 : *
760 : * Deprecated in favor of outerProduct()
761 : *
762 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
763 : * RealVectorValue a(1, 2, 3);
764 : * RealVectorValue b(4, 5, 6);
765 : * RankTwoTensor A;
766 : * A.vectorOuterProduct(a, b);
767 : * // A = [ 4 5 6
768 : * // 8 10 12
769 : * // 12 15 18 ]
770 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
771 : */
772 : void vectorOuterProduct(const libMesh::TypeVector<T> &, const libMesh::TypeVector<T> &);
773 :
774 : /**
775 : * @brief Assign values to a specific row of the second order tensor.
776 : *
777 : * @param r The row number, r = 0, 1, 2
778 : * @param v The values to be set
779 : *
780 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
781 : * RealVectorValue a(1, 2, 3);
782 : * RankTwoTensor A;
783 : * // A = [ 0 0 0
784 : * // 0 0 0
785 : * // 0 0 0 ]
786 : * A.fillRow(1, a);
787 : * // A = [ 0 0 0
788 : * // 1 2 3
789 : * // 0 0 0 ]
790 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
791 : */
792 : void fillRow(unsigned int r, const libMesh::TypeVector<T> & v);
793 :
794 : /**
795 : * @brief Assign values to a specific column of the second order tensor.
796 : *
797 : * @param c The column number, c = 0, 1, 2
798 : * @param v The values to be set
799 : *
800 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
801 : * RealVectorValue a(1, 2, 3);
802 : * RankTwoTensor A;
803 : * // A = [ 0 0 0
804 : * // 0 0 0
805 : * // 0 0 0 ]
806 : * A.fillColumn(1, a);
807 : * // A = [ 0 1 0
808 : * // 0 2 0
809 : * // 0 3 0 ]
810 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
811 : */
812 : void fillColumn(unsigned int c, const libMesh::TypeVector<T> & v);
813 :
814 : /**
815 : * @brief Set the tensor to identity.
816 : *
817 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
818 : * RankTwoTensor A;
819 : * // A = [ 0 0 0
820 : * // 0 0 0
821 : * // 0 0 0 ]
822 : * A.setToIdentity();
823 : * // A = [ 1 0 0
824 : * // 0 1 0
825 : * // 0 0 1 ]
826 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
827 : */
828 : void setToIdentity();
829 :
830 : /// Add identity times a to _coords
831 : /**
832 : * @brief Add a scalar to diagonal components \f$ A_{ij} + a\delta_{ij} \f$
833 : *
834 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
835 : * RankTwoTensor A;
836 : * // A = [ 0 0 0
837 : * // 0 0 0
838 : * // 0 0 0 ]
839 : * A.addIa(1.5);
840 : * // A = [ 1.5 0 0
841 : * // 0 1.5 0
842 : * // 0 0 1.5 ]
843 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
844 : */
845 : void addIa(const T & a);
846 :
847 : /**
848 : * Rotate the tensor in-place given a rotation tensor
849 : * \f$ A_{ij} \leftarrow R_{ij} A_{jk} R_{jk} \f$
850 : * @param R The rotation tensor
851 : *
852 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
853 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
854 : * // A = [ 1 4 7
855 : * // 2 5 8
856 : * // 3 6 9 ]
857 : * RankTwoTensor R(0, 1, 0, 1, 0, 0, 0, 0, 1);
858 : * // R = [ 0 1 0
859 : * // 1 0 0
860 : * // 0 0 1 ]
861 : * A.rotate(R);
862 : * // A = [ 5 2 8
863 : * // 4 1 7
864 : * // 6 3 9 ]
865 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
866 : */
867 : void rotate(const RankTwoTensorTempl<T> & R);
868 :
869 : /// @}
870 :
871 : /// @{
872 :
873 : /**
874 : * @brief Return \f$ A_{ij} = A_{ik}A_{kj} \f$
875 : *
876 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
877 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
878 : * // A = [ 1 4 7
879 : * // 2 5 8
880 : * // 3 6 9 ]
881 : * RankTwoTensor B = A.square();
882 : * // B = [ 30 66 102
883 : * // 36 81 126
884 : * // 42 96 150 ]
885 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
886 : */
887 : RankTwoTensorTempl<T> square() const;
888 :
889 : /**
890 : * Return the rotated tensor given a rotation tensor
891 : * \f$ A'_{ij} = R_{ij} A_{jk} R_{jk} \f$
892 : * @param R The rotation tensor
893 : *
894 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
895 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
896 : * // A = [ 1 4 7
897 : * // 2 5 8
898 : * // 3 6 9 ]
899 : * RankTwoTensor R(0, 1, 0, 1, 0, 0, 0, 0, 1);
900 : * // R = [ 0 1 0
901 : * // 1 0 0
902 : * // 0 0 1 ]
903 : * RankTwoTensor A_rotated = A.rotated(R);
904 : * // A_rotated = [ 5 2 8
905 : * // 4 1 7
906 : * // 6 3 9 ]
907 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
908 : */
909 : RankTwoTensorTempl<T> rotated(const RankTwoTensorTempl<T> & R) const;
910 :
911 : /**
912 : * Rotate the tensor about the z-axis
913 : * @param a The rotation angle in radians
914 : */
915 : RankTwoTensorTempl<T> rotateXyPlane(T a);
916 :
917 : /**
918 : * Return the tensor transposed
919 : *
920 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
921 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
922 : * // A = [ 1 4 7
923 : * // 2 5 8
924 : * // 3 6 9 ]
925 : * RankTwoTensor At = A.transpose();
926 : * // At = [ 1 2 3
927 : * // 4 5 6
928 : * // 7 8 9 ]
929 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
930 : */
931 : RankTwoTensorTempl<T> transpose() const;
932 :
933 : /**
934 : * Return the sum of two second order tensors
935 : *
936 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
937 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
938 : * // A = [ 1 2 3
939 : * // 2 4 6
940 : * // 3 6 9 ]
941 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
942 : * // B = [ 9 6 3
943 : * // 8 5 2
944 : * // 7 4 1 ]
945 : * RankTwoTensor C = A + B;
946 : * // C = [ 10 8 6
947 : * // 10 9 8
948 : * // 10 10 10 ]
949 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
950 : */
951 : template <typename T2>
952 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
953 : operator+(const libMesh::TypeTensor<T2> & a) const;
954 :
955 : /**
956 : * Return the subtraction of two second order tensors
957 : *
958 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
959 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
960 : * // A = [ 1 2 3
961 : * // 2 4 6
962 : * // 3 6 9 ]
963 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
964 : * // B = [ 9 6 3
965 : * // 8 5 2
966 : * // 7 4 1 ]
967 : * RankTwoTensor C = A - B;
968 : * // C = [ -8 -4 0
969 : * // -6 -1 4
970 : * // -4 2 8 ]
971 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
972 : */
973 : template <typename T2>
974 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
975 : operator-(const libMesh::TypeTensor<T2> & a) const;
976 :
977 : /**
978 : * Return the negation of this tensor
979 : *
980 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
981 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
982 : * // A = [ 1 2 3
983 : * // 2 4 6
984 : * // 3 6 9 ]
985 : * RankTwoTensor B = -A;
986 : * // B = [ -1 -2 -3
987 : * // -2 -4 -6
988 : * // -3 -6 -9 ]
989 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
990 : */
991 : RankTwoTensorTempl<T> operator-() const;
992 :
993 : /**
994 : * Return this tensor multiplied by a scalar (component-wise)
995 : *
996 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
997 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
998 : * // A = [ 1 2 3
999 : * // 2 4 6
1000 : * // 3 6 9 ]
1001 : * RankTwoTensor B = A * 2;
1002 : * // B = [ 2 4 6
1003 : * // 4 8 12
1004 : * // 6 12 18 ]
1005 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1006 : */
1007 : template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type = 0>
1008 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1009 : operator*(const T2 & a) const;
1010 :
1011 : /**
1012 : * Return this tensor divided by a scalar (component-wise)
1013 : *
1014 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1015 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
1016 : * // A = [ 1 2 3
1017 : * // 2 4 6
1018 : * // 3 6 9 ]
1019 : * RankTwoTensor B = A / 2;
1020 : * // B = [ 0.5 1.0 1.5
1021 : * // 1.0 2.0 3.0
1022 : * // 1.5 3.0 4.5 ]
1023 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1024 : */
1025 : template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type = 0>
1026 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1027 : operator/(const T2 & a) const;
1028 :
1029 : /**
1030 : * Return this tensor multiplied by a vector. \f$ b_i = A_{ij} a_j \f$
1031 : *
1032 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1033 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
1034 : * // A = [ 1 2 3
1035 : * // 2 4 6
1036 : * // 3 6 9 ]
1037 : * RealVectorValue a(1, 2, 3);
1038 : * // a = [ 1
1039 : * // 2
1040 : * // 3 ]
1041 : * RealVectorValue b = A * a;
1042 : * // b = [ 30
1043 : * // 36
1044 : * // 42 ]
1045 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1046 : */
1047 : template <typename T2>
1048 : libMesh::TypeVector<typename libMesh::CompareTypes<T, T2>::supertype>
1049 : operator*(const libMesh::TypeVector<T2> & a) const;
1050 :
1051 : /**
1052 : * Multiplication with another second order tensor
1053 : *
1054 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1055 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
1056 : * // A = [ 1 2 3
1057 : * // 2 4 6
1058 : * // 3 6 9 ]
1059 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
1060 : * // B = [ 9 6 3
1061 : * // 8 5 2
1062 : * // 7 4 1 ]
1063 : * RankTwoTensor C = A * B;
1064 : * // C = [ 90 54 18
1065 : * // 114 69 24
1066 : * // 138 84 30 ]
1067 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1068 : */
1069 : template <typename T2>
1070 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1071 : operator*(const libMesh::TypeTensor<T2> & a) const;
1072 :
1073 : /**
1074 : * Return the double contraction with another second order tensor \f$ A_{ij} B_{ij} \f$
1075 : *
1076 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1077 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
1078 : * // A = [ 1 2 3
1079 : * // 2 4 6
1080 : * // 3 6 9 ]
1081 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
1082 : * // B = [ 9 6 3
1083 : * // 8 5 2
1084 : * // 7 4 1 ]
1085 : * Real result = A.doubleContraction(B);
1086 : * // result = 143
1087 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1088 : */
1089 :
1090 : T doubleContraction(const RankTwoTensorTempl<T> & a) const;
1091 :
1092 : /**
1093 : * Return the general tensor product of this second order tensor and another second order tensor
1094 : * defined as \f$ C_{ijkl} = A_{\mathcal{M}(n)\mathcal{M}(o)} B_{\mathcal{M}(p)\mathcal{M}(q)} \f$
1095 : * where the multiplication order is defined by the index map \f$ \mathcal{M}: \{n,o,p,q\} \to
1096 : * \{i,j,k,l\} \f$. The index map is specified using the template parameters. See examples below
1097 : * for detailed explanation.
1098 : *
1099 : * Suppose we have two second order tensors A and B, and we denote the output indices as i, j, k,
1100 : * l:
1101 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1102 : * usingTensorIndices(i, j, k, l);
1103 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1104 : *
1105 : * The outer product of A and B is defined as \f$ A_{ij} B_{kl} \f$, hence the template
1106 : * specialization should be `times<i, j, k, l>`
1107 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1108 : * RankFourTensor C = A.times<i, j, k, l>(B);
1109 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1110 : *
1111 : * The tensor product of A and B, i.e. \f$ A_{ik} B_{jl} \f$ can be expressed using the template
1112 : * specialization `times<i, k, j, l>`
1113 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1114 : * RankFourTensor C = A.times<i, k, j, l>(B);
1115 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1116 : *
1117 : * Similarly, another tensor product of A and B, i.e. \f$ A_{il} B_{jk} \f$ can be expressed using
1118 : * the template specialization `times<i, l, j, k>`
1119 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1120 : * RankFourTensor C = A.times<i, l, j, k>(B);
1121 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1122 : *
1123 : * The combination goes on...
1124 : */
1125 : template <int n, int o, int p, int q>
1126 : RankFourTensorTempl<T> times(const RankTwoTensorTempl<T> & b) const;
1127 :
1128 : /**
1129 : * Return the single contraction of this second order tensor with a fourth order tensor defined as
1130 : * \f$ C_{ijkl} = A_{\mathcal{M}(m)\mathcal{M}(n)}
1131 : * B_{\mathcal{M}(p)\mathcal{M}(q)\mathcal{M}(r)\mathcal{M}(s)} \f$ where the multiplication order
1132 : * is defined by the index map \f$ \mathcal{M}: \{m,n,p,q,r,s\} \to \{i,j,k,l\} \f$. The index map
1133 : * is specified using the template parameters. See examples below for detailed explanation.
1134 : *
1135 : * Suppose we have a second order tensors A and a fourth order tensor B, and we denote the indices
1136 : * (four output indices and a dummy index to be contracted) as i, j, k, l, m:
1137 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1138 : * usingTensorIndices(i, j, k, l, m);
1139 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1140 : *
1141 : * The single contraction of A and B defined as \f$ A_{im} B_{mjkl} \f$ can be expressed using the
1142 : * template specialization `times<i, m, m, j, k, l>`
1143 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1144 : * RankFourTensor C = A.times<i, m, m, j, k, l>(B);
1145 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1146 : *
1147 : * Similarly, another single contraction of A and B, i.e. \f$ A_{m, i} A_{j, k, m, l} \f$ can be
1148 : * expressed using the template specialization `times<m, i, j, k, m, l>`
1149 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1150 : * RankFourTensor C = A.times<m, i, j, k, m, l>(B);
1151 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1152 : *
1153 : * The combination goes on. Note that this method assumes exactly one repeated index.
1154 : */
1155 : template <int n, int o, int p, int q, int r, int s>
1156 : RankFourTensorTempl<T> times(const RankFourTensorTempl<T> & b) const;
1157 :
1158 : /**
1159 : * @brief Return the outer product \f$ C_{ijkl} = A_{ij} B_{kl} \f$
1160 : */
1161 19 : RankFourTensorTempl<T> outerProduct(const RankTwoTensorTempl<T> & b) const
1162 : {
1163 : usingTensorIndices(i_, j_, k_, l_);
1164 19 : return times<i_, j_, k_, l_>(b);
1165 : }
1166 :
1167 : /**
1168 : * @brief Return the single contraction of this second order tensor with a third order tensor \f$
1169 : * C_{ikl} = A_{ij} B_{jkl} \f$
1170 : */
1171 : RankThreeTensorTempl<T> contraction(const RankThreeTensorTempl<T> & b) const;
1172 :
1173 : /**
1174 : * @brief Return the tensor product of this second order tensor with a vector \f$ C_{ijk} = A_{jk}
1175 : * b_{i} \f$
1176 : */
1177 : RankThreeTensorTempl<T> mixedProductJkI(const libMesh::VectorValue<T> & b) const;
1178 :
1179 : /**
1180 : * @brief Return the positive projection tensor
1181 : *
1182 : * Consider the eigenvalue decomposition of this second order tensor \f$ A = V D V^T \f$, the part
1183 : * of this tensor that lies on the positive spectrum is defined as \f$ A_+ = V \left<D\right> V^T
1184 : * \f$ where the angled brackets are the Macaulay brackets. The positive projection tensor is the
1185 : * linear projection from the full spectrum to the positive spectrum, i.e. \f$ A_+ = P A \f$. The
1186 : * derivation of this positive projection tensor can be found in C. Miehe and M. Lambrecht,
1187 : * Commun. Numer. Meth. Engng 2001; 17:337~353
1188 : *
1189 : * @param eigvals The three eigenvalues of this second order tensor will be filled into this
1190 : * vector.
1191 : * @param eigvecs The three eigenvectors of this second order tensor will be filled into this
1192 : * tensor.
1193 : * @return The fourth order positive projection tensor.
1194 : */
1195 : RankFourTensorTempl<T> positiveProjectionEigenDecomposition(std::vector<T> &,
1196 : RankTwoTensorTempl<T> &) const;
1197 :
1198 : /**
1199 : * @brief Return the deviatoric part of this tensor \f$ A_{ij} - \frac{1}{3} A_{kk} \delta_{ij}
1200 : * \f$
1201 : */
1202 : RankTwoTensorTempl<T> deviatoric() const;
1203 :
1204 : /**
1205 : * @brief A wrapper for tr()
1206 : * @see tr()
1207 : */
1208 : T trace() const;
1209 :
1210 : /**
1211 : * Return the derivative of the trace w.r.t. this second order tensor itself \f$ \frac{\partial
1212 : * A_{kk}}{\partial A_{ij}} = \delta_{ij} \f$
1213 : */
1214 : RankTwoTensorTempl<T> dtrace() const;
1215 :
1216 : /**
1217 : * @brief Return the inverse of this second order tensor
1218 : */
1219 : RankTwoTensorTempl<T> inverse() const;
1220 :
1221 : /**
1222 : * @brief Return the principal second invariant of this second order tensor
1223 : *
1224 : * \f$ I_2 = \frac{1}{2} \left( A_{kk}^2 - A_{ij}A_{ij} \right) \f$
1225 : */
1226 : T generalSecondInvariant() const;
1227 :
1228 : /**
1229 : * @brief Return the main second invariant of this second order tensor
1230 : *
1231 : * \f$ J_2 = \frac{1}{2} \left( S_{ij}S_{ij} \right) \f$, where \f$ S_{ij} = A_{ij} -
1232 : * \frac{1}{3}A_{kk}\delta_{ij} \f$
1233 : */
1234 : T secondInvariant() const;
1235 :
1236 : /**
1237 : * Return the derivative of the main second invariant w.r.t. this second order tensor itself \f$
1238 : * \frac{\partial J_2}{\partial A_{ij}} = S_{ij} = A_{ij} -
1239 : * \frac{1}{3}A_{kk}\delta_{ij} \f$
1240 : */
1241 : RankTwoTensorTempl<T> dsecondInvariant() const;
1242 :
1243 : /**
1244 : * Return the second derivative of the main second invariant w.r.t. this second order tensor
1245 : * itself \f$ \frac{\partial^2 J_2}{\partial A_{ij}A_{kl}} = \delta_{ik}\delta_{jl} -
1246 : * \frac{1}{3}\delta_{ij}\delta_{kl} \f$
1247 : */
1248 : RankFourTensorTempl<T> d2secondInvariant() const;
1249 :
1250 : /**
1251 : * @brief Sin(3*Lode_angle)
1252 : *
1253 : * If secondInvariant() <= r0 then return r0_value
1254 : * This is to gaurd against precision-loss errors.
1255 : * Note that sin(3*Lode_angle) is not defined for secondInvariant() = 0
1256 : */
1257 : T sin3Lode(const T & r0, const T & r0_value) const;
1258 :
1259 : /**
1260 : * d(sin3Lode)/dA_ij
1261 : * If secondInvariant() <= r0 then return zero
1262 : * This is to gaurd against precision-loss errors.
1263 : * Note that sin(3*Lode_angle) is not defined for secondInvariant() = 0
1264 : */
1265 : RankTwoTensorTempl<T> dsin3Lode(const T & r0) const;
1266 :
1267 : /**
1268 : * d^2(sin3Lode)/dA_ij/dA_kl
1269 : * If secondInvariant() <= r0 then return zero
1270 : * This is to gaurd against precision-loss errors.
1271 : * Note that sin(3*Lode_angle) is not defined for secondInvariant() = 0
1272 : */
1273 : RankFourTensorTempl<T> d2sin3Lode(const T & r0) const;
1274 :
1275 : /**
1276 : * Denote the _coords[i][j] by A_ij, then
1277 : * S_ij = A_ij - de_ij*tr(A)/3
1278 : * Then this returns det(S + S.transpose())/2
1279 : * Note the explicit symmeterisation
1280 : */
1281 : T thirdInvariant() const;
1282 :
1283 : /**
1284 : * Denote the _coords[i][j] by A_ij, then
1285 : * this returns d(thirdInvariant()/dA_ij
1286 : */
1287 : RankTwoTensorTempl<T> dthirdInvariant() const;
1288 :
1289 : /**
1290 : * Denote the _coords[i][j] by A_ij, then this returns
1291 : * d^2(thirdInvariant)/dA_ij/dA_kl
1292 : */
1293 : RankFourTensorTempl<T> d2thirdInvariant() const;
1294 :
1295 : /**
1296 : * Denote the _coords[i][j] by A_ij, then this returns
1297 : * d(det)/dA_ij
1298 : */
1299 : RankTwoTensorTempl<T> ddet() const;
1300 :
1301 : /// Sqrt(_coords[i][j]*_coords[i][j])
1302 : T L2norm() const;
1303 :
1304 : /**
1305 : * computes eigenvalues, assuming tens is symmetric, and places them
1306 : * in ascending order in eigvals
1307 : */
1308 : void symmetricEigenvalues(std::vector<T> & eigvals) const;
1309 :
1310 : /**
1311 : * computes eigenvalues and eigenvectors, assuming tens is symmetric, and places them
1312 : * in ascending order in eigvals. eigvecs is a matrix with the first column
1313 : * being the first eigenvector, the second column being the second, etc.
1314 : */
1315 : void symmetricEigenvaluesEigenvectors(std::vector<T> & eigvals,
1316 : RankTwoTensorTempl<T> & eigvecs) const;
1317 :
1318 : /**
1319 : * computes eigenvalues, and their symmetric derivatives wrt vals,
1320 : * assuming tens is symmetric
1321 : * @param eigvals are the eigenvalues of the matrix, in ascending order
1322 : * @param deigvals Here digvals[i](j,k) = (1/2)*(d(eigvals[i])/dA_jk + d(eigvals[i]/dA_kj))
1323 : * Note the explicit symmeterisation here.
1324 : * For equal eigenvalues, these derivatives are not gauranteed to
1325 : * be the ones you expect, since the derivatives in this case are
1326 : * often defined by continuation from the un-equal case, and that is
1327 : * too sophisticated for this routine.
1328 : */
1329 : void dsymmetricEigenvalues(std::vector<T> & eigvals,
1330 : std::vector<RankTwoTensorTempl<T>> & deigvals) const;
1331 :
1332 : /**
1333 : * Computes second derivatives of Eigenvalues of a rank two tensor
1334 : * @param deriv store second derivative of the current tensor in here
1335 : */
1336 : void d2symmetricEigenvalues(std::vector<RankFourTensorTempl<T>> & deriv) const;
1337 :
1338 : /**
1339 : * Uses the petscblaslapack.h LAPACKsyev_ routine to perform RU decomposition and obtain the
1340 : * rotation tensor.
1341 : */
1342 : void getRUDecompositionRotation(RankTwoTensorTempl<T> & rot) const;
1343 :
1344 : /// returns this_ij * b_ijkl
1345 : RankTwoTensorTempl<T> initialContraction(const RankFourTensorTempl<T> & b) const;
1346 :
1347 : /// @}
1348 :
1349 : /// @{
1350 :
1351 : /// Defines logical equality with another RankTwoTensorTempl<T>
1352 : bool operator==(const RankTwoTensorTempl<T> & a) const;
1353 :
1354 : /// Test for symmetry
1355 : bool isSymmetric() const;
1356 :
1357 : /// @}
1358 :
1359 : protected:
1360 : /**
1361 : * Uses the petscblaslapack.h LAPACKsyev_ routine to find, for symmetric _coords:
1362 : * (1) the eigenvalues (if calculation_type == "N")
1363 : * (2) the eigenvalues and eigenvectors (if calculation_type == "V")
1364 : * @param calculation_type If "N" then calculation eigenvalues only
1365 : * @param eigvals Eigenvalues are placed in this array, in ascending order
1366 : * @param a Eigenvectors are placed in this array if calculation_type == "V".
1367 : * See code in dsymmetricEigenvalues for extracting eigenvectors from the a output.
1368 : */
1369 : void syev(const char * calculation_type, std::vector<T> & eigvals, std::vector<T> & a) const;
1370 :
1371 : private:
1372 : static constexpr Real identityCoords[N2] = {1, 0, 0, 0, 1, 0, 0, 0, 1};
1373 :
1374 : template <class T2>
1375 : friend void dataStore(std::ostream &, RankTwoTensorTempl<T2> &, void *);
1376 :
1377 : using libMesh::TensorValue<T>::_coords;
1378 :
1379 : template <class T2>
1380 : friend void dataLoad(std::istream &, RankTwoTensorTempl<T2> &, void *);
1381 : template <class T2>
1382 : friend class RankFourTensorTempl;
1383 : template <class T2>
1384 : friend class RankThreeTensorTempl;
1385 : };
1386 :
1387 : namespace MetaPhysicL
1388 : {
1389 : template <typename T>
1390 : struct RawType<RankTwoTensorTempl<T>>
1391 : {
1392 : typedef RankTwoTensorTempl<typename RawType<T>::value_type> value_type;
1393 :
1394 26 : static value_type value(const RankTwoTensorTempl<T> & in)
1395 : {
1396 26 : value_type ret;
1397 104 : for (auto i : libMesh::make_range(RankTwoTensorTempl<T>::N))
1398 312 : for (auto j : libMesh::make_range(RankTwoTensorTempl<T>::N))
1399 234 : ret(i, j) = raw_value(in(i, j));
1400 :
1401 26 : return ret;
1402 0 : }
1403 : };
1404 : }
1405 :
1406 : template <typename T>
1407 : template <typename T2>
1408 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1409 1348 : RankTwoTensorTempl<T>::operator+(const libMesh::TypeTensor<T2> & b) const
1410 : {
1411 1348 : return libMesh::TensorValue<T>::operator+(b);
1412 : }
1413 :
1414 : template <typename T>
1415 : template <typename T2>
1416 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1417 792400 : RankTwoTensorTempl<T>::operator-(const libMesh::TypeTensor<T2> & b) const
1418 : {
1419 792400 : return libMesh::TensorValue<T>::operator-(b);
1420 : }
1421 :
1422 : template <typename T>
1423 : template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type>
1424 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1425 1341 : RankTwoTensorTempl<T>::operator*(const T2 & b) const
1426 : {
1427 1341 : return libMesh::TensorValue<T>::operator*(b);
1428 : }
1429 :
1430 : template <typename T>
1431 : template <typename T2>
1432 : libMesh::TypeVector<typename libMesh::CompareTypes<T, T2>::supertype>
1433 14 : RankTwoTensorTempl<T>::operator*(const libMesh::TypeVector<T2> & b) const
1434 : {
1435 14 : return libMesh::TensorValue<T>::operator*(b);
1436 : }
1437 :
1438 : template <typename T>
1439 : template <typename T2>
1440 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1441 25 : RankTwoTensorTempl<T>::operator*(const libMesh::TypeTensor<T2> & b) const
1442 : {
1443 25 : return libMesh::TensorValue<T>::operator*(b);
1444 : }
1445 :
1446 : template <typename T>
1447 : template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type>
1448 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1449 254 : RankTwoTensorTempl<T>::operator/(const T2 & b) const
1450 : {
1451 254 : return libMesh::TensorValue<T>::operator/(b);
1452 : }
1453 :
1454 : template <typename T>
1455 : RankFourTensorTempl<T>
1456 1 : RankTwoTensorTempl<T>::positiveProjectionEigenDecomposition(std::vector<T> & eigval,
1457 : RankTwoTensorTempl<T> & eigvec) const
1458 : {
1459 : if constexpr (MooseUtils::IsLikeReal<T>::value)
1460 : {
1461 : // Compute eigenvectors and eigenvalues of this tensor
1462 1 : this->symmetricEigenvaluesEigenvectors(eigval, eigvec);
1463 :
1464 : // Separate out positive and negative eigen values
1465 0 : std::array<T, N> epos;
1466 0 : std::array<T, N> d;
1467 4 : for (auto i : libMesh::make_range(N))
1468 : {
1469 3 : epos[i] = (std::abs(eigval[i]) + eigval[i]) / 2.0;
1470 3 : d[i] = 0 < eigval[i] ? 1.0 : 0.0;
1471 : }
1472 :
1473 : // projection tensor
1474 1 : RankFourTensorTempl<T> proj_pos;
1475 1 : RankFourTensorTempl<T> Gab, Gba;
1476 :
1477 4 : for (auto a : libMesh::make_range(N))
1478 : {
1479 3 : const auto Ma = RankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(a));
1480 3 : proj_pos += d[a] * Ma.outerProduct(Ma);
1481 : }
1482 :
1483 : usingTensorIndices(i_, j_, k_, l_);
1484 4 : for (const auto a : libMesh::make_range(N))
1485 6 : for (const auto b : libMesh::make_range(a))
1486 : {
1487 3 : const auto Ma = RankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(a));
1488 3 : const auto Mb = RankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(b));
1489 :
1490 3 : Gab = Ma.template times<i_, k_, j_, l_>(Mb) + Ma.template times<i_, l_, j_, k_>(Mb);
1491 3 : Gba = Mb.template times<i_, k_, j_, l_>(Ma) + Mb.template times<i_, l_, j_, k_>(Ma);
1492 :
1493 0 : T theta_ab;
1494 3 : if (!MooseUtils::absoluteFuzzyEqual(eigval[a], eigval[b]))
1495 3 : theta_ab = 0.5 * (epos[a] - epos[b]) / (eigval[a] - eigval[b]);
1496 : else
1497 0 : theta_ab = 0.25 * (d[a] + d[b]);
1498 :
1499 3 : proj_pos += theta_ab * (Gab + Gba);
1500 : }
1501 1 : return proj_pos;
1502 0 : }
1503 : else
1504 : mooseError("positiveProjectionEigenDecomposition is only available for ordered tensor "
1505 : "component types");
1506 : }
1507 :
1508 : template <typename T>
1509 : T
1510 32 : RankTwoTensorTempl<T>::sin3Lode(const T & r0, const T & r0_value) const
1511 : {
1512 : if constexpr (MooseUtils::IsLikeReal<T>::value)
1513 : {
1514 32 : T bar = secondInvariant();
1515 32 : if (bar <= r0)
1516 : // in this case the Lode angle is not defined
1517 0 : return r0_value;
1518 : else
1519 : // the min and max here gaurd against precision-loss when bar is tiny but nonzero.
1520 32 : return std::max(std::min(-1.5 * std::sqrt(3.0) * thirdInvariant() / std::pow(bar, 1.5), 1.0),
1521 64 : -1.0);
1522 0 : }
1523 : else
1524 : mooseError("sin3Lode is only available for ordered tensor component types");
1525 : }
1526 :
1527 : template <typename T>
1528 : RankTwoTensorTempl<T>
1529 248 : RankTwoTensorTempl<T>::dsin3Lode(const T & r0) const
1530 : {
1531 : if constexpr (MooseUtils::IsLikeReal<T>::value)
1532 : {
1533 248 : T bar = secondInvariant();
1534 248 : if (bar <= r0)
1535 0 : return RankTwoTensorTempl<T>();
1536 : else
1537 0 : return -1.5 * std::sqrt(3.0) *
1538 0 : (dthirdInvariant() / std::pow(bar, 1.5) -
1539 248 : 1.5 * dsecondInvariant() * thirdInvariant() / std::pow(bar, 2.5));
1540 0 : }
1541 : else
1542 : mooseError("dsin3Lode is only available for ordered tensor component types");
1543 : }
1544 :
1545 : template <typename T>
1546 : RankFourTensorTempl<T>
1547 2 : RankTwoTensorTempl<T>::d2sin3Lode(const T & r0) const
1548 : {
1549 : if constexpr (MooseUtils::IsLikeReal<T>::value)
1550 : {
1551 2 : T bar = secondInvariant();
1552 2 : if (bar <= r0)
1553 0 : return RankFourTensorTempl<T>();
1554 :
1555 2 : T J3 = thirdInvariant();
1556 2 : RankTwoTensorTempl<T> dII = dsecondInvariant();
1557 2 : RankTwoTensorTempl<T> dIII = dthirdInvariant();
1558 2 : RankFourTensorTempl<T> deriv = d2thirdInvariant() / std::pow(bar, 1.5) -
1559 2 : 1.5 * d2secondInvariant() * J3 / std::pow(bar, 2.5);
1560 :
1561 8 : for (unsigned i = 0; i < N; ++i)
1562 24 : for (unsigned j = 0; j < N; ++j)
1563 72 : for (unsigned k = 0; k < N; ++k)
1564 216 : for (unsigned l = 0; l < N; ++l)
1565 324 : deriv(i, j, k, l) += (-1.5 * dII(i, j) * dIII(k, l) - 1.5 * dIII(i, j) * dII(k, l)) /
1566 324 : std::pow(bar, 2.5) +
1567 162 : 1.5 * 2.5 * dII(i, j) * dII(k, l) * J3 / std::pow(bar, 3.5);
1568 :
1569 2 : deriv *= -1.5 * std::sqrt(3.0);
1570 2 : return deriv;
1571 2 : }
1572 : else
1573 : mooseError("d2sin3Lode is only available for ordered tensor component types");
1574 : }
1575 :
1576 : template <typename T>
1577 : template <int n, int o, int p, int q>
1578 : RankFourTensorTempl<T>
1579 91 : RankTwoTensorTempl<T>::times(const RankTwoTensorTempl<T> & b) const
1580 : {
1581 91 : RankFourTensorTempl<T> result;
1582 : std::size_t x[4];
1583 364 : for (x[0] = 0; x[0] < N; ++x[0])
1584 1092 : for (x[1] = 0; x[1] < N; ++x[1])
1585 3276 : for (x[2] = 0; x[2] < N; ++x[2])
1586 9828 : for (x[3] = 0; x[3] < N; ++x[3])
1587 7371 : result(x[0], x[1], x[2], x[3]) = (*this)(x[n], x[o]) * b(x[p], x[q]);
1588 :
1589 182 : return result;
1590 0 : }
1591 :
1592 : template <typename T>
1593 : template <int n, int o, int p, int q, int r, int s>
1594 : RankFourTensorTempl<T>
1595 3 : RankTwoTensorTempl<T>::times(const RankFourTensorTempl<T> & b) const
1596 : {
1597 3 : RankFourTensorTempl<T> result;
1598 : std::size_t x[5];
1599 12 : for (x[0] = 0; x[0] < N; ++x[0])
1600 36 : for (x[1] = 0; x[1] < N; ++x[1])
1601 108 : for (x[2] = 0; x[2] < N; ++x[2])
1602 324 : for (x[3] = 0; x[3] < N; ++x[3])
1603 972 : for (x[4] = 0; x[4] < N; ++x[4])
1604 729 : result(x[0], x[1], x[2], x[3]) += (*this)(x[n], x[o]) * b(x[p], x[q], x[r], x[s]);
1605 :
1606 6 : return result;
1607 : }
|
