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 18 : 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 44 : RankTwoTensorTempl(const RankTwoTensorTempl<T> & a) = default;
266 :
267 : /**
268 : * @brief The conversion operator from a `libMesh::TensorValue`
269 : */
270 274992 : RankTwoTensorTempl(const libMesh::TensorValue<T> & a) : libMesh::TensorValue<T>(a) {}
271 :
272 : /**
273 : * @brief The conversion operator from a `libMesh::TypeTensor`
274 : */
275 802938 : 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 4 : RankTwoTensorTempl(const SymmetricRankTwoTensorTempl<T2> & a)
282 0 : : RankTwoTensorTempl<T>(a(0),
283 4 : a(1),
284 4 : a(2),
285 4 : a(3) / MathUtils::sqrt2,
286 4 : a(4) / MathUtils::sqrt2,
287 8 : a(5) / MathUtils::sqrt2)
288 : {
289 4 : }
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 39 : RankTwoTensorTempl(const RankTwoTensorTempl<T2> & a) : libMesh::TensorValue<T>(a)
297 : {
298 39 : }
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 100472 : 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 std::enable_if<libMesh::ScalarTraits<Scalar>::value, RankTwoTensorTempl &>::type
553 : operator=(const Scalar & libmesh_dbg_var(p))
554 : {
555 : libmesh_assert_equal_to(p, Scalar(0));
556 : this->zero();
557 : return *this;
558 : }
559 :
560 : /**
561 : * Add another second order tensor to this one
562 : *
563 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
564 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
565 : * // A = [ 1 2 3
566 : * // 2 4 6
567 : * // 3 6 9 ]
568 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
569 : * // B = [ 9 6 3
570 : * // 8 5 2
571 : * // 7 4 1 ]
572 : * A += B;
573 : * // A = [ 10 8 6
574 : * // 10 9 8
575 : * // 10 10 10 ]
576 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
577 : */
578 : RankTwoTensorTempl<T> & operator+=(const RankTwoTensorTempl<T> & a);
579 :
580 : /**
581 : * Subtract another second order tensor from this one
582 : *
583 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
584 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
585 : * // A = [ 1 2 3
586 : * // 2 4 6
587 : * // 3 6 9 ]
588 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
589 : * // B = [ 9 6 3
590 : * // 8 5 2
591 : * // 7 4 1 ]
592 : * A -= B;
593 : * // A = [ -8 -4 0
594 : * // -6 -1 4
595 : * // -4 2 8 ]
596 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
597 : */
598 : RankTwoTensorTempl<T> & operator-=(const RankTwoTensorTempl<T> & a);
599 :
600 : /**
601 : * Multiply this tensor by a scalar (component-wise)
602 : *
603 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
604 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
605 : * // A = [ 1 2 3
606 : * // 2 4 6
607 : * // 3 6 9 ]
608 : * A *= 2;
609 : * // A = [ 2 4 6
610 : * // 4 8 12
611 : * // 6 12 18 ]
612 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
613 : */
614 : RankTwoTensorTempl<T> & operator*=(const T & a);
615 :
616 : /**
617 : * Divide this tensor by a scalar (component-wise)
618 : *
619 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
620 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
621 : * // A = [ 1 2 3
622 : * // 2 4 6
623 : * // 3 6 9 ]
624 : * A /= 2;
625 : * // A = [ 0.5 1.0 1.5
626 : * // 1.0 2.0 3.0
627 : * // 1.5 3.0 4.5 ]
628 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
629 : */
630 : RankTwoTensorTempl<T> & operator/=(const T & a);
631 :
632 : /**
633 : * Multiplication with another second order tensor
634 : *
635 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
636 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
637 : * // A = [ 1 2 3
638 : * // 2 4 6
639 : * // 3 6 9 ]
640 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
641 : * // B = [ 9 6 3
642 : * // 8 5 2
643 : * // 7 4 1 ]
644 : * A *= B;
645 : * // A = [ 90 54 18
646 : * // 114 69 24
647 : * // 138 84 30 ]
648 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
649 : */
650 : RankTwoTensorTempl<T> & operator*=(const libMesh::TypeTensor<T> & a);
651 :
652 : /**
653 : * @brief The smart mutator that determines how to fill the second order tensor based on the size
654 : * of the input vector.
655 : *
656 : * @param input The input vector, can be of size 1, 3, 6, or 9
657 : * @param fill_method The fill method, default to autodetect.
658 : *
659 : * When `input.size() == 1`, the vector value is used to fill the diagonal components of the
660 : * second order tensor:
661 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
662 : * RankTwoTensor A;
663 : * A.fillFromInputVector({1.5});
664 : * // A = [ 1.5 0 0
665 : * // 0 1.5 0
666 : * // 0 0 1.5 ]
667 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
668 : *
669 : * When `input.size() == 3`, the vector values are used to fill the diagonal components of the
670 : * second order tensor:
671 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
672 : * RankTwoTensor A;
673 : * A.fillFromInputVector({1, 2, 3});
674 : * // A = [ 1 0 0
675 : * // 0 2 0
676 : * // 0 0 3 ]
677 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
678 : *
679 : * When `input.size() == 6`, the second order tensor is filled symmetrically using the Voigt
680 : * notation:
681 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
682 : * RankTwoTensor A;
683 : * A.fillFromInputVector({1, 2, 3, 4, 5, 6});
684 : * // A = [ 1 6 5
685 : * // 6 2 4
686 : * // 5 4 3 ]
687 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
688 : *
689 : * When `input.size() == 9`, all components of the second order tensor are filled in a
690 : * column-major fashion:
691 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
692 : * RankTwoTensor A;
693 : * A.fillFromInputVector({1, 2, 3, 4, 5, 6, 7, 8, 9});
694 : * // A = [ 1 4 7
695 : * // 2 5 8
696 : * // 3 6 9 ]
697 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
698 : */
699 : void fillFromInputVector(const std::vector<T> & input, FillMethod fill_method = autodetect);
700 :
701 : /**
702 : * @brief The smart mutator that determines how to fill the second order tensor based on the order
703 : * of the scalar_variable.
704 : *
705 : * @param scalar_variable The input scalar variable. Supported orders are FIRST, THIRD, and SIXTH.
706 : *
707 : * When `scalar_variable.size() == 1`, the scalar value is used to fill the very first component
708 : * of the second order tensor:
709 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
710 : & // Suppose v[0] = 1
711 : * RankTwoTensor A;
712 : * A.fillFromScalarVariable(v);
713 : * // A = [ 1 0 0
714 : * // 0 0 0
715 : * // 0 0 0 ]
716 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
717 : *
718 : * When `scalar_variable.size() == 3`, the scalar values are used to fill the in-plane components
719 : * of the second order tensor using the Voigt notation:
720 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
721 : * // Suppose v[0] = 1
722 : * // v[1] = 2
723 : * // v[2] = 3
724 : * RankTwoTensor A;
725 : * A.fillFromScalarVariable(v);
726 : * // A = [ 1 3 0
727 : * // 3 2 0
728 : * // 0 0 0 ]
729 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
730 : *
731 : * When `scalar_variable.size() == 6`, the second order tensor is filled symmetrically using the
732 : * Voigt notation:
733 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
734 : * // Suppose v[0] = 1
735 : * // v[1] = 2
736 : * // v[2] = 3
737 : * // v[3] = 4
738 : * // v[4] = 5
739 : * // v[5] = 6
740 : * RankTwoTensor A;
741 : * A.fillFromScalarVariable(v);
742 : * // A = [ 1 6 5
743 : * // 6 2 4
744 : * // 5 4 3 ]
745 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
746 : */
747 : void fillFromScalarVariable(const VariableValue & scalar_variable);
748 :
749 : /**
750 : * sets _coords[0][0], _coords[0][1], _coords[1][0], _coords[1][1] to input,
751 : * and the remainder to zero
752 : */
753 : void surfaceFillFromInputVector(const std::vector<T> & input);
754 :
755 : /**
756 : * @brief Set the values of the second order tensor to be the outer product of two vectors, i.e.
757 : * \f$ A_{ij} = a_i b_j \f$.
758 : *
759 : * Deprecated in favor of outerProduct()
760 : *
761 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
762 : * RealVectorValue a(1, 2, 3);
763 : * RealVectorValue b(4, 5, 6);
764 : * RankTwoTensor A;
765 : * A.vectorOuterProduct(a, b);
766 : * // A = [ 4 5 6
767 : * // 8 10 12
768 : * // 12 15 18 ]
769 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
770 : */
771 : void vectorOuterProduct(const libMesh::TypeVector<T> &, const libMesh::TypeVector<T> &);
772 :
773 : /**
774 : * @brief Assign values to a specific row of the second order tensor.
775 : *
776 : * @param r The row number, r = 0, 1, 2
777 : * @param v The values to be set
778 : *
779 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
780 : * RealVectorValue a(1, 2, 3);
781 : * RankTwoTensor A;
782 : * // A = [ 0 0 0
783 : * // 0 0 0
784 : * // 0 0 0 ]
785 : * A.fillRow(1, a);
786 : * // A = [ 0 0 0
787 : * // 1 2 3
788 : * // 0 0 0 ]
789 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
790 : */
791 : void fillRow(unsigned int r, const libMesh::TypeVector<T> & v);
792 :
793 : /**
794 : * @brief Assign values to a specific column of the second order tensor.
795 : *
796 : * @param c The column number, c = 0, 1, 2
797 : * @param v The values to be set
798 : *
799 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
800 : * RealVectorValue a(1, 2, 3);
801 : * RankTwoTensor A;
802 : * // A = [ 0 0 0
803 : * // 0 0 0
804 : * // 0 0 0 ]
805 : * A.fillColumn(1, a);
806 : * // A = [ 0 1 0
807 : * // 0 2 0
808 : * // 0 3 0 ]
809 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
810 : */
811 : void fillColumn(unsigned int c, const libMesh::TypeVector<T> & v);
812 :
813 : /**
814 : * @brief Set the tensor to identity.
815 : *
816 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
817 : * RankTwoTensor A;
818 : * // A = [ 0 0 0
819 : * // 0 0 0
820 : * // 0 0 0 ]
821 : * A.setToIdentity();
822 : * // A = [ 1 0 0
823 : * // 0 1 0
824 : * // 0 0 1 ]
825 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
826 : */
827 : void setToIdentity();
828 :
829 : /// Add identity times a to _coords
830 : /**
831 : * @brief Add a scalar to diagonal components \f$ A_{ij} + a\delta_{ij} \f$
832 : *
833 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
834 : * RankTwoTensor A;
835 : * // A = [ 0 0 0
836 : * // 0 0 0
837 : * // 0 0 0 ]
838 : * A.addIa(1.5);
839 : * // A = [ 1.5 0 0
840 : * // 0 1.5 0
841 : * // 0 0 1.5 ]
842 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
843 : */
844 : void addIa(const T & a);
845 :
846 : /**
847 : * Rotate the tensor in-place given a rotation tensor
848 : * \f$ A_{ij} \leftarrow R_{ij} A_{jk} R_{jk} \f$
849 : * @param R The rotation tensor
850 : *
851 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
852 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
853 : * // A = [ 1 4 7
854 : * // 2 5 8
855 : * // 3 6 9 ]
856 : * RankTwoTensor R(0, 1, 0, 1, 0, 0, 0, 0, 1);
857 : * // R = [ 0 1 0
858 : * // 1 0 0
859 : * // 0 0 1 ]
860 : * A.rotate(R);
861 : * // A = [ 5 2 8
862 : * // 4 1 7
863 : * // 6 3 9 ]
864 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
865 : */
866 : void rotate(const RankTwoTensorTempl<T> & R);
867 :
868 : /// @}
869 :
870 : /// @{
871 :
872 : /**
873 : * @brief Return \f$ A_{ij} = A_{ik}A_{kj} \f$
874 : *
875 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
876 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
877 : * // A = [ 1 4 7
878 : * // 2 5 8
879 : * // 3 6 9 ]
880 : * RankTwoTensor B = A.square();
881 : * // B = [ 30 66 102
882 : * // 36 81 126
883 : * // 42 96 150 ]
884 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
885 : */
886 : RankTwoTensorTempl<T> square() const;
887 :
888 : /**
889 : * Return the rotated tensor given a rotation tensor
890 : * \f$ A'_{ij} = R_{ij} A_{jk} R_{jk} \f$
891 : * @param R The rotation tensor
892 : *
893 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
894 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
895 : * // A = [ 1 4 7
896 : * // 2 5 8
897 : * // 3 6 9 ]
898 : * RankTwoTensor R(0, 1, 0, 1, 0, 0, 0, 0, 1);
899 : * // R = [ 0 1 0
900 : * // 1 0 0
901 : * // 0 0 1 ]
902 : * RankTwoTensor A_rotated = A.rotated(R);
903 : * // A_rotated = [ 5 2 8
904 : * // 4 1 7
905 : * // 6 3 9 ]
906 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
907 : */
908 : RankTwoTensorTempl<T> rotated(const RankTwoTensorTempl<T> & R) const;
909 :
910 : /**
911 : * Rotate the tensor about the z-axis
912 : * @param a The rotation angle in radians
913 : */
914 : RankTwoTensorTempl<T> rotateXyPlane(T a);
915 :
916 : /**
917 : * Return the tensor transposed
918 : *
919 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
920 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
921 : * // A = [ 1 4 7
922 : * // 2 5 8
923 : * // 3 6 9 ]
924 : * RankTwoTensor At = A.transpose();
925 : * // At = [ 1 2 3
926 : * // 4 5 6
927 : * // 7 8 9 ]
928 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
929 : */
930 : RankTwoTensorTempl<T> transpose() const;
931 :
932 : /**
933 : * Return the sum of two second order tensors
934 : *
935 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
936 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
937 : * // A = [ 1 2 3
938 : * // 2 4 6
939 : * // 3 6 9 ]
940 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
941 : * // B = [ 9 6 3
942 : * // 8 5 2
943 : * // 7 4 1 ]
944 : * RankTwoTensor C = A + B;
945 : * // C = [ 10 8 6
946 : * // 10 9 8
947 : * // 10 10 10 ]
948 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
949 : */
950 : template <typename T2>
951 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
952 : operator+(const libMesh::TypeTensor<T2> & a) const;
953 :
954 : /**
955 : * Return the subtraction of two second order tensors
956 : *
957 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
958 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
959 : * // A = [ 1 2 3
960 : * // 2 4 6
961 : * // 3 6 9 ]
962 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
963 : * // B = [ 9 6 3
964 : * // 8 5 2
965 : * // 7 4 1 ]
966 : * RankTwoTensor C = A - B;
967 : * // C = [ -8 -4 0
968 : * // -6 -1 4
969 : * // -4 2 8 ]
970 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
971 : */
972 : template <typename T2>
973 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
974 : operator-(const libMesh::TypeTensor<T2> & a) const;
975 :
976 : /**
977 : * Return the negation of this tensor
978 : *
979 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
980 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
981 : * // A = [ 1 2 3
982 : * // 2 4 6
983 : * // 3 6 9 ]
984 : * RankTwoTensor B = -A;
985 : * // B = [ -1 -2 -3
986 : * // -2 -4 -6
987 : * // -3 -6 -9 ]
988 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
989 : */
990 : RankTwoTensorTempl<T> operator-() const;
991 :
992 : /**
993 : * Return this tensor multiplied by a scalar (component-wise)
994 : *
995 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
996 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
997 : * // A = [ 1 2 3
998 : * // 2 4 6
999 : * // 3 6 9 ]
1000 : * RankTwoTensor B = A * 2;
1001 : * // B = [ 2 4 6
1002 : * // 4 8 12
1003 : * // 6 12 18 ]
1004 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1005 : */
1006 : template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type = 0>
1007 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1008 : operator*(const T2 & a) const;
1009 :
1010 : /**
1011 : * Return this tensor divided by a scalar (component-wise)
1012 : *
1013 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1014 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
1015 : * // A = [ 1 2 3
1016 : * // 2 4 6
1017 : * // 3 6 9 ]
1018 : * RankTwoTensor B = A / 2;
1019 : * // B = [ 0.5 1.0 1.5
1020 : * // 1.0 2.0 3.0
1021 : * // 1.5 3.0 4.5 ]
1022 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1023 : */
1024 : template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type = 0>
1025 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1026 : operator/(const T2 & a) const;
1027 :
1028 : /**
1029 : * Return this tensor multiplied by a vector. \f$ b_i = A_{ij} a_j \f$
1030 : *
1031 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1032 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
1033 : * // A = [ 1 2 3
1034 : * // 2 4 6
1035 : * // 3 6 9 ]
1036 : * RealVectorValue a(1, 2, 3);
1037 : * // a = [ 1
1038 : * // 2
1039 : * // 3 ]
1040 : * RealVectorValue b = A * a;
1041 : * // b = [ 30
1042 : * // 36
1043 : * // 42 ]
1044 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1045 : */
1046 : template <typename T2>
1047 : libMesh::TypeVector<typename libMesh::CompareTypes<T, T2>::supertype>
1048 : operator*(const libMesh::TypeVector<T2> & a) const;
1049 :
1050 : /**
1051 : * Multiplication with another second order tensor
1052 : *
1053 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1054 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
1055 : * // A = [ 1 2 3
1056 : * // 2 4 6
1057 : * // 3 6 9 ]
1058 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
1059 : * // B = [ 9 6 3
1060 : * // 8 5 2
1061 : * // 7 4 1 ]
1062 : * RankTwoTensor C = A * B;
1063 : * // C = [ 90 54 18
1064 : * // 114 69 24
1065 : * // 138 84 30 ]
1066 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1067 : */
1068 : template <typename T2>
1069 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1070 : operator*(const libMesh::TypeTensor<T2> & a) const;
1071 :
1072 : /**
1073 : * Return the double contraction with another second order tensor \f$ A_{ij} B_{ij} \f$
1074 : *
1075 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1076 : * RankTwoTensor A(1, 2, 3, 4, 5, 6, 7, 8, 9);
1077 : * // A = [ 1 2 3
1078 : * // 2 4 6
1079 : * // 3 6 9 ]
1080 : * RankTwoTensor B(9, 8, 7, 6, 5, 4, 3, 2, 1);
1081 : * // B = [ 9 6 3
1082 : * // 8 5 2
1083 : * // 7 4 1 ]
1084 : * Real result = A.doubleContraction(B);
1085 : * // result = 143
1086 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1087 : */
1088 :
1089 : T doubleContraction(const RankTwoTensorTempl<T> & a) const;
1090 :
1091 : /**
1092 : * Return the general tensor product of this second order tensor and another second order tensor
1093 : * defined as \f$ C_{ijkl} = A_{\mathcal{M}(n)\mathcal{M}(o)} B_{\mathcal{M}(p)\mathcal{M}(q)} \f$
1094 : * where the multiplication order is defined by the index map \f$ \mathcal{M}: \{n,o,p,q\} \to
1095 : * \{i,j,k,l\} \f$. The index map is specified using the template parameters. See examples below
1096 : * for detailed explanation.
1097 : *
1098 : * Suppose we have two second order tensors A and B, and we denote the output indices as i, j, k,
1099 : * l:
1100 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1101 : * usingTensorIndices(i, j, k, l);
1102 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1103 : *
1104 : * The outer product of A and B is defined as \f$ A_{ij} B_{kl} \f$, hence the template
1105 : * specialization should be `times<i, j, k, l>`
1106 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1107 : * RankFourTensor C = A.times<i, j, k, l>(B);
1108 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1109 : *
1110 : * The tensor product of A and B, i.e. \f$ A_{ik} B_{jl} \f$ can be expressed using the template
1111 : * specialization `times<i, k, j, l>`
1112 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1113 : * RankFourTensor C = A.times<i, k, j, l>(B);
1114 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1115 : *
1116 : * Similarly, another tensor product of A and B, i.e. \f$ A_{il} B_{jk} \f$ can be expressed using
1117 : * the template specialization `times<i, l, j, k>`
1118 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1119 : * RankFourTensor C = A.times<i, l, j, k>(B);
1120 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1121 : *
1122 : * The combination goes on...
1123 : */
1124 : template <int n, int o, int p, int q>
1125 : RankFourTensorTempl<T> times(const RankTwoTensorTempl<T> & b) const;
1126 :
1127 : /**
1128 : * Return the single contraction of this second order tensor with a fourth order tensor defined as
1129 : * \f$ C_{ijkl} = A_{\mathcal{M}(m)\mathcal{M}(n)}
1130 : * B_{\mathcal{M}(p)\mathcal{M}(q)\mathcal{M}(r)\mathcal{M}(s)} \f$ where the multiplication order
1131 : * is defined by the index map \f$ \mathcal{M}: \{m,n,p,q,r,s\} \to \{i,j,k,l\} \f$. The index map
1132 : * is specified using the template parameters. See examples below for detailed explanation.
1133 : *
1134 : * Suppose we have a second order tensors A and a fourth order tensor B, and we denote the indices
1135 : * (four output indices and a dummy index to be contracted) as i, j, k, l, m:
1136 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1137 : * usingTensorIndices(i, j, k, l, m);
1138 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1139 : *
1140 : * The single contraction of A and B defined as \f$ A_{im} B_{mjkl} \f$ can be expressed using the
1141 : * template specialization `times<i, m, m, j, k, l>`
1142 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1143 : * RankFourTensor C = A.times<i, m, m, j, k, l>(B);
1144 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1145 : *
1146 : * Similarly, another single contraction of A and B, i.e. \f$ A_{m, i} A_{j, k, m, l} \f$ can be
1147 : * expressed using the template specialization `times<m, i, j, k, m, l>`
1148 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp}
1149 : * RankFourTensor C = A.times<m, i, j, k, m, l>(B);
1150 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1151 : *
1152 : * The combination goes on. Note that this method assumes exactly one repeated index.
1153 : */
1154 : template <int n, int o, int p, int q, int r, int s>
1155 : RankFourTensorTempl<T> times(const RankFourTensorTempl<T> & b) const;
1156 :
1157 : /**
1158 : * @brief Return the outer product \f$ C_{ijkl} = A_{ij} B_{kl} \f$
1159 : */
1160 38 : RankFourTensorTempl<T> outerProduct(const RankTwoTensorTempl<T> & b) const
1161 : {
1162 : usingTensorIndices(i_, j_, k_, l_);
1163 38 : return times<i_, j_, k_, l_>(b);
1164 : }
1165 :
1166 : /**
1167 : * @brief Return the single contraction of this second order tensor with a third order tensor \f$
1168 : * C_{ikl} = A_{ij} B_{jkl} \f$
1169 : */
1170 : RankThreeTensorTempl<T> contraction(const RankThreeTensorTempl<T> & b) const;
1171 :
1172 : /**
1173 : * @brief Return the tensor product of this second order tensor with a vector \f$ C_{ijk} = A_{jk}
1174 : * b_{i} \f$
1175 : */
1176 : RankThreeTensorTempl<T> mixedProductJkI(const libMesh::VectorValue<T> & b) const;
1177 :
1178 : /**
1179 : * @brief Return the positive projection tensor
1180 : *
1181 : * Consider the eigenvalue decomposition of this second order tensor \f$ A = V D V^T \f$, the part
1182 : * of this tensor that lies on the positive spectrum is defined as \f$ A_+ = V \left<D\right> V^T
1183 : * \f$ where the angled brackets are the Macaulay brackets. The positive projection tensor is the
1184 : * linear projection from the full spectrum to the positive spectrum, i.e. \f$ A_+ = P A \f$. The
1185 : * derivation of this positive projection tensor can be found in C. Miehe and M. Lambrecht,
1186 : * Commun. Numer. Meth. Engng 2001; 17:337~353
1187 : *
1188 : * @param eigvals The three eigenvalues of this second order tensor will be filled into this
1189 : * vector.
1190 : * @param eigvecs The three eigenvectors of this second order tensor will be filled into this
1191 : * tensor.
1192 : * @return The fourth order positive projection tensor.
1193 : */
1194 : RankFourTensorTempl<T> positiveProjectionEigenDecomposition(std::vector<T> &,
1195 : RankTwoTensorTempl<T> &) const;
1196 :
1197 : /**
1198 : * @brief Return the deviatoric part of this tensor \f$ A_{ij} - \frac{1}{3} A_{kk} \delta_{ij}
1199 : * \f$
1200 : */
1201 : RankTwoTensorTempl<T> deviatoric() const;
1202 :
1203 : /**
1204 : * @brief A wrapper for tr()
1205 : * @see tr()
1206 : */
1207 : T trace() const;
1208 :
1209 : /**
1210 : * Return the derivative of the trace w.r.t. this second order tensor itself \f$ \frac{\partial
1211 : * A_{kk}}{\partial A_{ij}} = \delta_{ij} \f$
1212 : */
1213 : RankTwoTensorTempl<T> dtrace() const;
1214 :
1215 : /**
1216 : * @brief Return the inverse of this second order tensor
1217 : */
1218 : RankTwoTensorTempl<T> inverse() const;
1219 :
1220 : /**
1221 : * @brief Return the principal second invariant of this second order tensor
1222 : *
1223 : * \f$ I_2 = \frac{1}{2} \left( A_{kk}^2 - A_{ij}A_{ij} \right) \f$
1224 : */
1225 : T generalSecondInvariant() const;
1226 :
1227 : /**
1228 : * @brief Return the main second invariant of this second order tensor
1229 : *
1230 : * \f$ J_2 = \frac{1}{2} \left( S_{ij}S_{ij} \right) \f$, where \f$ S_{ij} = A_{ij} -
1231 : * \frac{1}{3}A_{kk}\delta_{ij} \f$
1232 : */
1233 : T secondInvariant() const;
1234 :
1235 : /**
1236 : * Return the derivative of the main second invariant w.r.t. this second order tensor itself \f$
1237 : * \frac{\partial J_2}{\partial A_{ij}} = S_{ij} = A_{ij} -
1238 : * \frac{1}{3}A_{kk}\delta_{ij} \f$
1239 : */
1240 : RankTwoTensorTempl<T> dsecondInvariant() const;
1241 :
1242 : /**
1243 : * Return the second derivative of the main second invariant w.r.t. this second order tensor
1244 : * itself \f$ \frac{\partial^2 J_2}{\partial A_{ij}A_{kl}} = \delta_{ik}\delta_{jl} -
1245 : * \frac{1}{3}\delta_{ij}\delta_{kl} \f$
1246 : */
1247 : RankFourTensorTempl<T> d2secondInvariant() const;
1248 :
1249 : /**
1250 : * @brief Sin(3*Lode_angle)
1251 : *
1252 : * If secondInvariant() <= r0 then return r0_value
1253 : * This is to gaurd against precision-loss errors.
1254 : * Note that sin(3*Lode_angle) is not defined for secondInvariant() = 0
1255 : */
1256 : T sin3Lode(const T & r0, const T & r0_value) const;
1257 :
1258 : /**
1259 : * d(sin3Lode)/dA_ij
1260 : * If secondInvariant() <= r0 then return zero
1261 : * This is to gaurd against precision-loss errors.
1262 : * Note that sin(3*Lode_angle) is not defined for secondInvariant() = 0
1263 : */
1264 : RankTwoTensorTempl<T> dsin3Lode(const T & r0) const;
1265 :
1266 : /**
1267 : * d^2(sin3Lode)/dA_ij/dA_kl
1268 : * If secondInvariant() <= r0 then return zero
1269 : * This is to gaurd against precision-loss errors.
1270 : * Note that sin(3*Lode_angle) is not defined for secondInvariant() = 0
1271 : */
1272 : RankFourTensorTempl<T> d2sin3Lode(const T & r0) const;
1273 :
1274 : /**
1275 : * Denote the _coords[i][j] by A_ij, then
1276 : * S_ij = A_ij - de_ij*tr(A)/3
1277 : * Then this returns det(S + S.transpose())/2
1278 : * Note the explicit symmeterisation
1279 : */
1280 : T thirdInvariant() const;
1281 :
1282 : /**
1283 : * Denote the _coords[i][j] by A_ij, then
1284 : * this returns d(thirdInvariant()/dA_ij
1285 : */
1286 : RankTwoTensorTempl<T> dthirdInvariant() const;
1287 :
1288 : /**
1289 : * Denote the _coords[i][j] by A_ij, then this returns
1290 : * d^2(thirdInvariant)/dA_ij/dA_kl
1291 : */
1292 : RankFourTensorTempl<T> d2thirdInvariant() const;
1293 :
1294 : /**
1295 : * Denote the _coords[i][j] by A_ij, then this returns
1296 : * d(det)/dA_ij
1297 : */
1298 : RankTwoTensorTempl<T> ddet() const;
1299 :
1300 : /// Sqrt(_coords[i][j]*_coords[i][j])
1301 : T L2norm() const;
1302 :
1303 : /**
1304 : * computes eigenvalues, assuming tens is symmetric, and places them
1305 : * in ascending order in eigvals
1306 : */
1307 : void symmetricEigenvalues(std::vector<T> & eigvals) const;
1308 :
1309 : /**
1310 : * computes eigenvalues and eigenvectors, assuming tens is symmetric, and places them
1311 : * in ascending order in eigvals. eigvecs is a matrix with the first column
1312 : * being the first eigenvector, the second column being the second, etc.
1313 : */
1314 : void symmetricEigenvaluesEigenvectors(std::vector<T> & eigvals,
1315 : RankTwoTensorTempl<T> & eigvecs) const;
1316 :
1317 : /**
1318 : * computes eigenvalues, and their symmetric derivatives wrt vals,
1319 : * assuming tens is symmetric
1320 : * @param eigvals are the eigenvalues of the matrix, in ascending order
1321 : * @param deigvals Here digvals[i](j,k) = (1/2)*(d(eigvals[i])/dA_jk + d(eigvals[i]/dA_kj))
1322 : * Note the explicit symmeterisation here.
1323 : * For equal eigenvalues, these derivatives are not gauranteed to
1324 : * be the ones you expect, since the derivatives in this case are
1325 : * often defined by continuation from the un-equal case, and that is
1326 : * too sophisticated for this routine.
1327 : */
1328 : void dsymmetricEigenvalues(std::vector<T> & eigvals,
1329 : std::vector<RankTwoTensorTempl<T>> & deigvals) const;
1330 :
1331 : /**
1332 : * Computes second derivatives of Eigenvalues of a rank two tensor
1333 : * @param deriv store second derivative of the current tensor in here
1334 : */
1335 : void d2symmetricEigenvalues(std::vector<RankFourTensorTempl<T>> & deriv) const;
1336 :
1337 : /**
1338 : * Uses the petscblaslapack.h LAPACKsyev_ routine to perform RU decomposition and obtain the
1339 : * rotation tensor.
1340 : */
1341 : void getRUDecompositionRotation(RankTwoTensorTempl<T> & rot) const;
1342 :
1343 : /// returns this_ij * b_ijkl
1344 : RankTwoTensorTempl<T> initialContraction(const RankFourTensorTempl<T> & b) const;
1345 :
1346 : /// @}
1347 :
1348 : /// @{
1349 :
1350 : /// Defines logical equality with another RankTwoTensorTempl<T>
1351 : bool operator==(const RankTwoTensorTempl<T> & a) const;
1352 :
1353 : /// Test for symmetry
1354 : bool isSymmetric() const;
1355 :
1356 : /// @}
1357 :
1358 : protected:
1359 : /**
1360 : * Uses the petscblaslapack.h LAPACKsyev_ routine to find, for symmetric _coords:
1361 : * (1) the eigenvalues (if calculation_type == "N")
1362 : * (2) the eigenvalues and eigenvectors (if calculation_type == "V")
1363 : * @param calculation_type If "N" then calculation eigenvalues only
1364 : * @param eigvals Eigenvalues are placed in this array, in ascending order
1365 : * @param a Eigenvectors are placed in this array if calculation_type == "V".
1366 : * See code in dsymmetricEigenvalues for extracting eigenvectors from the a output.
1367 : */
1368 : void syev(const char * calculation_type, std::vector<T> & eigvals, std::vector<T> & a) const;
1369 :
1370 : private:
1371 : static constexpr Real identityCoords[N2] = {1, 0, 0, 0, 1, 0, 0, 0, 1};
1372 :
1373 : template <class T2>
1374 : friend void dataStore(std::ostream &, RankTwoTensorTempl<T2> &, void *);
1375 :
1376 : using libMesh::TensorValue<T>::_coords;
1377 :
1378 : template <class T2>
1379 : friend void dataLoad(std::istream &, RankTwoTensorTempl<T2> &, void *);
1380 : template <class T2>
1381 : friend class RankFourTensorTempl;
1382 : template <class T2>
1383 : friend class RankThreeTensorTempl;
1384 : };
1385 :
1386 : namespace MetaPhysicL
1387 : {
1388 : template <typename T>
1389 : struct RawType<RankTwoTensorTempl<T>>
1390 : {
1391 : typedef RankTwoTensorTempl<typename RawType<T>::value_type> value_type;
1392 :
1393 52 : static value_type value(const RankTwoTensorTempl<T> & in)
1394 : {
1395 52 : value_type ret;
1396 208 : for (auto i : libMesh::make_range(RankTwoTensorTempl<T>::N))
1397 624 : for (auto j : libMesh::make_range(RankTwoTensorTempl<T>::N))
1398 468 : ret(i, j) = raw_value(in(i, j));
1399 :
1400 52 : return ret;
1401 0 : }
1402 : };
1403 : }
1404 :
1405 : template <typename T>
1406 : template <typename T2>
1407 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1408 2696 : RankTwoTensorTempl<T>::operator+(const libMesh::TypeTensor<T2> & b) const
1409 : {
1410 2696 : return libMesh::TensorValue<T>::operator+(b);
1411 : }
1412 :
1413 : template <typename T>
1414 : template <typename T2>
1415 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1416 793312 : RankTwoTensorTempl<T>::operator-(const libMesh::TypeTensor<T2> & b) const
1417 : {
1418 793312 : return libMesh::TensorValue<T>::operator-(b);
1419 : }
1420 :
1421 : template <typename T>
1422 : template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type>
1423 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1424 2682 : RankTwoTensorTempl<T>::operator*(const T2 & b) const
1425 : {
1426 2682 : return libMesh::TensorValue<T>::operator*(b);
1427 : }
1428 :
1429 : template <typename T>
1430 : template <typename T2>
1431 : libMesh::TypeVector<typename libMesh::CompareTypes<T, T2>::supertype>
1432 28 : RankTwoTensorTempl<T>::operator*(const libMesh::TypeVector<T2> & b) const
1433 : {
1434 28 : return libMesh::TensorValue<T>::operator*(b);
1435 : }
1436 :
1437 : template <typename T>
1438 : template <typename T2>
1439 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1440 50 : RankTwoTensorTempl<T>::operator*(const libMesh::TypeTensor<T2> & b) const
1441 : {
1442 50 : return libMesh::TensorValue<T>::operator*(b);
1443 : }
1444 :
1445 : template <typename T>
1446 : template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type>
1447 : RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
1448 508 : RankTwoTensorTempl<T>::operator/(const T2 & b) const
1449 : {
1450 508 : return libMesh::TensorValue<T>::operator/(b);
1451 : }
1452 :
1453 : template <typename T>
1454 : RankFourTensorTempl<T>
1455 2 : RankTwoTensorTempl<T>::positiveProjectionEigenDecomposition(std::vector<T> & eigval,
1456 : RankTwoTensorTempl<T> & eigvec) const
1457 : {
1458 : using std::abs;
1459 : if constexpr (MooseUtils::IsLikeReal<T>::value)
1460 : {
1461 : // Compute eigenvectors and eigenvalues of this tensor
1462 2 : 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 8 : for (auto i : libMesh::make_range(N))
1468 : {
1469 6 : epos[i] = (abs(eigval[i]) + eigval[i]) / 2.0;
1470 6 : d[i] = 0 < eigval[i] ? 1.0 : 0.0;
1471 : }
1472 :
1473 : // projection tensor
1474 2 : RankFourTensorTempl<T> proj_pos;
1475 2 : RankFourTensorTempl<T> Gab, Gba;
1476 :
1477 8 : for (auto a : libMesh::make_range(N))
1478 : {
1479 6 : const auto Ma = RankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(a));
1480 6 : proj_pos += d[a] * Ma.outerProduct(Ma);
1481 : }
1482 :
1483 : usingTensorIndices(i_, j_, k_, l_);
1484 8 : for (const auto a : libMesh::make_range(N))
1485 12 : for (const auto b : libMesh::make_range(a))
1486 : {
1487 6 : const auto Ma = RankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(a));
1488 6 : const auto Mb = RankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(b));
1489 :
1490 6 : Gab = Ma.template times<i_, k_, j_, l_>(Mb) + Ma.template times<i_, l_, j_, k_>(Mb);
1491 6 : Gba = Mb.template times<i_, k_, j_, l_>(Ma) + Mb.template times<i_, l_, j_, k_>(Ma);
1492 :
1493 0 : T theta_ab;
1494 6 : if (!MooseUtils::absoluteFuzzyEqual(eigval[a], eigval[b]))
1495 6 : 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 6 : proj_pos += theta_ab * (Gab + Gba);
1500 : }
1501 2 : 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 64 : RankTwoTensorTempl<T>::sin3Lode(const T & r0, const T & r0_value) const
1511 : {
1512 : using std::pow, std::max, std::sqrt, std::min;
1513 : if constexpr (MooseUtils::IsLikeReal<T>::value)
1514 : {
1515 64 : T bar = secondInvariant();
1516 64 : if (bar <= r0)
1517 : // in this case the Lode angle is not defined
1518 0 : return r0_value;
1519 : else
1520 : // the min and max here gaurd against precision-loss when bar is tiny but nonzero.
1521 64 : return max(min(-1.5 * sqrt(3.0) * thirdInvariant() / pow(bar, 1.5), 1.0), -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 496 : RankTwoTensorTempl<T>::dsin3Lode(const T & r0) const
1530 : {
1531 : using std::sqrt, std::pow;
1532 : if constexpr (MooseUtils::IsLikeReal<T>::value)
1533 : {
1534 496 : T bar = secondInvariant();
1535 496 : if (bar <= r0)
1536 0 : return RankTwoTensorTempl<T>();
1537 : else
1538 0 : return -1.5 * sqrt(3.0) *
1539 0 : (dthirdInvariant() / pow(bar, 1.5) -
1540 496 : 1.5 * dsecondInvariant() * thirdInvariant() / pow(bar, 2.5));
1541 0 : }
1542 : else
1543 : mooseError("dsin3Lode is only available for ordered tensor component types");
1544 : }
1545 :
1546 : template <typename T>
1547 : RankFourTensorTempl<T>
1548 4 : RankTwoTensorTempl<T>::d2sin3Lode(const T & r0) const
1549 : {
1550 : using std::pow, std::sqrt;
1551 : if constexpr (MooseUtils::IsLikeReal<T>::value)
1552 : {
1553 4 : T bar = secondInvariant();
1554 4 : if (bar <= r0)
1555 0 : return RankFourTensorTempl<T>();
1556 :
1557 4 : T J3 = thirdInvariant();
1558 4 : RankTwoTensorTempl<T> dII = dsecondInvariant();
1559 4 : RankTwoTensorTempl<T> dIII = dthirdInvariant();
1560 0 : RankFourTensorTempl<T> deriv =
1561 4 : d2thirdInvariant() / pow(bar, 1.5) - 1.5 * d2secondInvariant() * J3 / pow(bar, 2.5);
1562 :
1563 16 : for (unsigned i = 0; i < N; ++i)
1564 48 : for (unsigned j = 0; j < N; ++j)
1565 144 : for (unsigned k = 0; k < N; ++k)
1566 432 : for (unsigned l = 0; l < N; ++l)
1567 324 : deriv(i, j, k, l) +=
1568 324 : (-1.5 * dII(i, j) * dIII(k, l) - 1.5 * dIII(i, j) * dII(k, l)) / pow(bar, 2.5) +
1569 324 : 1.5 * 2.5 * dII(i, j) * dII(k, l) * J3 / pow(bar, 3.5);
1570 :
1571 4 : deriv *= -1.5 * sqrt(3.0);
1572 4 : return deriv;
1573 4 : }
1574 : else
1575 : mooseError("d2sin3Lode is only available for ordered tensor component types");
1576 : }
1577 :
1578 : template <typename T>
1579 : template <int n, int o, int p, int q>
1580 : RankFourTensorTempl<T>
1581 182 : RankTwoTensorTempl<T>::times(const RankTwoTensorTempl<T> & b) const
1582 : {
1583 182 : RankFourTensorTempl<T> result;
1584 : std::size_t x[4];
1585 728 : for (x[0] = 0; x[0] < N; ++x[0])
1586 2184 : for (x[1] = 0; x[1] < N; ++x[1])
1587 6552 : for (x[2] = 0; x[2] < N; ++x[2])
1588 19656 : for (x[3] = 0; x[3] < N; ++x[3])
1589 14742 : result(x[0], x[1], x[2], x[3]) = (*this)(x[n], x[o]) * b(x[p], x[q]);
1590 :
1591 364 : return result;
1592 0 : }
1593 :
1594 : template <typename T>
1595 : template <int n, int o, int p, int q, int r, int s>
1596 : RankFourTensorTempl<T>
1597 6 : RankTwoTensorTempl<T>::times(const RankFourTensorTempl<T> & b) const
1598 : {
1599 6 : RankFourTensorTempl<T> result;
1600 : std::size_t x[5];
1601 24 : for (x[0] = 0; x[0] < N; ++x[0])
1602 72 : for (x[1] = 0; x[1] < N; ++x[1])
1603 216 : for (x[2] = 0; x[2] < N; ++x[2])
1604 648 : for (x[3] = 0; x[3] < N; ++x[3])
1605 1944 : for (x[4] = 0; x[4] < N; ++x[4])
1606 1458 : result(x[0], x[1], x[2], x[3]) += (*this)(x[n], x[o]) * b(x[p], x[q], x[r], x[s]);
1607 :
1608 12 : return result;
1609 : }
|
