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