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 "ADRankFourTensorForward.h"
15 : #include "ADSymmetricRankTwoTensorForward.h"
16 : #include "ADSymmetricRankFourTensorForward.h"
17 : #include "MooseUtils.h"
18 : #include "MathUtils.h"
19 :
20 : #include "libmesh/libmesh.h"
21 : #include "libmesh/tuple_of.h"
22 : #include "libmesh/int_range.h"
23 :
24 : #include "metaphysicl/raw_type.h"
25 :
26 : #include <petscsys.h>
27 :
28 : // Eigen includes
29 : #include <Eigen/Core>
30 : #include <Eigen/Dense>
31 :
32 : #include <array>
33 :
34 : using libMesh::Real;
35 : using libMesh::tuple_of;
36 : namespace libMesh
37 : {
38 : template <typename>
39 : class TensorValue;
40 : template <typename>
41 : class TypeTensor;
42 : template <typename>
43 : class VectorValue;
44 : }
45 :
46 : // Forward declarations
47 : class MooseEnum;
48 :
49 : namespace MathUtils
50 : {
51 : /**
52 : * Helper function template specialization to set an object to zero.
53 : * Needed by DerivativeMaterialInterface
54 : */
55 : template <>
56 : void mooseSetToZero<SymmetricRankFourTensor>(SymmetricRankFourTensor & v);
57 : template <>
58 : void mooseSetToZero<ADSymmetricRankFourTensor>(ADSymmetricRankFourTensor & v);
59 : }
60 :
61 : /**
62 : * SymmetricRankFourTensorTempl is designed to handle an N-dimensional fourth order tensor with
63 : * minor symmetry, C. Since N is hard-coded to 3, SymmetricRankFourTensorTempl holds 36 separate
64 : * C_ij entries. Within the code i,j = 0, .., 5.
65 : */
66 : template <typename T>
67 : class SymmetricRankFourTensorTempl
68 : {
69 : public:
70 : ///@{ tensor dimension, Mandel matrix dimension, and Mandel matrix size
71 : static constexpr unsigned int Ndim = LIBMESH_DIM;
72 : static constexpr unsigned int N = Ndim + Ndim * (Ndim - 1) / 2;
73 : static constexpr unsigned int N2 = N * N;
74 : ///@}
75 :
76 : // Full tensor indices in the Mandel/Voigt representation
77 : static constexpr unsigned int full_index[6][6][4] = {
78 : {{0, 0, 0, 0}, {0, 0, 1, 1}, {0, 0, 2, 2}, {0, 0, 1, 2}, {0, 0, 0, 2}, {0, 0, 0, 1}},
79 : {{1, 1, 0, 0}, {1, 1, 1, 1}, {1, 1, 2, 2}, {1, 1, 1, 2}, {1, 1, 0, 2}, {1, 1, 0, 1}},
80 : {{2, 2, 0, 0}, {2, 2, 1, 1}, {2, 2, 2, 2}, {2, 2, 1, 2}, {2, 2, 0, 2}, {2, 2, 0, 1}},
81 : {{1, 2, 0, 0}, {1, 2, 1, 1}, {1, 2, 2, 2}, {1, 2, 1, 2}, {1, 2, 0, 2}, {1, 2, 0, 1}},
82 : {{0, 2, 0, 0}, {0, 2, 1, 1}, {0, 2, 2, 2}, {0, 2, 1, 2}, {0, 2, 0, 2}, {0, 2, 0, 1}},
83 : {{0, 1, 0, 0}, {0, 1, 1, 1}, {0, 1, 2, 2}, {0, 1, 1, 2}, {0, 1, 0, 2}, {0, 1, 0, 1}}};
84 :
85 : /// returns the 1, sqrt(2), or 2 prefactor in the Mandel notation for the indices i,j ranging from 0-5.
86 12384 : static constexpr Real mandelFactor(unsigned int i, unsigned int j)
87 : {
88 12384 : return SymmetricRankTwoTensorTempl<T>::mandelFactor(i) *
89 12384 : SymmetricRankTwoTensorTempl<T>::mandelFactor(j);
90 : }
91 :
92 : /// Initialization method
93 : enum InitMethod
94 : {
95 : initNone,
96 : initIdentity,
97 : initIdentitySymmetricFour
98 : };
99 :
100 : /**
101 : * To fill up the 36 entries in the 4th-order tensor, fillFromInputVector
102 : * is called with one of the following fill_methods.
103 : * See the fill*FromInputVector functions for more details
104 : */
105 : enum FillMethod
106 : {
107 : symmetric9, // fillSymmetric9FromInputVector
108 : symmetric21, // fillSymmetric21FromInputVector
109 : symmetric_isotropic, // fillSymmetricIsotropicFromInputVector
110 : symmetric_isotropic_E_nu, // fillSymmetricIsotropicEandNu
111 : axisymmetric_rz, // fillAxisymmetricRZFromInputVector
112 : principal, // fillPrincipalFromInputVector
113 : orthotropic // fillGeneralOrthotropicFromInputVector
114 : };
115 :
116 : template <template <typename> class Tensor, typename Scalar>
117 : struct TwoTensorMultTraits
118 : {
119 : static const bool value = false;
120 : };
121 : template <typename Scalar>
122 : struct TwoTensorMultTraits<SymmetricRankTwoTensorTempl, Scalar>
123 : {
124 : static const bool value = libMesh::ScalarTraits<Scalar>::value;
125 : };
126 : template <typename Scalar>
127 : struct TwoTensorMultTraits<TensorValue, Scalar>
128 : {
129 : static const bool value = libMesh::ScalarTraits<Scalar>::value;
130 : };
131 : template <typename Scalar>
132 : struct TwoTensorMultTraits<TypeTensor, Scalar>
133 : {
134 : static const bool value = libMesh::ScalarTraits<Scalar>::value;
135 : };
136 :
137 : /// Default constructor; fills to zero
138 : SymmetricRankFourTensorTempl();
139 :
140 : /// Select specific initialization pattern
141 : SymmetricRankFourTensorTempl(const InitMethod);
142 :
143 : /// Fill from vector
144 : SymmetricRankFourTensorTempl(const std::vector<T> &, FillMethod);
145 :
146 : /// Copy assignment operator must be defined if used
147 0 : SymmetricRankFourTensorTempl(const SymmetricRankFourTensorTempl<T> & a) = default;
148 :
149 : /**
150 : * Copy constructor
151 : */
152 : template <typename T2>
153 : SymmetricRankFourTensorTempl(const SymmetricRankFourTensorTempl<T2> & copy);
154 :
155 : /// Copy constructor from RankFourTensorTempl<T>
156 : explicit SymmetricRankFourTensorTempl(const RankFourTensorTempl<T> & a);
157 :
158 : /// The conversion operator to `RankFourTensorTempl`
159 : explicit operator RankFourTensorTempl<T>();
160 :
161 : // Named constructors
162 0 : static SymmetricRankFourTensorTempl<T> identity()
163 : {
164 0 : return SymmetricRankFourTensorTempl<T>(initIdentity);
165 : }
166 0 : static SymmetricRankFourTensorTempl<T> identitySymmetricFour()
167 : {
168 0 : return SymmetricRankFourTensorTempl<T>(initIdentitySymmetricFour);
169 : };
170 :
171 : /// Gets the value for the indices specified. Takes indices ranging from 0-5 for i and j.
172 14189928 : inline T & operator()(unsigned int i, unsigned int j) { return _vals[i * N + j]; }
173 :
174 : /**
175 : * Gets the value for the indices specified. Takes indices ranging from 0-5 for i and j.
176 : * used for const
177 : */
178 3042 : inline const T & operator()(unsigned int i, unsigned int j) const { return _vals[i * N + j]; }
179 :
180 : /// Zeros out the tensor.
181 : void zero();
182 :
183 : /// Print the rank four tensor
184 : void print(std::ostream & stm = Moose::out) const;
185 :
186 : /// Print the values of the rank four tensor
187 : void printReal(std::ostream & stm = Moose::out) const;
188 :
189 : friend std::ostream & operator<<(std::ostream & os, const SymmetricRankFourTensorTempl<T> & t)
190 : {
191 : t.print(os);
192 : return os;
193 : }
194 :
195 : /// copies values from a into this tensor
196 360002 : SymmetricRankFourTensorTempl<T> & operator=(const SymmetricRankFourTensorTempl<T> & a) = default;
197 :
198 : /**
199 : * Assignment-from-scalar operator. Used only to zero out the tensor.
200 : *
201 : * \returns A reference to *this.
202 : */
203 : template <typename Scalar>
204 : typename std::enable_if<libMesh::ScalarTraits<Scalar>::value,
205 : SymmetricRankFourTensorTempl &>::type
206 : operator=(const Scalar & libmesh_dbg_var(p))
207 : {
208 : libmesh_assert_equal_to(p, Scalar(0));
209 : this->zero();
210 : return *this;
211 : }
212 :
213 : /// C_ijkl*a_kl
214 : template <typename T2>
215 : auto operator*(const SymmetricRankTwoTensorTempl<T2> & b) const
216 : -> SymmetricRankTwoTensorTempl<decltype(T() * T2())>;
217 :
218 : /// C_ijkl*a
219 : template <typename T2>
220 : auto operator*(const T2 & a) const ->
221 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
222 : SymmetricRankFourTensorTempl<decltype(T() * T2())>>::type;
223 :
224 : /// C_ijkl *= a
225 : SymmetricRankFourTensorTempl<T> & operator*=(const T & a);
226 :
227 : /// C_ijkl/a
228 : template <typename T2>
229 : auto operator/(const T2 & a) const ->
230 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
231 : SymmetricRankFourTensorTempl<decltype(T() / T2())>>::type;
232 :
233 : /// C_ijkl /= a for all i, j, k, l
234 : SymmetricRankFourTensorTempl<T> & operator/=(const T & a);
235 :
236 : /// C_ijkl += a_ijkl for all i, j, k, l
237 : SymmetricRankFourTensorTempl<T> & operator+=(const SymmetricRankFourTensorTempl<T> & a);
238 :
239 : /// C_ijkl + a_ijkl
240 : template <typename T2>
241 : auto operator+(const SymmetricRankFourTensorTempl<T2> & a) const
242 : -> SymmetricRankFourTensorTempl<decltype(T() + T2())>;
243 :
244 : /// C_ijkl -= a_ijkl
245 : SymmetricRankFourTensorTempl<T> & operator-=(const SymmetricRankFourTensorTempl<T> & a);
246 :
247 : /// C_ijkl - a_ijkl
248 : template <typename T2>
249 : auto operator-(const SymmetricRankFourTensorTempl<T2> & a) const
250 : -> SymmetricRankFourTensorTempl<decltype(T() - T2())>;
251 :
252 : /// -C_ijkl
253 : SymmetricRankFourTensorTempl<T> operator-() const;
254 :
255 : /// C_ijpq*a_pqkl
256 : template <typename T2>
257 : auto operator*(const SymmetricRankFourTensorTempl<T2> & a) const
258 : -> SymmetricRankFourTensorTempl<decltype(T() * T2())>;
259 :
260 : /// sqrt(C_ijkl*C_ijkl)
261 : T L2norm() const;
262 :
263 : /**
264 : * This returns A_ijkl such that C_ijkl*A_klmn = 0.5*(de_im de_jn + de_in de_jm)
265 : * This routine assumes that C_ijkl = C_jikl = C_ijlk
266 : */
267 : SymmetricRankFourTensorTempl<T> invSymm() const;
268 :
269 : /**
270 : * Build a 6x6 rotation matrix
271 : * MEHRABADI, MORTEZA M.; COWIN, STEPHEN C. (1990). EIGENTENSORS OF LINEAR ANISOTROPIC ELASTIC
272 : * MATERIALS. The Quarterly Journal of Mechanics and Applied Mathematics, 43(1), 15-41.
273 : * doi:10.1093/qjmam/43.1.15
274 : */
275 : static SymmetricRankFourTensorTempl<T> rotationMatrix(const TypeTensor<T> & R);
276 :
277 : /**
278 : * Rotate the tensor using
279 : * C_ijkl = R_im R_jn R_ko R_lp C_mnop
280 : */
281 : void rotate(const TypeTensor<T> & R);
282 :
283 : /**
284 : * Transpose the tensor by swapping the first pair with the second pair of indices
285 : * This amounts to a regular transpose of the 6x6 matrix
286 : * @return C_klji
287 : */
288 : SymmetricRankFourTensorTempl<T> transposeMajor() const;
289 :
290 : /**
291 : * Transpose the tensor by swapping the first two indices - a no-op
292 : */
293 0 : SymmetricRankFourTensorTempl<T> transposeIj() const { return *this; }
294 :
295 : /// Static method for use in validParams for getting the "fill_method"
296 : static MooseEnum fillMethodEnum();
297 :
298 : /**
299 : * fillFromInputVector takes some number of inputs to fill
300 : * the Rank-4 tensor.
301 : * @param input the numbers that will be placed in the tensor
302 : * @param fill_method See FillMethod
303 : */
304 : void fillFromInputVector(const std::vector<T> & input, FillMethod fill_method);
305 :
306 : ///@{ Vector-less fill API functions. See docs of the corresponding ...FromInputVector methods
307 : void fillSymmetricIsotropic(const T & i0, const T & i1);
308 : void fillSymmetricIsotropicEandNu(const T & E, const T & nu);
309 : ///@}
310 :
311 : /**
312 : * fillSymmetric9FromInputVector takes 9 inputs to fill in
313 : * the Rank-4 tensor with the appropriate crystal symmetries maintained. I.e., C_ijkl = C_klij,
314 : * C_ijkl = C_ijlk, C_ijkl = C_jikl
315 : * @param input is:
316 : * C1111 C1122 C1133 C2222 C2233 C3333 C2323 C1313 C1212
317 : * In the isotropic case this is (la is first Lame constant, mu is second (shear)
318 : * Lame constant)
319 : * la+2mu la la la+2mu la la+2mu mu mu mu
320 : */
321 : template <typename T2>
322 : void fillSymmetric9FromInputVector(const T2 & input);
323 :
324 : /**
325 : * fillSymmetric21FromInputVector takes 21 inputs to fill in the Rank-4 tensor with the
326 : * appropriate crystal symmetries maintained.
327 : * I.e., C_ijkl = C_klij, C_ijkl = C_ijlk, C_ijkl = C_jikl
328 : *
329 : * @param input is C1111 C1122 C1133 C1123 C1113 C1112 C2222 C2233 C2223 C2213 C2212 C3333 C3323
330 : * C3313 C3312 C2323 C2313 C2312 C1313 C1312 C1212
331 : */
332 : template <typename T2>
333 : void fillSymmetric21FromInputVector(const T2 & input);
334 :
335 : /// Calculates the sum of Ciijj for i and j varying from 0 to 2
336 : T sum3x3() const;
337 :
338 : /// Calculates the vector a[i] = sum over j Ciijj for i and j varying from 0 to 2
339 : libMesh::VectorValue<T> sum3x1() const;
340 :
341 : /// checks if the tensor is symmetric
342 : bool isSymmetric() const;
343 :
344 : /// checks if the tensor is isotropic
345 : bool isIsotropic() const;
346 :
347 : protected:
348 : /// The values of the rank-four tensor
349 : std::array<T, N2> _vals;
350 :
351 : /**
352 : * fillSymmetricIsotropicFromInputVector takes 2 inputs to fill the
353 : * the symmetric Rank-4 tensor with the appropriate symmetries maintained.
354 : * C_ijkl = lambda*de_ij*de_kl + mu*(de_ik*de_jl + de_il*de_jk)
355 : * where lambda is the first Lame modulus, mu is the second (shear) Lame modulus,
356 : * @param input this is lambda and mu in the above formula
357 : */
358 : void fillSymmetricIsotropicFromInputVector(const std::vector<T> & input);
359 :
360 : /**
361 : * fillSymmetricIsotropicEandNuFromInputVector is a variation of the
362 : * fillSymmetricIsotropicFromInputVector which takes as inputs the
363 : * more commonly used Young's modulus (E) and Poisson's ratio (nu)
364 : * constants to fill the isotropic elasticity tensor. Using well-known formulas,
365 : * E and nu are used to calculate lambda and mu and then the vector is passed
366 : * to fillSymmetricIsotropicFromInputVector.
367 : * @param input Young's modulus (E) and Poisson's ratio (nu)
368 : */
369 : void fillSymmetricIsotropicEandNuFromInputVector(const std::vector<T> & input);
370 :
371 : /**
372 : * fillAxisymmetricRZFromInputVector takes 5 inputs to fill the axisymmetric
373 : * Rank-4 tensor with the appropriate symmetries maintatined for use with
374 : * axisymmetric problems using coord_type = RZ.
375 : * I.e. C1111 = C2222, C1133 = C2233, C2323 = C3131 and C1212 = 0.5*(C1111-C1122)
376 : * @param input this is C1111, C1122, C1133, C3333, C2323.
377 : */
378 : void fillAxisymmetricRZFromInputVector(const std::vector<T> & input);
379 :
380 : /**
381 : * fillPrincipalFromInputVector takes 9 inputs to fill a Rank-4 tensor
382 : * C1111 = input0
383 : * C1122 = input1
384 : * C1133 = input2
385 : * C2211 = input3
386 : * C2222 = input4
387 : * C2233 = input5
388 : * C3311 = input6
389 : * C3322 = input7
390 : * C3333 = input8
391 : * with all other components being zero
392 : */
393 :
394 : void fillPrincipalFromInputVector(const std::vector<T> & input);
395 :
396 : /**
397 : * fillGeneralOrhotropicFromInputVector takes 10 inputs to fill the Rank-4 tensor
398 : * It defines a general orthotropic tensor for which some constraints among
399 : * elastic parameters exist
400 : * @param input Ea, Eb, Ec, Gab, Gbc, Gca, nuba, nuca, nucb, nuab, nuac, nubc
401 : */
402 : void fillGeneralOrthotropicFromInputVector(const std::vector<T> & input);
403 :
404 : template <class T2>
405 : friend void dataStore(std::ostream &, SymmetricRankFourTensorTempl<T2> &, void *);
406 :
407 : template <class T2>
408 : friend void dataLoad(std::istream &, SymmetricRankFourTensorTempl<T2> &, void *);
409 :
410 : template <typename T2>
411 : friend class SymmetricRankTwoTensorTempl;
412 : template <typename T2>
413 : friend class SymmetricRankFourTensorTempl;
414 : template <typename T2>
415 : friend class RankThreeTensorTempl;
416 : };
417 :
418 : namespace MetaPhysicL
419 : {
420 : template <typename T>
421 : struct RawType<SymmetricRankFourTensorTempl<T>>
422 : {
423 : typedef SymmetricRankFourTensorTempl<typename RawType<T>::value_type> value_type;
424 :
425 8 : static value_type value(const SymmetricRankFourTensorTempl<T> & in)
426 : {
427 8 : value_type ret;
428 56 : for (unsigned int i = 0; i < SymmetricRankFourTensorTempl<T>::N; ++i)
429 336 : for (unsigned int j = 0; j < SymmetricRankFourTensorTempl<T>::N; ++j)
430 288 : ret(i, j) = raw_value(in(i, j));
431 :
432 8 : return ret;
433 : }
434 : };
435 : }
436 :
437 : template <typename T1, typename T2>
438 : inline auto
439 12 : operator*(const T1 & a, const SymmetricRankFourTensorTempl<T2> & b) ->
440 : typename std::enable_if<libMesh::ScalarTraits<T1>::value,
441 : SymmetricRankFourTensorTempl<decltype(T1() * T2())>>::type
442 : {
443 12 : return b * a;
444 : }
445 :
446 : template <typename T>
447 : template <typename T2>
448 6 : SymmetricRankFourTensorTempl<T>::SymmetricRankFourTensorTempl(
449 6 : const SymmetricRankFourTensorTempl<T2> & copy)
450 : {
451 222 : for (const auto i : make_range(N2))
452 216 : _vals[i] = copy._vals[i];
453 6 : }
454 :
455 : template <typename T>
456 : template <typename T2>
457 : auto
458 14 : SymmetricRankFourTensorTempl<T>::operator*(const T2 & b) const ->
459 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
460 : SymmetricRankFourTensorTempl<decltype(T() * T2())>>::type
461 : {
462 : typedef decltype(T() * T2()) ValueType;
463 14 : SymmetricRankFourTensorTempl<ValueType> result;
464 :
465 518 : for (const auto i : make_range(N2))
466 504 : result._vals[i] = _vals[i] * b;
467 :
468 14 : return result;
469 0 : }
470 :
471 : template <typename T>
472 : template <typename T2>
473 : auto
474 4 : SymmetricRankFourTensorTempl<T>::operator*(const SymmetricRankTwoTensorTempl<T2> & b) const
475 : -> SymmetricRankTwoTensorTempl<decltype(T() * T2())>
476 : {
477 : typedef decltype(T() * T2()) ValueType;
478 4 : SymmetricRankTwoTensorTempl<ValueType> result;
479 :
480 4 : std::size_t index = 0;
481 28 : for (const auto i : make_range(N))
482 : {
483 24 : ValueType tmp = 0.0;
484 168 : for (const auto j : make_range(N))
485 144 : tmp += _vals[index++] * b._vals[j];
486 24 : result._vals[i] = tmp;
487 : }
488 :
489 4 : return result;
490 0 : }
491 :
492 : template <typename T>
493 : template <typename T2>
494 : auto
495 2 : SymmetricRankFourTensorTempl<T>::operator/(const T2 & b) const ->
496 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
497 : SymmetricRankFourTensorTempl<decltype(T() / T2())>>::type
498 : {
499 2 : SymmetricRankFourTensorTempl<decltype(T() / T2())> result;
500 74 : for (const auto i : make_range(N2))
501 72 : result._vals[i] = _vals[i] / b;
502 2 : return result;
503 : }
504 :
505 : template <typename T>
506 : template <typename T2>
507 : void
508 42 : SymmetricRankFourTensorTempl<T>::fillSymmetric9FromInputVector(const T2 & input)
509 : {
510 : mooseAssert(LIBMESH_DIM == 3, "This method assumes LIBMESH_DIM == 3");
511 : mooseAssert(input.size() == 9,
512 : "To use fillSymmetric9FromInputVector, your input must have size 9.");
513 42 : zero();
514 :
515 42 : _vals[0] = input[0]; // C1111
516 42 : _vals[7] = input[3]; // C2222
517 42 : _vals[14] = input[5]; // C3333
518 :
519 42 : _vals[1] = _vals[6] = input[1]; // C1122
520 42 : _vals[2] = _vals[12] = input[2]; // C1133
521 42 : _vals[8] = _vals[13] = input[4]; // C2233
522 :
523 : static constexpr std::size_t C2323 = 21;
524 : static constexpr std::size_t C1313 = 28;
525 : static constexpr std::size_t C1212 = 35;
526 :
527 42 : _vals[C2323] = 2.0 * input[6];
528 42 : _vals[C1313] = 2.0 * input[7];
529 42 : _vals[C1212] = 2.0 * input[8];
530 42 : }
531 :
532 : template <typename T>
533 : template <typename T2>
534 : void
535 48 : SymmetricRankFourTensorTempl<T>::fillSymmetric21FromInputVector(const T2 & input)
536 : {
537 : // C1111 C1122 C1133 C1123 C1113 C1112
538 : // C2222 C2233 C2223 C2213 C2212
539 : // C3333 C3323 C3313 C3312
540 : // C2323 C2313 C2312
541 : // C1313 C1312
542 : // C1212
543 :
544 : mooseAssert(LIBMESH_DIM == 3, "This method assumes LIBMESH_DIM == 3");
545 : mooseAssert(input.size() == 21,
546 : "To use fillSymmetric21FromInputVector, your input must have size 21.");
547 48 : std::size_t index = 0;
548 336 : for (const auto i : make_range(N))
549 1296 : for (const auto j : make_range(i, N))
550 : {
551 1008 : _vals[i + N * j] = mandelFactor(i, j) * input[index];
552 1008 : _vals[j + N * i] = mandelFactor(j, i) * input[index];
553 1008 : index++;
554 : }
555 48 : }
556 :
557 : template <typename T>
558 : SymmetricRankFourTensorTempl<T>
559 8 : SymmetricRankFourTensorTempl<T>::invSymm() const
560 : {
561 8 : SymmetricRankFourTensorTempl<T> result(initNone);
562 :
563 : if constexpr (SymmetricRankFourTensorTempl<T>::N2 * sizeof(T) > EIGEN_STACK_ALLOCATION_LIMIT)
564 : {
565 2 : Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> mat(N, N);
566 14 : for (const auto i : make_range(N))
567 84 : for (const auto j : make_range(N))
568 72 : mat(i, j) = (*this)(i, j);
569 2 : mat = mat.inverse();
570 14 : for (const auto i : make_range(N))
571 84 : for (const auto j : make_range(N))
572 72 : result(i, j) = mat(i, j);
573 2 : }
574 : else
575 : {
576 6 : const Eigen::Map<const Eigen::Matrix<T, N, N, Eigen::RowMajor>> mat(&_vals[0]);
577 6 : Eigen::Map<Eigen::Matrix<T, N, N, Eigen::RowMajor>> res(&result._vals[0]);
578 6 : res = mat.inverse();
579 : }
580 :
581 8 : return result;
582 0 : }
|
