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 "ADRankTwoTensorForward.h"
14 : #include "ADSymmetricRankTwoTensorForward.h"
15 : #include "ADSymmetricRankFourTensorForward.h"
16 : #include "MooseUtils.h"
17 : #include "MathUtils.h"
18 :
19 : #include "libmesh/libmesh.h"
20 : #include "libmesh/tensor_value.h"
21 :
22 : #include "metaphysicl/raw_type.h"
23 :
24 : #include <petscsys.h>
25 : #include <vector>
26 :
27 : // Forward declarations
28 : class MooseEnum;
29 : template <typename T>
30 : class MooseArray;
31 : typedef MooseArray<Real> VariableValue;
32 : template <typename>
33 : class ColumnMajorMatrixTempl;
34 : namespace libMesh
35 : {
36 : template <typename>
37 : class VectorValue;
38 : template <typename>
39 : class TypeVector;
40 : template <typename>
41 : class TypeTensor;
42 : template <typename>
43 : class TensorValue;
44 : }
45 :
46 : namespace MathUtils
47 : {
48 : template <typename T>
49 : void mooseSetToZero(T & v);
50 :
51 : /**
52 : * Helper function template specialization to set an object to zero.
53 : * Needed by DerivativeMaterialInterface
54 : */
55 : template <>
56 : void mooseSetToZero<SymmetricRankTwoTensor>(SymmetricRankTwoTensor & v);
57 :
58 : /**
59 : * Helper function template specialization to set an object to zero.
60 : * Needed by DerivativeMaterialInterface
61 : */
62 : template <>
63 : void mooseSetToZero<ADSymmetricRankTwoTensor>(ADSymmetricRankTwoTensor & v);
64 : }
65 :
66 : /**
67 : * SymmetricRankTwoTensorTempl is designed to handle the Stress or Strain Tensor for
68 : * an anisotropic material. It is designed to reduce the redundancies of the
69 : * Complete tensor classes for regular mechanics problems and to enable use of the
70 : * Mandel notation.
71 : */
72 : template <typename T>
73 : class SymmetricRankTwoTensorTempl
74 : {
75 : public:
76 : /// For generic programming
77 : typedef T value_type;
78 :
79 : ///@{ tensor dimension and Mandel vector length
80 : static constexpr unsigned int Ndim = Moose::dim;
81 : static constexpr unsigned int N = Ndim + (Ndim * (Ndim - 1)) / 2;
82 : ///@}
83 :
84 : // Full tensor indices in the Mandel/Voigt representation
85 : static constexpr unsigned int full_index[6][2] = {{0, 0}, {1, 1}, {2, 2}, {1, 2}, {0, 2}, {0, 1}};
86 :
87 : // Reverse tensor indices in the Mandel/Voigt representation
88 : static constexpr unsigned int reverse_index[3][3] = {{0, 5, 4}, {5, 1, 3}, {4, 3, 2}};
89 :
90 : /// returns the 1 or sqrt(2) prefactor in the Mandel notation for the index i ranging from 0-5.
91 28952 : static constexpr Real mandelFactor(unsigned int i) { return i < Ndim ? 1.0 : MathUtils::sqrt2; }
92 :
93 : // Select initialization
94 : enum InitMethod
95 : {
96 : initNone,
97 : initIdentity
98 : };
99 :
100 : /// Default constructor; fills to zero
101 : SymmetricRankTwoTensorTempl();
102 :
103 : /// Select specific initialization pattern
104 : SymmetricRankTwoTensorTempl(const InitMethod);
105 :
106 : /**
107 : * To fill up the 6 entries in the 2nd-order tensor, fillFromInputVector
108 : * is called with one of the following fill_methods.
109 : * See the fill*FromInputVector functions for more details
110 : */
111 : enum FillMethod
112 : {
113 : autodetect = 0,
114 : isotropic1 = 1,
115 : diagonal3 = 3,
116 : symmetric6 = 6
117 : };
118 :
119 : /// Constructor that proxies the fillFromInputVector method
120 0 : SymmetricRankTwoTensorTempl(const std::vector<T> & input) { this->fillFromInputVector(input); };
121 :
122 : /// Initialization list replacement constructors, 6 arguments
123 : SymmetricRankTwoTensorTempl(
124 : const T & S11, const T & S22, const T & S33, const T & S23, const T & S13, const T & S12);
125 :
126 : // explicit cast to a full tensor
127 : explicit operator RankTwoTensorTempl<T>();
128 :
129 : private:
130 : /// Initialization list replacement constructors, 9 arguments (for internal use only)
131 : SymmetricRankTwoTensorTempl(const T & S11,
132 : const T & S21,
133 : const T & S31,
134 : const T & S12,
135 : const T & S22,
136 : const T & S32,
137 : const T & S13,
138 : const T & S23,
139 : const T & S33);
140 :
141 : // internal use named constructor that does not apply the sqrt(2) factors
142 : static SymmetricRankTwoTensorTempl fromRawComponents(
143 : const T & S11, const T & S22, const T & S33, const T & S23, const T & S13, const T & S12);
144 :
145 : public:
146 : /// Copy constructor from TensorValue<T>
147 : explicit SymmetricRankTwoTensorTempl(const TensorValue<T> & a);
148 :
149 : /// Copy constructor from TypeTensor<T>
150 : explicit SymmetricRankTwoTensorTempl(const TypeTensor<T> & a);
151 :
152 : /// Construct from other template
153 : template <typename T2>
154 : SymmetricRankTwoTensorTempl(const SymmetricRankTwoTensorTempl<T2> & a)
155 : {
156 : for (const auto i : make_range(N))
157 : _vals[i] = a(i);
158 : }
159 :
160 : // Named constructors
161 0 : static SymmetricRankTwoTensorTempl identity()
162 : {
163 0 : return SymmetricRankTwoTensorTempl(initIdentity);
164 : }
165 :
166 : /// named constructor for initializing symmetrically
167 : static SymmetricRankTwoTensorTempl
168 : initializeSymmetric(const TypeVector<T> & v0, const TypeVector<T> & v1, const TypeVector<T> & v2);
169 :
170 : /// Static method for use in validParams for getting the "fill_method"
171 : static MooseEnum fillMethodEnum();
172 :
173 : /// fillFromInputVector takes 1, 3, or 6 inputs to fill in the symmmetric Rank-2 tensor.
174 : void fillFromInputVector(const std::vector<T> & input, FillMethod fill_method = autodetect);
175 :
176 : /**
177 : * fillFromScalarVariable takes FIRST/THIRD/SIXTH order scalar variable to fill in the Rank-2
178 : * tensor.
179 : */
180 : void fillFromScalarVariable(const VariableValue & scalar_variable);
181 :
182 : /// Gets the raw value for the index specified. Takes index = 0,1,2,3,4,5
183 2363388 : inline T & operator()(unsigned int i) { return _vals[i]; }
184 :
185 : /**
186 : * Gets the raw value for the index specified. Takes index = 0,1,2
187 : * used for const
188 : */
189 936 : inline T operator()(unsigned int i) const { return _vals[i]; }
190 :
191 : /**
192 : * Gets the value for the index specified. Takes index = 0,1,2
193 : */
194 1746 : inline T operator()(unsigned int i, unsigned int j) const
195 : {
196 1746 : const auto a = reverse_index[i][j];
197 1746 : return _vals[a] / mandelFactor(a);
198 : }
199 :
200 : /// get the specified row of the tensor
201 : libMesh::VectorValue<T> row(const unsigned int n) const;
202 :
203 : /// get the specified column of the tensor
204 0 : libMesh::VectorValue<T> column(const unsigned int n) const { return row(n); }
205 :
206 : /// return the matrix multiplied with its transpose A*A^T (guaranteed symmetric)
207 : [[nodiscard]] static SymmetricRankTwoTensorTempl<T> timesTranspose(const RankTwoTensorTempl<T> &);
208 : [[nodiscard]] static SymmetricRankTwoTensorTempl<T>
209 : timesTranspose(const SymmetricRankTwoTensorTempl<T> &);
210 :
211 : /// return the matrix multiplied with its transpose A^T*A (guaranteed symmetric)
212 : [[nodiscard]] static SymmetricRankTwoTensorTempl<T> transposeTimes(const RankTwoTensorTempl<T> &);
213 : [[nodiscard]] static SymmetricRankTwoTensorTempl<T>
214 : transposeTimes(const SymmetricRankTwoTensorTempl<T> &);
215 :
216 : /// return the matrix plus its transpose A-A^T (guaranteed symmetric)
217 : [[nodiscard]] static SymmetricRankTwoTensorTempl<T> plusTranspose(const RankTwoTensorTempl<T> &);
218 : [[nodiscard]] static SymmetricRankTwoTensorTempl<T>
219 : plusTranspose(const SymmetricRankTwoTensorTempl<T> &);
220 :
221 : /// Returns the matrix squared
222 : SymmetricRankTwoTensorTempl<T> square() const;
223 :
224 : /// Returns the trace
225 30 : T tr() const { return _vals[0] + _vals[1] + _vals[2]; }
226 :
227 : /// Set all components to zero
228 : void zero();
229 :
230 : /**
231 : * rotates the tensor data given a rank two tensor rotation tensor
232 : * _vals[i][j] = R_ij * R_jl * _vals[k][l]
233 : * @param R rotation matrix as a TypeTensor
234 : */
235 : void rotate(const TypeTensor<T> & R);
236 :
237 : /**
238 : * Returns a matrix that is the transpose of the matrix this
239 : * was called on. This is a non-operation.
240 : */
241 : SymmetricRankTwoTensorTempl<T> transpose() const;
242 :
243 : /// sets _vals to a, and returns _vals
244 : template <typename T2>
245 : SymmetricRankTwoTensorTempl<T> & operator=(const SymmetricRankTwoTensorTempl<T2> & a);
246 :
247 : /**
248 : * Assignment-from-scalar operator. Used only to zero out vectors.
249 : */
250 : template <typename Scalar>
251 : typename std::enable_if<libMesh::ScalarTraits<Scalar>::value, SymmetricRankTwoTensorTempl &>::type
252 : operator=(const Scalar & libmesh_dbg_var(p))
253 : {
254 : libmesh_assert_equal_to(p, Scalar(0));
255 : this->zero();
256 : return *this;
257 : }
258 :
259 : /// adds a to _vals
260 : SymmetricRankTwoTensorTempl<T> & operator+=(const SymmetricRankTwoTensorTempl<T> & a);
261 :
262 : /// returns _vals + a
263 : template <typename T2>
264 : SymmetricRankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
265 : operator+(const SymmetricRankTwoTensorTempl<T2> & a) const;
266 :
267 : /// sets _vals -= a and returns vals
268 : SymmetricRankTwoTensorTempl<T> & operator-=(const SymmetricRankTwoTensorTempl<T> & a);
269 :
270 : /// returns _vals - a
271 : template <typename T2>
272 : SymmetricRankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
273 : operator-(const SymmetricRankTwoTensorTempl<T2> & a) const;
274 :
275 : /// returns -_vals
276 : SymmetricRankTwoTensorTempl<T> operator-() const;
277 :
278 : /// performs _vals *= a
279 : SymmetricRankTwoTensorTempl<T> & operator*=(const T & a);
280 :
281 : /// returns _vals*a
282 : template <typename T2>
283 : auto operator*(const T2 & a) const ->
284 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
285 : SymmetricRankTwoTensorTempl<decltype(T() * T2())>>::type;
286 :
287 : /// performs _vals /= a
288 : SymmetricRankTwoTensorTempl<T> & operator/=(const T & a);
289 :
290 : /// returns _vals/a
291 : template <typename T2>
292 : auto operator/(const T2 & a) const ->
293 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
294 : SymmetricRankTwoTensorTempl<decltype(T() / T2())>>::type;
295 :
296 : /// Defines multiplication with a vector to get a vector
297 : template <typename T2>
298 : TypeVector<typename libMesh::CompareTypes<T, T2>::supertype>
299 : operator*(const TypeVector<T2> & a) const;
300 :
301 : /// Defines logical equality with another SymmetricRankTwoTensorTempl<T2>
302 : template <typename T2>
303 : bool operator==(const SymmetricRankTwoTensorTempl<T2> & a) const;
304 :
305 : /// Test for symmetry. Surprisingly this is always true.
306 0 : bool isSymmetric() const { return true; }
307 :
308 : /// Defines logical inequality with another SymmetricRankTwoTensorTempl<T2>
309 : template <typename T2>
310 : bool operator!=(const SymmetricRankTwoTensorTempl<T2> & a) const;
311 :
312 : /// Sets _vals to the values in a ColumnMajorMatrix (must be 3x3)
313 : SymmetricRankTwoTensorTempl<T> & operator=(const ColumnMajorMatrixTempl<T> & a);
314 :
315 : /// returns _vals_ij * a_ij (sum on i, j)
316 : T doubleContraction(const SymmetricRankTwoTensorTempl<T> & a) const;
317 :
318 : /// returns C_ijkl = a_ij * b_kl
319 : SymmetricRankFourTensorTempl<T> outerProduct(const SymmetricRankTwoTensorTempl<T> & a) const;
320 :
321 : /// return positive projection tensor of eigen-decomposition
322 : SymmetricRankFourTensorTempl<T>
323 : positiveProjectionEigenDecomposition(std::vector<T> &, RankTwoTensorTempl<T> &) const;
324 :
325 : /// returns A_ij - de_ij*tr(A)/3, where A are the _vals
326 : SymmetricRankTwoTensorTempl<T> deviatoric() const;
327 :
328 : /// returns the trace of the tensor, ie _vals[i][i] (sum i = 0, 1, 2)
329 : T trace() const;
330 :
331 : /// retuns the inverse of the tensor
332 : SymmetricRankTwoTensorTempl<T> inverse() const;
333 :
334 : /**
335 : * Denote the _vals[i][j] by A_ij, then this returns
336 : * d(trace)/dA_ij
337 : */
338 : SymmetricRankTwoTensorTempl<T> dtrace() const;
339 :
340 : /**
341 : * Calculates the second invariant (I2) of a tensor
342 : */
343 : T generalSecondInvariant() const;
344 :
345 : /**
346 : * Denote the _vals[i][j] by A_ij, then
347 : * S_ij = A_ij - de_ij*tr(A)/3
348 : * Then this returns (S_ij + S_ji)*(S_ij + S_ji)/8
349 : * Note the explicit symmeterisation
350 : */
351 : T secondInvariant() const;
352 :
353 : /**
354 : * Denote the _vals[i][j] by A_ij, then this returns
355 : * d(secondInvariant)/dA_ij
356 : */
357 : SymmetricRankTwoTensorTempl<T> dsecondInvariant() const;
358 :
359 : /**
360 : * Denote the _vals[i][j] by A_ij, then this returns
361 : * d^2(secondInvariant)/dA_ij/dA_kl
362 : */
363 : SymmetricRankFourTensorTempl<T> d2secondInvariant() const;
364 :
365 : /**
366 : * Sin(3*Lode_angle)
367 : * If secondInvariant() <= r0 then return r0_value
368 : * This is to gaurd against precision-loss errors.
369 : * Note that sin(3*Lode_angle) is not defined for secondInvariant() = 0
370 : */
371 : T sin3Lode(const T & r0, const T & r0_value) const;
372 :
373 : /**
374 : * d(sin3Lode)/dA_ij
375 : * If secondInvariant() <= r0 then return zero
376 : * This is to gaurd against precision-loss errors.
377 : * Note that sin(3*Lode_angle) is not defined for secondInvariant() = 0
378 : */
379 : SymmetricRankTwoTensorTempl<T> dsin3Lode(const T & r0) const;
380 :
381 : /**
382 : * Denote the _vals[i][j] by A_ij, then
383 : * S_ij = A_ij - de_ij*tr(A)/3
384 : * Then this returns det(S + S.transpose())/2
385 : * Note the explicit symmeterisation
386 : */
387 : T thirdInvariant() const;
388 :
389 : /**
390 : * Denote the _vals[i][j] by A_ij, then
391 : * this returns d(thirdInvariant()/dA_ij
392 : */
393 : SymmetricRankTwoTensorTempl<T> dthirdInvariant() const;
394 :
395 : /**
396 : * Denote the _vals[i][j] by A_ij, then this returns
397 : * d^2(thirdInvariant)/dA_ij/dA_kl
398 : */
399 : SymmetricRankFourTensorTempl<T> d2thirdInvariant() const;
400 :
401 : // determinant of the tensor
402 : T det() const;
403 :
404 : /**
405 : * Denote the _vals[i][j] by A_ij, then this returns
406 : * d(det)/dA_ij
407 : */
408 : SymmetricRankTwoTensorTempl<T> ddet() const;
409 :
410 : /// Print the rank two tensor
411 : void print(std::ostream & stm = Moose::out) const;
412 :
413 : /// Print the Real part of the ADReal rank two tensor
414 : void printReal(std::ostream & stm = Moose::out) const;
415 :
416 : /// Print the Real part of the ADReal rank two tensor along with its first nDual dual numbers
417 : void printADReal(unsigned int nDual, std::ostream & stm = Moose::out) const;
418 :
419 : friend std::ostream & operator<<(std::ostream & os, const SymmetricRankTwoTensorTempl<T> & t)
420 : {
421 : t.print(os);
422 : return os;
423 : }
424 :
425 : /// Add identity times a to _vals
426 : void addIa(const T & a);
427 :
428 : /// Sqrt(_vals[i][j]*_vals[i][j])
429 : T L2norm() const;
430 :
431 : /**
432 : * sets _vals[0][0], _vals[0][1], _vals[1][0], _vals[1][1] to input,
433 : * and the remainder to zero
434 : */
435 : void surfaceFillFromInputVector(const std::vector<T> & input);
436 :
437 : /**
438 : * computes eigenvalues, assuming tens is symmetric, and places them
439 : * in ascending order in eigvals
440 : */
441 : void symmetricEigenvalues(std::vector<T> & eigvals) const;
442 :
443 : /**
444 : * computes eigenvalues and eigenvectors, assuming tens is symmetric, and places them
445 : * in ascending order in eigvals. eigvecs is a matrix with the first column
446 : * being the first eigenvector, the second column being the second, etc.
447 : */
448 : void symmetricEigenvaluesEigenvectors(std::vector<T> & eigvals,
449 : RankTwoTensorTempl<T> & eigvecs) const;
450 :
451 : /**
452 : * This function initializes random seed based on a user-defined number.
453 : */
454 : static void initRandom(unsigned int);
455 :
456 : /**
457 : * This function generates a random symmetric rank two tensor.
458 : * The first real scales the random number.
459 : * The second real offsets the uniform random number
460 : */
461 : [[nodiscard]] static SymmetricRankTwoTensorTempl<T> genRandomSymmTensor(T, T);
462 :
463 : /// SymmetricRankTwoTensorTempl<T> from outer product of a vector with itself
464 : [[nodiscard]] static SymmetricRankTwoTensorTempl<T> selfOuterProduct(const TypeVector<T> &);
465 :
466 : /// returns this_ij * b_ijkl
467 : SymmetricRankTwoTensorTempl<T>
468 : initialContraction(const SymmetricRankFourTensorTempl<T> & b) const;
469 :
470 : /// set the tensor to the identity matrix
471 : void setToIdentity();
472 :
473 : protected:
474 : /**
475 : * Uses the petscblaslapack.h LAPACKsyev_ routine to find, for symmetric _vals:
476 : * (1) the eigenvalues (if calculation_type == "N")
477 : * (2) the eigenvalues and eigenvectors (if calculation_type == "V")
478 : * @param calculation_type If "N" then calculation eigenvalues only
479 : * @param eigvals Eigenvalues are placed in this array, in ascending order
480 : * @param a Eigenvectors are placed in this array if calculation_type == "V".
481 : * See code in dsymmetricEigenvalues for extracting eigenvectors from the a output.
482 : */
483 : void syev(const char * calculation_type, std::vector<T> & eigvals, std::vector<T> & a) const;
484 :
485 : private:
486 : static constexpr std::array<Real, N> identityCoords = {{1, 1, 1, 0, 0, 0}};
487 :
488 : // tensor components
489 : std::array<T, N> _vals;
490 :
491 : template <class T2>
492 : friend void dataStore(std::ostream &, SymmetricRankTwoTensorTempl<T2> &, void *);
493 :
494 : template <class T2>
495 : friend void dataLoad(std::istream &, SymmetricRankTwoTensorTempl<T2> &, void *);
496 : template <class T2>
497 : friend class SymmetricRankFourTensorTempl;
498 : };
499 :
500 : namespace MetaPhysicL
501 : {
502 : template <typename T>
503 : struct RawType<SymmetricRankTwoTensorTempl<T>>
504 : {
505 : typedef SymmetricRankTwoTensorTempl<typename RawType<T>::value_type> value_type;
506 :
507 2 : static value_type value(const SymmetricRankTwoTensorTempl<T> & in)
508 : {
509 2 : value_type ret;
510 14 : for (unsigned int i = 0; i < SymmetricRankTwoTensorTempl<T>::N; ++i)
511 12 : ret(i) = raw_value(in(i));
512 :
513 2 : return ret;
514 : }
515 : };
516 : }
517 :
518 : template <typename T>
519 : template <typename T2>
520 : auto
521 50 : SymmetricRankTwoTensorTempl<T>::operator*(const T2 & a) const ->
522 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
523 : SymmetricRankTwoTensorTempl<decltype(T() * T2())>>::type
524 : {
525 50 : SymmetricRankTwoTensorTempl<decltype(T() * T2())> result;
526 350 : for (const auto i : make_range(N))
527 300 : result._vals[i] = _vals[i] * a;
528 50 : return result;
529 0 : }
530 :
531 : template <typename T>
532 : template <typename T2>
533 : auto
534 8 : SymmetricRankTwoTensorTempl<T>::operator/(const T2 & a) const ->
535 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
536 : SymmetricRankTwoTensorTempl<decltype(T() / T2())>>::type
537 : {
538 8 : SymmetricRankTwoTensorTempl<decltype(T() / T2())> result;
539 56 : for (const auto i : make_range(N))
540 48 : result._vals[i] = _vals[i] / a;
541 8 : return result;
542 0 : }
543 :
544 : /// Defines multiplication with a vector to get a vector
545 : template <typename T>
546 : template <typename T2>
547 : TypeVector<typename libMesh::CompareTypes<T, T2>::supertype>
548 : SymmetricRankTwoTensorTempl<T>::operator*(const TypeVector<T2> & a) const
549 : {
550 : TypeVector<typename libMesh::CompareTypes<T, T2>::supertype> ret;
551 : ret(0) = a(0) * _vals[0] + a(1) * _vals[5] + a(2) * _vals[4];
552 : ret(1) = a(0) * _vals[5] + a(1) * _vals[1] + a(2) * _vals[3];
553 : ret(2) = a(0) * _vals[4] + a(1) * _vals[3] + a(2) * _vals[2];
554 : }
555 :
556 : template <typename T>
557 : template <typename T2>
558 : bool
559 0 : SymmetricRankTwoTensorTempl<T>::operator==(const SymmetricRankTwoTensorTempl<T2> & a) const
560 : {
561 0 : for (std::size_t i = 0; i < N; ++i)
562 0 : if (_vals[i] != a._vals[i])
563 0 : return false;
564 0 : return true;
565 : }
566 :
567 : template <typename T>
568 : template <typename T2>
569 : bool
570 : SymmetricRankTwoTensorTempl<T>::operator!=(const SymmetricRankTwoTensorTempl<T2> & a) const
571 : {
572 : for (std::size_t i = 0; i < N; ++i)
573 : if (_vals[i] != a._vals[i])
574 : return true;
575 : return false;
576 : }
577 :
578 : template <typename T>
579 : SymmetricRankFourTensorTempl<T>
580 2 : SymmetricRankTwoTensorTempl<T>::positiveProjectionEigenDecomposition(
581 : std::vector<T> & eigval, RankTwoTensorTempl<T> & eigvec) const
582 : {
583 : using std::abs;
584 :
585 : // The calculate of projection tensor follows
586 : // C. Miehe and M. Lambrecht, Commun. Numer. Meth. Engng 2001; 17:337~353
587 :
588 : // Compute eigenvectors and eigenvalues of this tensor
589 : if (MooseUtils::IsLikeReal<T>::value)
590 : {
591 2 : this->symmetricEigenvaluesEigenvectors(eigval, eigvec);
592 :
593 : // Separate out positive and negative eigen values
594 0 : std::array<T, N> epos;
595 0 : std::array<T, N> d;
596 8 : for (unsigned int i = 0; i < Ndim; ++i)
597 : {
598 6 : epos[i] = (abs(eigval[i]) + eigval[i]) / 2.0;
599 6 : d[i] = 0 < eigval[i] ? 1.0 : 0.0;
600 : }
601 :
602 : // projection tensor
603 2 : SymmetricRankFourTensorTempl<T> proj_pos;
604 :
605 8 : for (unsigned int a = 0; a < Ndim; ++a)
606 : {
607 6 : const auto Ma = SymmetricRankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(a));
608 6 : proj_pos += d[a] * Ma.outerProduct(Ma);
609 : }
610 :
611 8 : for (const auto a : make_range(Ndim))
612 12 : for (const auto b : make_range(a))
613 : {
614 6 : const auto Ma = SymmetricRankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(a));
615 6 : const auto Mb = SymmetricRankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(b));
616 :
617 6 : SymmetricRankFourTensorTempl<T> Gabba(SymmetricRankFourTensorTempl<T>::initNone);
618 42 : for (const auto aa : make_range(N))
619 252 : for (const auto bb : make_range(N))
620 : {
621 216 : const auto i = SymmetricRankFourTensorTempl<T>::full_index[aa][bb][0];
622 216 : const auto j = SymmetricRankFourTensorTempl<T>::full_index[aa][bb][1];
623 216 : const auto k = SymmetricRankFourTensorTempl<T>::full_index[aa][bb][2];
624 216 : const auto l = SymmetricRankFourTensorTempl<T>::full_index[aa][bb][3];
625 :
626 432 : Gabba(aa, bb) = (Ma(i, k) * Mb(j, l) + Ma(i, l) * Mb(j, k) + Ma(j, l) * Mb(i, k) +
627 216 : Ma(j, k) * Mb(i, l)) *
628 216 : SymmetricRankFourTensorTempl<T>::mandelFactor(aa, bb);
629 : }
630 :
631 0 : T theta_ab;
632 6 : if (!MooseUtils::relativeFuzzyEqual(eigval[a], eigval[b]))
633 6 : theta_ab = 0.5 * (epos[a] - epos[b]) / (eigval[a] - eigval[b]);
634 : else
635 0 : theta_ab = 0.25 * (d[a] + d[b]);
636 :
637 6 : proj_pos += theta_ab * Gabba;
638 : }
639 2 : return proj_pos;
640 0 : }
641 : else
642 : mooseError("positiveProjectionEigenDecomposition is only available for ordered tensor "
643 : "component types");
644 : }
645 :
646 : template <typename T>
647 : T
648 4 : SymmetricRankTwoTensorTempl<T>::sin3Lode(const T & r0, const T & r0_value) const
649 : {
650 : using std::max, std::min, std::sqrt, std::pow;
651 : if (MooseUtils::IsLikeReal<T>::value)
652 : {
653 4 : T bar = secondInvariant();
654 4 : if (bar <= r0)
655 : // in this case the Lode angle is not defined
656 2 : return r0_value;
657 : else
658 : // the min and max here gaurd against precision-loss when bar is tiny but nonzero.
659 2 : return max(min(thirdInvariant() * -1.5 * sqrt(3.0) / pow(bar, 1.5), 1.0), -1.0);
660 0 : }
661 : else
662 : mooseError("sin3Lode is only available for ordered tensor component types");
663 : }
664 :
665 : template <typename T>
666 : SymmetricRankTwoTensorTempl<T>
667 2 : SymmetricRankTwoTensorTempl<T>::dsin3Lode(const T & r0) const
668 : {
669 : using std::sqrt, std::pow;
670 : if (MooseUtils::IsLikeReal<T>::value)
671 : {
672 2 : T bar = secondInvariant();
673 2 : if (bar <= r0)
674 0 : return SymmetricRankTwoTensorTempl<T>();
675 : else
676 2 : return -1.5 * sqrt(3.0) *
677 2 : (dthirdInvariant() / pow(bar, 1.5) -
678 4 : 1.5 * dsecondInvariant() * thirdInvariant() / pow(bar, 2.5));
679 0 : }
680 : else
681 : mooseError("dsin3Lode is only available for ordered tensor component types");
682 : }
683 :
684 : template <typename T>
685 : template <typename T2>
686 : SymmetricRankTwoTensorTempl<T> &
687 35 : SymmetricRankTwoTensorTempl<T>::operator=(const SymmetricRankTwoTensorTempl<T2> & a)
688 : {
689 245 : for (const auto i : make_range(N))
690 210 : (*this)(i) = a(i);
691 35 : return *this;
692 : }
693 :
694 : // non-member operators
695 :
696 : template <typename T, typename Scalar>
697 : inline typename std::enable_if_t<
698 : libMesh::ScalarTraits<Scalar>::value,
699 : SymmetricRankTwoTensorTempl<typename libMesh::CompareTypes<T, Scalar>::supertype>>
700 4 : operator*(const Scalar & factor, const SymmetricRankTwoTensorTempl<T> & t)
701 : {
702 4 : return t * factor;
703 : }
|
