Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "Moose.h"
13 : #include "MooseTypes.h"
14 : #include "ADRankTwoTensorForward.h"
15 : #include "ADRankFourTensorForward.h"
16 : #include "ADRankThreeTensorForward.h"
17 : #include "MooseError.h"
18 :
19 : #include "libmesh/libmesh.h"
20 : #include "libmesh/tuple_of.h"
21 : #include "libmesh/int_range.h"
22 :
23 : #include "metaphysicl/raw_type.h"
24 :
25 : #include <petscsys.h>
26 :
27 : #include <Eigen/Core>
28 : #include <Eigen/Dense>
29 :
30 : using libMesh::tuple_of;
31 : namespace libMesh
32 : {
33 : template <typename>
34 : class TensorValue;
35 : template <typename>
36 : class TypeTensor;
37 : template <typename>
38 : class VectorValue;
39 : }
40 :
41 : // Forward declarations
42 : class MooseEnum;
43 :
44 : namespace MathUtils
45 : {
46 : template <typename T>
47 : void mooseSetToZero(T & v);
48 :
49 : /**
50 : * Helper function template specialization to set an object to zero.
51 : * Needed by DerivativeMaterialInterface
52 : */
53 : template <>
54 : void mooseSetToZero<RankFourTensor>(RankFourTensor & v);
55 : template <>
56 : void mooseSetToZero<ADRankFourTensor>(ADRankFourTensor & v);
57 : }
58 :
59 : /**
60 : * RankFourTensorTempl is designed to handle any N-dimensional fourth order tensor, C.
61 : * Since N is hard-coded to 3, RankFourTensorTempl holds 81 separate C_ijkl entries,
62 : * with i,j,k,l = 0, 1, 2.
63 : */
64 : template <typename T>
65 : class RankFourTensorTempl
66 : {
67 : public:
68 : ///@{ tensor dimension and powers of the dimension
69 : static constexpr unsigned int N = Moose::dim;
70 : static constexpr unsigned int N2 = N * N;
71 : static constexpr unsigned int N3 = N * N * N;
72 : static constexpr unsigned int N4 = N * N * N * N;
73 : ///@}
74 :
75 : typedef tuple_of<4, unsigned int> index_type;
76 : typedef T value_type;
77 :
78 : /// Initialization method
79 : enum InitMethod
80 : {
81 : initNone,
82 : initIdentity,
83 : initIdentityFour,
84 : initIdentitySymmetricFour,
85 : initIdentityDeviatoric
86 : };
87 :
88 : /**
89 : * To fill up the 81 entries in the 4th-order tensor, fillFromInputVector
90 : * is called with one of the following fill_methods.
91 : * See the fill*FromInputVector functions for more details
92 : */
93 : enum FillMethod
94 : {
95 : antisymmetric,
96 : symmetric9,
97 : symmetric21,
98 : general_isotropic,
99 : symmetric_isotropic,
100 : symmetric_isotropic_E_nu,
101 : antisymmetric_isotropic,
102 : axisymmetric_rz,
103 : general,
104 : principal,
105 : orthotropic
106 : };
107 :
108 : template <template <typename> class Tensor, typename Scalar>
109 : struct TwoTensorMultTraits
110 : {
111 : static const bool value = false;
112 : };
113 : template <typename Scalar>
114 : struct TwoTensorMultTraits<RankTwoTensorTempl, Scalar>
115 : {
116 : static const bool value = libMesh::ScalarTraits<Scalar>::value;
117 : };
118 : template <typename Scalar>
119 : struct TwoTensorMultTraits<TensorValue, Scalar>
120 : {
121 : static const bool value = libMesh::ScalarTraits<Scalar>::value;
122 : };
123 : template <typename Scalar>
124 : struct TwoTensorMultTraits<TypeTensor, Scalar>
125 : {
126 : static const bool value = libMesh::ScalarTraits<Scalar>::value;
127 : };
128 :
129 : /// Default constructor; fills to zero
130 : RankFourTensorTempl();
131 :
132 : /// Select specific initialization pattern
133 : RankFourTensorTempl(const InitMethod);
134 :
135 : /// Fill from vector
136 : RankFourTensorTempl(const std::vector<T> &, FillMethod);
137 :
138 : /// Copy assignment operator must be defined if used
139 0 : RankFourTensorTempl(const RankFourTensorTempl<T> & a) = default;
140 :
141 : /**
142 : * Copy constructor
143 : */
144 : template <typename T2>
145 : RankFourTensorTempl(const RankFourTensorTempl<T2> & copy);
146 :
147 : /**
148 : * The conversion operator from a `SymmetricRankFourTensorTempl`
149 : */
150 : template <typename T2>
151 : RankFourTensorTempl(const SymmetricRankFourTensorTempl<T2> & t)
152 : {
153 : for (const auto a : make_range(SymmetricRankFourTensorTempl<T2>::N))
154 : for (const auto b : make_range(SymmetricRankFourTensorTempl<T2>::N))
155 : {
156 : const auto & idx = SymmetricRankFourTensorTempl<T2>::full_index[a][b];
157 : const auto i = idx[0];
158 : const auto j = idx[1];
159 : const auto k = idx[2];
160 : const auto l = idx[3];
161 :
162 : // Rijkl = Rjikl = Rijlk = Rjilk
163 : (*this)(i, j, k, l) = t(a, b) / SymmetricRankFourTensorTempl<T2>::mandelFactor(a, b);
164 : (*this)(j, i, k, l) = t(a, b) / SymmetricRankFourTensorTempl<T2>::mandelFactor(a, b);
165 : (*this)(i, j, l, k) = t(a, b) / SymmetricRankFourTensorTempl<T2>::mandelFactor(a, b);
166 : (*this)(j, i, l, k) = t(a, b) / SymmetricRankFourTensorTempl<T2>::mandelFactor(a, b);
167 : }
168 : }
169 :
170 : // Named constructors
171 0 : static RankFourTensorTempl<T> Identity() { return RankFourTensorTempl<T>(initIdentity); }
172 2 : static RankFourTensorTempl<T> IdentityFour() { return RankFourTensorTempl<T>(initIdentityFour); };
173 : /// Identity of type \delta_{ik} \delta_{jl} - \delta_{ij} \delta_{kl} / 3
174 0 : static RankFourTensorTempl<T> IdentityDeviatoric()
175 : {
176 0 : return RankFourTensorTempl<T>(initIdentityDeviatoric);
177 : };
178 :
179 : /// Gets the value for the indices specified. Takes indices ranging from 0-2 for i, j, k, and l.
180 160790861 : inline T & operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
181 : {
182 160790861 : return _vals[i * N3 + j * N2 + k * N + l];
183 : }
184 :
185 : /**
186 : * Gets the value for the indices specified. Takes indices ranging from 0-2 for i, j, k, and l.
187 : * used for const
188 : */
189 1721668 : inline const T & operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
190 : {
191 1721668 : return _vals[i * N3 + j * N2 + k * N + l];
192 : }
193 :
194 : /// Zeros out the tensor.
195 : void zero();
196 :
197 : /// Print the rank four tensor
198 : void print(std::ostream & stm = Moose::out) const;
199 :
200 : friend std::ostream & operator<<(std::ostream & os, const RankFourTensorTempl<T> & t)
201 : {
202 : t.print(os);
203 : return os;
204 : }
205 :
206 : /// Print the values of the rank four tensor
207 : void printReal(std::ostream & stm = Moose::out) const;
208 :
209 : /// copies values from a into this tensor
210 66420492 : RankFourTensorTempl<T> & operator=(const RankFourTensorTempl<T> & a) = default;
211 :
212 : /**
213 : * Assignment-from-scalar operator. Used only to zero out the tensor.
214 : *
215 : * \returns A reference to *this.
216 : */
217 : template <typename Scalar>
218 : typename std::enable_if<libMesh::ScalarTraits<Scalar>::value, RankFourTensorTempl &>::type
219 : operator=(const Scalar & libmesh_dbg_var(p))
220 : {
221 : libmesh_assert_equal_to(p, Scalar(0));
222 : this->zero();
223 : return *this;
224 : }
225 :
226 : /// C_ijkl*a_kl
227 : template <template <typename> class Tensor, typename T2>
228 : auto operator*(const Tensor<T2> & a) const ->
229 : typename std::enable_if<TwoTensorMultTraits<Tensor, T2>::value,
230 : RankTwoTensorTempl<decltype(T() * T2())>>::type;
231 :
232 : /// C_ijkl*a
233 : template <typename T2>
234 : auto operator*(const T2 & a) const ->
235 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
236 : RankFourTensorTempl<decltype(T() * T2())>>::type;
237 :
238 : /// C_ijkl *= a
239 : RankFourTensorTempl<T> & operator*=(const T & a);
240 :
241 : /// C_ijkl/a
242 : template <typename T2>
243 : auto operator/(const T2 & a) const ->
244 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
245 : RankFourTensorTempl<decltype(T() / T2())>>::type;
246 :
247 : /// C_ijkl /= a for all i, j, k, l
248 : RankFourTensorTempl<T> & operator/=(const T & a);
249 :
250 : /// C_ijkl += a_ijkl for all i, j, k, l
251 : RankFourTensorTempl<T> & operator+=(const RankFourTensorTempl<T> & a);
252 :
253 : /// C_ijkl + a_ijkl
254 : template <typename T2>
255 : auto operator+(const RankFourTensorTempl<T2> & a) const
256 : -> RankFourTensorTempl<decltype(T() + T2())>;
257 :
258 : /// C_ijkl -= a_ijkl
259 : RankFourTensorTempl<T> & operator-=(const RankFourTensorTempl<T> & a);
260 :
261 : /// C_ijkl - a_ijkl
262 : template <typename T2>
263 : auto operator-(const RankFourTensorTempl<T2> & a) const
264 : -> RankFourTensorTempl<decltype(T() - T2())>;
265 :
266 : /// -C_ijkl
267 : RankFourTensorTempl<T> operator-() const;
268 :
269 : /// C_ijpq*a_pqkl
270 : template <typename T2>
271 : auto operator*(const RankFourTensorTempl<T2> & a) const
272 : -> RankFourTensorTempl<decltype(T() * T2())>;
273 :
274 : /// sqrt(C_ijkl*C_ijkl)
275 : T L2norm() const;
276 :
277 : /**
278 : * This returns A_ijkl such that C_ijkl*A_klmn = 0.5*(de_im de_jn + de_in de_jm)
279 : * This routine assumes that C_ijkl = C_jikl = C_ijlk
280 : */
281 : RankFourTensorTempl<T> invSymm() const;
282 :
283 : /**
284 : * This returns A_ijkl such that C_ijkl*A_klmn = de_im de_jn
285 : * i.e. the general rank four inverse
286 : */
287 : RankFourTensorTempl<T> inverse() const;
288 :
289 : /**
290 : * Rotate the tensor using
291 : * C_ijkl = R_im R_jn R_ko R_lp C_mnop
292 : */
293 : void rotate(const TypeTensor<T> & R);
294 :
295 : /**
296 : * Transpose the tensor by swapping the first pair with the second pair of indices
297 : * @return C_klji
298 : */
299 : RankFourTensorTempl<T> transposeMajor() const;
300 :
301 : /**
302 : * Transpose the tensor by swapping the first two indeces
303 : * @return C_jikl
304 : */
305 : RankFourTensorTempl<T> transposeIj() const;
306 :
307 : /**
308 : * Transpose the tensor by swapping the last two indeces
309 : * @return C_ijlk
310 : */
311 : RankFourTensorTempl<T> transposeKl() const;
312 :
313 : /**
314 : * single contraction of a RankFourTensor with a vector over index m
315 : * @return C_xxx = a_ijkl*b_m where m={i,j,k,l} and xxx the remaining indices
316 : */
317 : template <int m>
318 : RankThreeTensorTempl<T> contraction(const libMesh::VectorValue<T> & b) const;
319 :
320 : /**
321 : * Fills the tensor entries ignoring the last dimension (ie, C_ijkl=0 if any of i, j, k, or l =
322 : * 3).
323 : * Fill method depends on size of input
324 : * Input size = 2. Then C_1111 = C_2222 = input[0], and C_1122 = input[1], and C_1212 = (input[0]
325 : * - input[1])/2,
326 : * and C_ijkl = C_jikl = C_ijlk = C_klij, and C_1211 = C_1222 = 0.
327 : * Input size = 9. Then C_1111 = input[0], C_1112 = input[1], C_1122 = input[3],
328 : * C_1212 = input[4], C_1222 = input[5], C_1211 = input[6]
329 : * C_2211 = input[7], C_2212 = input[8], C_2222 = input[9]
330 : * and C_ijkl = C_jikl = C_ijlk
331 : */
332 : void surfaceFillFromInputVector(const std::vector<T> & input);
333 :
334 : /// Static method for use in validParams for getting the "fill_method"
335 : static MooseEnum fillMethodEnum();
336 :
337 : /**
338 : * fillFromInputVector takes some number of inputs to fill
339 : * the Rank-4 tensor.
340 : * @param input the numbers that will be placed in the tensor
341 : * @param fill_method this can be:
342 : * antisymmetric (use fillAntisymmetricFromInputVector)
343 : * symmetric9 (use fillSymmetric9FromInputVector)
344 : * symmetric21 (use fillSymmetric21FromInputVector)
345 : * general_isotropic (use fillGeneralIsotropicFrominputVector)
346 : * symmetric_isotropic (use fillSymmetricIsotropicFromInputVector)
347 : * antisymmetric_isotropic (use fillAntisymmetricIsotropicFromInputVector)
348 : * axisymmetric_rz (use fillAxisymmetricRZFromInputVector)
349 : * general (use fillGeneralFromInputVector)
350 : * principal (use fillPrincipalFromInputVector)
351 : */
352 : void fillFromInputVector(const std::vector<T> & input, FillMethod fill_method);
353 :
354 : ///@{ Vector-less fill API functions. See docs of the corresponding ...FromInputVector methods
355 : void fillGeneralIsotropic(const T & i0, const T & i1, const T & i2);
356 : void fillAntisymmetricIsotropic(const T & i0);
357 : void fillSymmetricIsotropic(const T & i0, const T & i1);
358 : void fillSymmetricIsotropicEandNu(const T & E, const T & nu);
359 : ///@}
360 :
361 : /**
362 : * fillSymmetric9FromInputVector takes 9 inputs to fill in
363 : * the Rank-4 tensor with the appropriate crystal symmetries maintained. I.e., C_ijkl = C_klij,
364 : * C_ijkl = C_ijlk, C_ijkl = C_jikl
365 : * @param input is:
366 : * C1111 C1122 C1133 C2222 C2233 C3333 C2323 C1313 C1212
367 : * In the isotropic case this is (la is first Lame constant, mu is second (shear)
368 : * Lame constant)
369 : * la+2mu la la la+2mu la la+2mu mu mu mu
370 : */
371 : template <typename T2>
372 : void fillSymmetric9FromInputVector(const T2 & input);
373 :
374 : /**
375 : * fillSymmetric21FromInputVector takes either 21 inputs to fill in
376 : * the Rank-4 tensor with the appropriate crystal symmetries maintained. I.e., C_ijkl = C_klij,
377 : * C_ijkl = C_ijlk, C_ijkl = C_jikl
378 : * @param input is
379 : * C1111 C1122 C1133 C1123 C1113 C1112 C2222 C2233 C2223 C2213 C2212 C3333 C3323
380 : * C3313 C3312 C2323 C2313 C2312 C1313 C1312 C1212
381 : */
382 : template <typename T2>
383 : void fillSymmetric21FromInputVector(const T2 & input);
384 :
385 : /// Inner product of the major transposed tensor with a rank two tensor
386 : RankTwoTensorTempl<T> innerProductTranspose(const RankTwoTensorTempl<T> &) const;
387 :
388 : /// Sum C_ijkl M_kl for a given i,j
389 : T contractionIj(unsigned int, unsigned int, const RankTwoTensorTempl<T> &) const;
390 :
391 : /// Sum M_ij C_ijkl for a given k,l
392 : T contractionKl(unsigned int, unsigned int, const RankTwoTensorTempl<T> &) const;
393 :
394 : /// Calculates the sum of Ciijj for i and j varying from 0 to 2
395 : T sum3x3() const;
396 :
397 : /// Calculates the vector a[i] = sum over j Ciijj for i and j varying from 0 to 2
398 : libMesh::VectorValue<T> sum3x1() const;
399 :
400 : /// Calculates C_imnt A_jm B_kn C_lt
401 : RankFourTensorTempl<T> tripleProductJkl(const RankTwoTensorTempl<T> &,
402 : const RankTwoTensorTempl<T> &,
403 : const RankTwoTensorTempl<T> &) const;
404 : /// Calculates C_mjnt A_im B_kn C_lt
405 : RankFourTensorTempl<T> tripleProductIkl(const RankTwoTensorTempl<T> &,
406 : const RankTwoTensorTempl<T> &,
407 : const RankTwoTensorTempl<T> &) const;
408 : /// Calculates C_mnkt A_im B_jn C_lt
409 : RankFourTensorTempl<T> tripleProductIjl(const RankTwoTensorTempl<T> &,
410 : const RankTwoTensorTempl<T> &,
411 : const RankTwoTensorTempl<T> &) const;
412 : /// Calculates C_mntl A_im B_jn C_kt
413 : RankFourTensorTempl<T> tripleProductIjk(const RankTwoTensorTempl<T> &,
414 : const RankTwoTensorTempl<T> &,
415 : const RankTwoTensorTempl<T> &) const;
416 :
417 : /// Calculates C_mjkl A_im
418 : RankFourTensorTempl<T> singleProductI(const RankTwoTensorTempl<T> &) const;
419 : /// Calculates C_imkl A_jm
420 : RankFourTensorTempl<T> singleProductJ(const RankTwoTensorTempl<T> &) const;
421 : /// Calculates C_ijml A_km
422 : RankFourTensorTempl<T> singleProductK(const RankTwoTensorTempl<T> &) const;
423 : /// Calculates C_ijkm A_lm
424 : RankFourTensorTempl<T> singleProductL(const RankTwoTensorTempl<T> &) const;
425 :
426 : /// checks if the tensor is symmetric
427 : bool isSymmetric() const;
428 :
429 : /// checks if the tensor is isotropic
430 : bool isIsotropic() const;
431 :
432 : protected:
433 : /// The values of the rank-four tensor stored by
434 : /// index=(((i * LIBMESH_DIM + j) * LIBMESH_DIM + k) * LIBMESH_DIM + l)
435 : T _vals[N4];
436 :
437 : /**
438 : * fillAntisymmetricFromInputVector takes 6 inputs to fill the
439 : * the antisymmetric Rank-4 tensor with the appropriate symmetries maintained.
440 : * I.e., B_ijkl = -B_jikl = -B_ijlk = B_klij
441 : * @param input this is B1212, B1213, B1223, B1313, B1323, B2323.
442 : */
443 : void fillAntisymmetricFromInputVector(const std::vector<T> & input);
444 :
445 : /**
446 : * fillGeneralIsotropicFromInputVector takes 3 inputs to fill the
447 : * Rank-4 tensor with symmetries C_ijkl = C_klij, and isotropy, ie
448 : * C_ijkl = la*de_ij*de_kl + mu*(de_ik*de_jl + de_il*de_jk) + a*ep_ijm*ep_klm
449 : * where la is the first Lame modulus, mu is the second (shear) Lame modulus,
450 : * and a is the antisymmetric shear modulus, and ep is the permutation tensor
451 : * @param input this is la, mu, a in the above formula
452 : */
453 : void fillGeneralIsotropicFromInputVector(const std::vector<T> & input);
454 :
455 : /**
456 : * fillAntisymmetricIsotropicFromInputVector takes 1 input to fill the
457 : * the antisymmetric Rank-4 tensor with the appropriate symmetries maintained.
458 : * I.e., C_ijkl = a * ep_ijm * ep_klm, where epsilon is the permutation tensor (and sum on m)
459 : * @param input this is a in the above formula
460 : */
461 : void fillAntisymmetricIsotropicFromInputVector(const std::vector<T> & input);
462 :
463 : /**
464 : * fillSymmetricIsotropicFromInputVector takes 2 inputs to fill the
465 : * the symmetric Rank-4 tensor with the appropriate symmetries maintained.
466 : * C_ijkl = lambda*de_ij*de_kl + mu*(de_ik*de_jl + de_il*de_jk)
467 : * where lambda is the first Lame modulus, mu is the second (shear) Lame modulus,
468 : * @param input this is lambda and mu in the above formula
469 : */
470 : void fillSymmetricIsotropicFromInputVector(const std::vector<T> & input);
471 :
472 : /**
473 : * fillSymmetricIsotropicEandNuFromInputVector is a variation of the
474 : * fillSymmetricIsotropicFromInputVector which takes as inputs the
475 : * more commonly used Young's modulus (E) and Poisson's ratio (nu)
476 : * constants to fill the isotropic elasticity tensor. Using well-known formulas,
477 : * E and nu are used to calculate lambda and mu and then the vector is passed
478 : * to fillSymmetricIsotropicFromInputVector.
479 : * @param input Young's modulus (E) and Poisson's ratio (nu)
480 : */
481 : void fillSymmetricIsotropicEandNuFromInputVector(const std::vector<T> & input);
482 :
483 : /**
484 : * fillAxisymmetricRZFromInputVector takes 5 inputs to fill the axisymmetric
485 : * Rank-4 tensor with the appropriate symmetries maintatined for use with
486 : * axisymmetric problems using coord_type = RZ.
487 : * I.e. C1111 = C2222, C1133 = C2233, C2323 = C3131 and C1212 = 0.5*(C1111-C1122)
488 : * @param input this is C1111, C1122, C1133, C3333, C2323.
489 : */
490 : void fillAxisymmetricRZFromInputVector(const std::vector<T> & input);
491 :
492 : /**
493 : * fillGeneralFromInputVector takes 81 inputs to fill the Rank-4 tensor
494 : * No symmetries are explicitly maintained
495 : * @param input C(i,j,k,l) = input[i*N*N*N + j*N*N + k*N + l]
496 : */
497 : void fillGeneralFromInputVector(const std::vector<T> & input);
498 :
499 : /**
500 : * fillPrincipalFromInputVector takes 9 inputs to fill a Rank-4 tensor
501 : * C1111 = input0
502 : * C1122 = input1
503 : * C1133 = input2
504 : * C2211 = input3
505 : * C2222 = input4
506 : * C2233 = input5
507 : * C3311 = input6
508 : * C3322 = input7
509 : * C3333 = input8
510 : * with all other components being zero
511 : */
512 :
513 : void fillPrincipalFromInputVector(const std::vector<T> & input);
514 :
515 : /**
516 : * fillGeneralOrhotropicFromInputVector takes 10 inputs to fill the Rank-4 tensor
517 : * It defines a general orthotropic tensor for which some constraints among
518 : * elastic parameters exist
519 : * @param input Ea, Eb, Ec, Gab, Gbc, Gca, nuba, nuca, nucb, nuab, nuac, nubc
520 : */
521 : void fillGeneralOrthotropicFromInputVector(const std::vector<T> & input);
522 :
523 : template <class T2>
524 : friend void dataStore(std::ostream &, RankFourTensorTempl<T2> &, void *);
525 :
526 : template <class T2>
527 : friend void dataLoad(std::istream &, RankFourTensorTempl<T2> &, void *);
528 :
529 : template <typename T2>
530 : friend class RankTwoTensorTempl;
531 : template <typename T2>
532 : friend class RankFourTensorTempl;
533 : template <typename T2>
534 : friend class RankThreeTensorTempl;
535 : };
536 :
537 : namespace MetaPhysicL
538 : {
539 : template <typename T>
540 : struct RawType<RankFourTensorTempl<T>>
541 : {
542 : typedef RankFourTensorTempl<typename RawType<T>::value_type> value_type;
543 :
544 12 : static value_type value(const RankFourTensorTempl<T> & in)
545 : {
546 12 : constexpr auto N = RankFourTensorTempl<T>::N;
547 12 : value_type ret;
548 48 : for (auto i : libMesh::make_range(N))
549 144 : for (auto j : libMesh::make_range(N))
550 432 : for (auto k : libMesh::make_range(N))
551 1296 : for (auto l : libMesh::make_range(N))
552 972 : ret(i, j, k, l) = raw_value(in(i, j, k, l));
553 :
554 12 : return ret;
555 : }
556 : };
557 : }
558 :
559 : template <typename T1, typename T2>
560 : inline auto
561 76 : operator*(const T1 & a, const RankFourTensorTempl<T2> & b) ->
562 : typename std::enable_if<libMesh::ScalarTraits<T1>::value,
563 : RankFourTensorTempl<decltype(T1() * T2())>>::type
564 : {
565 76 : return b * a;
566 : }
567 :
568 : template <typename T>
569 : template <typename T2>
570 984 : RankFourTensorTempl<T>::RankFourTensorTempl(const RankFourTensorTempl<T2> & copy)
571 : {
572 984 : for (auto i : libMesh::make_range(N4))
573 972 : _vals[i] = copy._vals[i];
574 12 : }
575 :
576 : template <typename T>
577 : template <typename T2>
578 : auto
579 82 : RankFourTensorTempl<T>::operator*(const T2 & b) const ->
580 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
581 : RankFourTensorTempl<decltype(T() * T2())>>::type
582 : {
583 : typedef decltype(T() * T2()) ValueType;
584 82 : RankFourTensorTempl<ValueType> result;
585 :
586 6724 : for (auto i : libMesh::make_range(N4))
587 6642 : result._vals[i] = _vals[i] * b;
588 :
589 82 : return result;
590 0 : }
591 :
592 : template <typename T>
593 : template <typename T2>
594 : auto
595 10 : RankFourTensorTempl<T>::operator/(const T2 & b) const ->
596 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
597 : RankFourTensorTempl<decltype(T() / T2())>>::type
598 : {
599 10 : RankFourTensorTempl<decltype(T() / T2())> result;
600 820 : for (auto i : libMesh::make_range(N4))
601 810 : result._vals[i] = _vals[i] / b;
602 10 : return result;
603 0 : }
604 :
605 : template <typename T>
606 : template <typename T2>
607 : void
608 4 : RankFourTensorTempl<T>::fillSymmetric9FromInputVector(const T2 & input)
609 : {
610 : mooseAssert(input.size() == 9,
611 : "To use fillSymmetric9FromInputVector, your input must have size 9.");
612 4 : zero();
613 :
614 4 : (*this)(0, 0, 0, 0) = input[0]; // C1111
615 4 : (*this)(1, 1, 1, 1) = input[3]; // C2222
616 4 : (*this)(2, 2, 2, 2) = input[5]; // C3333
617 :
618 4 : (*this)(0, 0, 1, 1) = input[1]; // C1122
619 4 : (*this)(1, 1, 0, 0) = input[1];
620 :
621 4 : (*this)(0, 0, 2, 2) = input[2]; // C1133
622 4 : (*this)(2, 2, 0, 0) = input[2];
623 :
624 4 : (*this)(1, 1, 2, 2) = input[4]; // C2233
625 4 : (*this)(2, 2, 1, 1) = input[4];
626 :
627 4 : (*this)(1, 2, 1, 2) = input[6]; // C2323
628 4 : (*this)(2, 1, 2, 1) = input[6];
629 4 : (*this)(2, 1, 1, 2) = input[6];
630 4 : (*this)(1, 2, 2, 1) = input[6];
631 :
632 4 : (*this)(0, 2, 0, 2) = input[7]; // C1313
633 4 : (*this)(2, 0, 2, 0) = input[7];
634 4 : (*this)(2, 0, 0, 2) = input[7];
635 4 : (*this)(0, 2, 2, 0) = input[7];
636 :
637 4 : (*this)(0, 1, 0, 1) = input[8]; // C1212
638 4 : (*this)(1, 0, 1, 0) = input[8];
639 4 : (*this)(1, 0, 0, 1) = input[8];
640 4 : (*this)(0, 1, 1, 0) = input[8];
641 4 : }
642 : template <typename T>
643 : template <typename T2>
644 : void
645 12 : RankFourTensorTempl<T>::fillSymmetric21FromInputVector(const T2 & input)
646 : {
647 : mooseAssert(input.size() == 21,
648 : "To use fillSymmetric21FromInputVector, your input must have size 21.");
649 :
650 12 : (*this)(0, 0, 0, 0) = input[0]; // C1111
651 12 : (*this)(1, 1, 1, 1) = input[6]; // C2222
652 12 : (*this)(2, 2, 2, 2) = input[11]; // C3333
653 :
654 12 : (*this)(0, 0, 1, 1) = input[1]; // C1122
655 12 : (*this)(1, 1, 0, 0) = input[1];
656 :
657 12 : (*this)(0, 0, 2, 2) = input[2]; // C1133
658 12 : (*this)(2, 2, 0, 0) = input[2];
659 :
660 12 : (*this)(1, 1, 2, 2) = input[7]; // C2233
661 12 : (*this)(2, 2, 1, 1) = input[7];
662 :
663 12 : (*this)(0, 0, 0, 2) = input[4]; // C1113
664 12 : (*this)(0, 0, 2, 0) = input[4];
665 12 : (*this)(0, 2, 0, 0) = input[4];
666 12 : (*this)(2, 0, 0, 0) = input[4];
667 :
668 12 : (*this)(0, 0, 0, 1) = input[5]; // C1112
669 12 : (*this)(0, 0, 1, 0) = input[5];
670 12 : (*this)(0, 1, 0, 0) = input[5];
671 12 : (*this)(1, 0, 0, 0) = input[5];
672 :
673 12 : (*this)(1, 1, 1, 2) = input[8]; // C2223
674 12 : (*this)(1, 1, 2, 1) = input[8];
675 12 : (*this)(1, 2, 1, 1) = input[8];
676 12 : (*this)(2, 1, 1, 1) = input[8];
677 :
678 12 : (*this)(1, 1, 1, 0) = input[10];
679 12 : (*this)(1, 1, 0, 1) = input[10];
680 12 : (*this)(1, 0, 1, 1) = input[10];
681 12 : (*this)(0, 1, 1, 1) = input[10]; // C2212 //flipped for filling purposes
682 :
683 12 : (*this)(2, 2, 2, 1) = input[12];
684 12 : (*this)(2, 2, 1, 2) = input[12];
685 12 : (*this)(2, 1, 2, 2) = input[12];
686 12 : (*this)(1, 2, 2, 2) = input[12]; // C3323 //flipped for filling purposes
687 :
688 12 : (*this)(2, 2, 2, 0) = input[13];
689 12 : (*this)(2, 2, 0, 2) = input[13];
690 12 : (*this)(2, 0, 2, 2) = input[13];
691 12 : (*this)(0, 2, 2, 2) = input[13]; // C3313 //flipped for filling purposes
692 :
693 12 : (*this)(0, 0, 1, 2) = input[3]; // C1123
694 12 : (*this)(0, 0, 2, 1) = input[3];
695 12 : (*this)(1, 2, 0, 0) = input[3];
696 12 : (*this)(2, 1, 0, 0) = input[3];
697 :
698 12 : (*this)(1, 1, 0, 2) = input[9];
699 12 : (*this)(1, 1, 2, 0) = input[9];
700 12 : (*this)(0, 2, 1, 1) = input[9]; // C2213 //flipped for filling purposes
701 12 : (*this)(2, 0, 1, 1) = input[9];
702 :
703 12 : (*this)(2, 2, 0, 1) = input[14];
704 12 : (*this)(2, 2, 1, 0) = input[14];
705 12 : (*this)(0, 1, 2, 2) = input[14]; // C3312 //flipped for filling purposes
706 12 : (*this)(1, 0, 2, 2) = input[14];
707 :
708 12 : (*this)(1, 2, 1, 2) = input[15]; // C2323
709 12 : (*this)(2, 1, 2, 1) = input[15];
710 12 : (*this)(2, 1, 1, 2) = input[15];
711 12 : (*this)(1, 2, 2, 1) = input[15];
712 :
713 12 : (*this)(0, 2, 0, 2) = input[18]; // C1313
714 12 : (*this)(2, 0, 2, 0) = input[18];
715 12 : (*this)(2, 0, 0, 2) = input[18];
716 12 : (*this)(0, 2, 2, 0) = input[18];
717 :
718 12 : (*this)(0, 1, 0, 1) = input[20]; // C1212
719 12 : (*this)(1, 0, 1, 0) = input[20];
720 12 : (*this)(1, 0, 0, 1) = input[20];
721 12 : (*this)(0, 1, 1, 0) = input[20];
722 :
723 12 : (*this)(1, 2, 0, 2) = input[16];
724 12 : (*this)(0, 2, 1, 2) = input[16]; // C2313 //flipped for filling purposes
725 12 : (*this)(2, 1, 0, 2) = input[16];
726 12 : (*this)(1, 2, 2, 0) = input[16];
727 12 : (*this)(2, 0, 1, 2) = input[16];
728 12 : (*this)(0, 2, 2, 1) = input[16];
729 12 : (*this)(2, 1, 2, 0) = input[16];
730 12 : (*this)(2, 0, 2, 1) = input[16];
731 :
732 12 : (*this)(1, 2, 0, 1) = input[17];
733 12 : (*this)(0, 1, 1, 2) = input[17]; // C2312 //flipped for filling purposes
734 12 : (*this)(2, 1, 0, 1) = input[17];
735 12 : (*this)(1, 2, 1, 0) = input[17];
736 12 : (*this)(1, 0, 1, 2) = input[17];
737 12 : (*this)(0, 1, 2, 1) = input[17];
738 12 : (*this)(2, 1, 1, 0) = input[17];
739 12 : (*this)(1, 0, 2, 1) = input[17];
740 :
741 12 : (*this)(0, 2, 0, 1) = input[19];
742 12 : (*this)(0, 1, 0, 2) = input[19]; // C1312 //flipped for filling purposes
743 12 : (*this)(2, 0, 0, 1) = input[19];
744 12 : (*this)(0, 2, 1, 0) = input[19];
745 12 : (*this)(1, 0, 0, 2) = input[19];
746 12 : (*this)(0, 1, 2, 0) = input[19];
747 12 : (*this)(2, 0, 1, 0) = input[19];
748 12 : (*this)(1, 0, 2, 0) = input[19];
749 12 : }
750 :
751 : template <typename T>
752 : RankFourTensorTempl<T>
753 8 : RankFourTensorTempl<T>::inverse() const
754 : {
755 :
756 : // The inverse of a 3x3x3x3 in the C_ijkl*A_klmn = de_im de_jn sense is
757 : // simply the inverse of the 9x9 matrix of the tensor entries.
758 : // So all we need to do is inverse _vals (with the appropriate row-major
759 : // storage)
760 :
761 8 : RankFourTensorTempl<T> result;
762 :
763 : if constexpr (RankFourTensorTempl<T>::N4 * sizeof(T) > EIGEN_STACK_ALLOCATION_LIMIT)
764 : {
765 : // Allocate on the heap if you're going to exceed the stack size limit
766 2 : Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> mat(9, 9);
767 164 : for (auto i : libMesh::make_range(9 * 9))
768 162 : mat(i) = _vals[i];
769 :
770 2 : mat = mat.inverse();
771 :
772 164 : for (auto i : libMesh::make_range(9 * 9))
773 162 : result._vals[i] = mat(i);
774 2 : }
775 : else
776 : {
777 : // Allocate on the stack if small enough
778 6 : const Eigen::Map<const Eigen::Matrix<T, 9, 9, Eigen::RowMajor>> mat(&_vals[0]);
779 6 : Eigen::Map<Eigen::Matrix<T, 9, 9, Eigen::RowMajor>> res(&result._vals[0]);
780 6 : res = mat.inverse();
781 : }
782 :
783 8 : return result;
784 0 : }
785 :
786 : template <typename T>
787 : template <int m>
788 : RankThreeTensorTempl<T>
789 4 : RankFourTensorTempl<T>::contraction(const libMesh::VectorValue<T> & b) const
790 : {
791 4 : RankThreeTensorTempl<T> result;
792 : static constexpr std::size_t z[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
793 : std::size_t x[4];
794 16 : for (x[0] = 0; x[0] < N; ++x[0])
795 48 : for (x[1] = 0; x[1] < N; ++x[1])
796 144 : for (x[2] = 0; x[2] < N; ++x[2])
797 432 : for (x[3] = 0; x[3] < N; ++x[3])
798 324 : result(x[z[m][0]], x[z[m][1]], x[z[m][2]]) += (*this)(x[0], x[1], x[2], x[3]) * b(x[m]);
799 :
800 8 : return result;
801 : }
|
