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