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 1 : 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 159177514 : inline T & operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
181 : {
182 159177514 : 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 1672104 : inline const T & operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
190 : {
191 1672104 : 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 63763446 : 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 libMesh::boostcopy::enable_if_c<libMesh::ScalarTraits<Scalar>::value,
219 : RankFourTensorTempl &>::type
220 : operator=(const Scalar & libmesh_dbg_var(p))
221 : {
222 : libmesh_assert_equal_to(p, Scalar(0));
223 : this->zero();
224 : return *this;
225 : }
226 :
227 : /// C_ijkl*a_kl
228 : template <template <typename> class Tensor, typename T2>
229 : auto operator*(const Tensor<T2> & a) const ->
230 : typename std::enable_if<TwoTensorMultTraits<Tensor, T2>::value,
231 : RankTwoTensorTempl<decltype(T() * T2())>>::type;
232 :
233 : /// C_ijkl*a
234 : template <typename T2>
235 : auto operator*(const T2 & a) const ->
236 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
237 : RankFourTensorTempl<decltype(T() * T2())>>::type;
238 :
239 : /// C_ijkl *= a
240 : RankFourTensorTempl<T> & operator*=(const T & a);
241 :
242 : /// C_ijkl/a
243 : template <typename T2>
244 : auto operator/(const T2 & a) const ->
245 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
246 : RankFourTensorTempl<decltype(T() / T2())>>::type;
247 :
248 : /// C_ijkl /= a for all i, j, k, l
249 : RankFourTensorTempl<T> & operator/=(const T & a);
250 :
251 : /// C_ijkl += a_ijkl for all i, j, k, l
252 : RankFourTensorTempl<T> & operator+=(const RankFourTensorTempl<T> & a);
253 :
254 : /// C_ijkl + a_ijkl
255 : template <typename T2>
256 : auto operator+(const RankFourTensorTempl<T2> & a) const
257 : -> RankFourTensorTempl<decltype(T() + T2())>;
258 :
259 : /// C_ijkl -= a_ijkl
260 : RankFourTensorTempl<T> & operator-=(const RankFourTensorTempl<T> & a);
261 :
262 : /// C_ijkl - a_ijkl
263 : template <typename T2>
264 : auto operator-(const RankFourTensorTempl<T2> & a) const
265 : -> RankFourTensorTempl<decltype(T() - T2())>;
266 :
267 : /// -C_ijkl
268 : RankFourTensorTempl<T> operator-() const;
269 :
270 : /// C_ijpq*a_pqkl
271 : template <typename T2>
272 : auto operator*(const RankFourTensorTempl<T2> & a) const
273 : -> RankFourTensorTempl<decltype(T() * T2())>;
274 :
275 : /// sqrt(C_ijkl*C_ijkl)
276 : T L2norm() const;
277 :
278 : /**
279 : * This returns A_ijkl such that C_ijkl*A_klmn = 0.5*(de_im de_jn + de_in de_jm)
280 : * This routine assumes that C_ijkl = C_jikl = C_ijlk
281 : */
282 : RankFourTensorTempl<T> invSymm() const;
283 :
284 : /**
285 : * This returns A_ijkl such that C_ijkl*A_klmn = de_im de_jn
286 : * i.e. the general rank four inverse
287 : */
288 : RankFourTensorTempl<T> inverse() const;
289 :
290 : /**
291 : * Rotate the tensor using
292 : * C_ijkl = R_im R_jn R_ko R_lp C_mnop
293 : */
294 : void rotate(const TypeTensor<T> & R);
295 :
296 : /**
297 : * Transpose the tensor by swapping the first pair with the second pair of indices
298 : * @return C_klji
299 : */
300 : RankFourTensorTempl<T> transposeMajor() const;
301 :
302 : /**
303 : * Transpose the tensor by swapping the first two indeces
304 : * @return C_jikl
305 : */
306 : RankFourTensorTempl<T> transposeIj() const;
307 :
308 : /**
309 : * Transpose the tensor by swapping the last two indeces
310 : * @return C_ijlk
311 : */
312 : RankFourTensorTempl<T> transposeKl() const;
313 :
314 : /**
315 : * single contraction of a RankFourTensor with a vector over index m
316 : * @return C_xxx = a_ijkl*b_m where m={i,j,k,l} and xxx the remaining indices
317 : */
318 : template <int m>
319 : RankThreeTensorTempl<T> contraction(const libMesh::VectorValue<T> & b) const;
320 :
321 : /**
322 : * Fills the tensor entries ignoring the last dimension (ie, C_ijkl=0 if any of i, j, k, or l =
323 : * 3).
324 : * Fill method depends on size of input
325 : * Input size = 2. Then C_1111 = C_2222 = input[0], and C_1122 = input[1], and C_1212 = (input[0]
326 : * - input[1])/2,
327 : * and C_ijkl = C_jikl = C_ijlk = C_klij, and C_1211 = C_1222 = 0.
328 : * Input size = 9. Then C_1111 = input[0], C_1112 = input[1], C_1122 = input[3],
329 : * C_1212 = input[4], C_1222 = input[5], C_1211 = input[6]
330 : * C_2211 = input[7], C_2212 = input[8], C_2222 = input[9]
331 : * and C_ijkl = C_jikl = C_ijlk
332 : */
333 : void surfaceFillFromInputVector(const std::vector<T> & input);
334 :
335 : /// Static method for use in validParams for getting the "fill_method"
336 : static MooseEnum fillMethodEnum();
337 :
338 : /**
339 : * fillFromInputVector takes some number of inputs to fill
340 : * the Rank-4 tensor.
341 : * @param input the numbers that will be placed in the tensor
342 : * @param fill_method this can be:
343 : * antisymmetric (use fillAntisymmetricFromInputVector)
344 : * symmetric9 (use fillSymmetric9FromInputVector)
345 : * symmetric21 (use fillSymmetric21FromInputVector)
346 : * general_isotropic (use fillGeneralIsotropicFrominputVector)
347 : * symmetric_isotropic (use fillSymmetricIsotropicFromInputVector)
348 : * antisymmetric_isotropic (use fillAntisymmetricIsotropicFromInputVector)
349 : * axisymmetric_rz (use fillAxisymmetricRZFromInputVector)
350 : * general (use fillGeneralFromInputVector)
351 : * principal (use fillPrincipalFromInputVector)
352 : */
353 : void fillFromInputVector(const std::vector<T> & input, FillMethod fill_method);
354 :
355 : ///@{ Vector-less fill API functions. See docs of the corresponding ...FromInputVector methods
356 : void fillGeneralIsotropic(const T & i0, const T & i1, const T & i2);
357 : void fillAntisymmetricIsotropic(const T & i0);
358 : void fillSymmetricIsotropic(const T & i0, const T & i1);
359 : void fillSymmetricIsotropicEandNu(const T & E, const T & nu);
360 : ///@}
361 :
362 : /**
363 : * fillSymmetric9FromInputVector takes 9 inputs to fill in
364 : * the Rank-4 tensor with the appropriate crystal symmetries maintained. I.e., C_ijkl = C_klij,
365 : * C_ijkl = C_ijlk, C_ijkl = C_jikl
366 : * @param input is:
367 : * C1111 C1122 C1133 C2222 C2233 C3333 C2323 C1313 C1212
368 : * In the isotropic case this is (la is first Lame constant, mu is second (shear)
369 : * Lame constant)
370 : * la+2mu la la la+2mu la la+2mu mu mu mu
371 : */
372 : template <typename T2>
373 : void fillSymmetric9FromInputVector(const T2 & input);
374 :
375 : /**
376 : * fillSymmetric21FromInputVector takes either 21 inputs to fill in
377 : * the Rank-4 tensor with the appropriate crystal symmetries maintained. I.e., C_ijkl = C_klij,
378 : * C_ijkl = C_ijlk, C_ijkl = C_jikl
379 : * @param input is
380 : * C1111 C1122 C1133 C1123 C1113 C1112 C2222 C2233 C2223 C2213 C2212 C3333 C3323
381 : * C3313 C3312 C2323 C2313 C2312 C1313 C1312 C1212
382 : */
383 : template <typename T2>
384 : void fillSymmetric21FromInputVector(const T2 & input);
385 :
386 : /// Inner product of the major transposed tensor with a rank two tensor
387 : RankTwoTensorTempl<T> innerProductTranspose(const RankTwoTensorTempl<T> &) const;
388 :
389 : /// Sum C_ijkl M_kl for a given i,j
390 : T contractionIj(unsigned int, unsigned int, const RankTwoTensorTempl<T> &) const;
391 :
392 : /// Sum M_ij C_ijkl for a given k,l
393 : T contractionKl(unsigned int, unsigned int, const RankTwoTensorTempl<T> &) const;
394 :
395 : /// Calculates the sum of Ciijj for i and j varying from 0 to 2
396 : T sum3x3() const;
397 :
398 : /// Calculates the vector a[i] = sum over j Ciijj for i and j varying from 0 to 2
399 : libMesh::VectorValue<T> sum3x1() const;
400 :
401 : /// Calculates C_imnt A_jm B_kn C_lt
402 : RankFourTensorTempl<T> tripleProductJkl(const RankTwoTensorTempl<T> &,
403 : const RankTwoTensorTempl<T> &,
404 : const RankTwoTensorTempl<T> &) const;
405 : /// Calculates C_mjnt A_im B_kn C_lt
406 : RankFourTensorTempl<T> tripleProductIkl(const RankTwoTensorTempl<T> &,
407 : const RankTwoTensorTempl<T> &,
408 : const RankTwoTensorTempl<T> &) const;
409 : /// Calculates C_mnkt A_im B_jn C_lt
410 : RankFourTensorTempl<T> tripleProductIjl(const RankTwoTensorTempl<T> &,
411 : const RankTwoTensorTempl<T> &,
412 : const RankTwoTensorTempl<T> &) const;
413 : /// Calculates C_mntl A_im B_jn C_kt
414 : RankFourTensorTempl<T> tripleProductIjk(const RankTwoTensorTempl<T> &,
415 : const RankTwoTensorTempl<T> &,
416 : const RankTwoTensorTempl<T> &) const;
417 :
418 : /// Calculates C_mjkl A_im
419 : RankFourTensorTempl<T> singleProductI(const RankTwoTensorTempl<T> &) const;
420 : /// Calculates C_imkl A_jm
421 : RankFourTensorTempl<T> singleProductJ(const RankTwoTensorTempl<T> &) const;
422 : /// Calculates C_ijml A_km
423 : RankFourTensorTempl<T> singleProductK(const RankTwoTensorTempl<T> &) const;
424 : /// Calculates C_ijkm A_lm
425 : RankFourTensorTempl<T> singleProductL(const RankTwoTensorTempl<T> &) const;
426 :
427 : /// checks if the tensor is symmetric
428 : bool isSymmetric() const;
429 :
430 : /// checks if the tensor is isotropic
431 : bool isIsotropic() const;
432 :
433 : protected:
434 : /// The values of the rank-four tensor stored by
435 : /// index=(((i * LIBMESH_DIM + j) * LIBMESH_DIM + k) * LIBMESH_DIM + l)
436 : T _vals[N4];
437 :
438 : /**
439 : * fillAntisymmetricFromInputVector takes 6 inputs to fill the
440 : * the antisymmetric Rank-4 tensor with the appropriate symmetries maintained.
441 : * I.e., B_ijkl = -B_jikl = -B_ijlk = B_klij
442 : * @param input this is B1212, B1213, B1223, B1313, B1323, B2323.
443 : */
444 : void fillAntisymmetricFromInputVector(const std::vector<T> & input);
445 :
446 : /**
447 : * fillGeneralIsotropicFromInputVector takes 3 inputs to fill the
448 : * Rank-4 tensor with symmetries C_ijkl = C_klij, and isotropy, ie
449 : * C_ijkl = la*de_ij*de_kl + mu*(de_ik*de_jl + de_il*de_jk) + a*ep_ijm*ep_klm
450 : * where la is the first Lame modulus, mu is the second (shear) Lame modulus,
451 : * and a is the antisymmetric shear modulus, and ep is the permutation tensor
452 : * @param input this is la, mu, a in the above formula
453 : */
454 : void fillGeneralIsotropicFromInputVector(const std::vector<T> & input);
455 :
456 : /**
457 : * fillAntisymmetricIsotropicFromInputVector takes 1 input to fill the
458 : * the antisymmetric Rank-4 tensor with the appropriate symmetries maintained.
459 : * I.e., C_ijkl = a * ep_ijm * ep_klm, where epsilon is the permutation tensor (and sum on m)
460 : * @param input this is a in the above formula
461 : */
462 : void fillAntisymmetricIsotropicFromInputVector(const std::vector<T> & input);
463 :
464 : /**
465 : * fillSymmetricIsotropicFromInputVector takes 2 inputs to fill the
466 : * the symmetric Rank-4 tensor with the appropriate symmetries maintained.
467 : * C_ijkl = lambda*de_ij*de_kl + mu*(de_ik*de_jl + de_il*de_jk)
468 : * where lambda is the first Lame modulus, mu is the second (shear) Lame modulus,
469 : * @param input this is lambda and mu in the above formula
470 : */
471 : void fillSymmetricIsotropicFromInputVector(const std::vector<T> & input);
472 :
473 : /**
474 : * fillSymmetricIsotropicEandNuFromInputVector is a variation of the
475 : * fillSymmetricIsotropicFromInputVector which takes as inputs the
476 : * more commonly used Young's modulus (E) and Poisson's ratio (nu)
477 : * constants to fill the isotropic elasticity tensor. Using well-known formulas,
478 : * E and nu are used to calculate lambda and mu and then the vector is passed
479 : * to fillSymmetricIsotropicFromInputVector.
480 : * @param input Young's modulus (E) and Poisson's ratio (nu)
481 : */
482 : void fillSymmetricIsotropicEandNuFromInputVector(const std::vector<T> & input);
483 :
484 : /**
485 : * fillAxisymmetricRZFromInputVector takes 5 inputs to fill the axisymmetric
486 : * Rank-4 tensor with the appropriate symmetries maintatined for use with
487 : * axisymmetric problems using coord_type = RZ.
488 : * I.e. C1111 = C2222, C1133 = C2233, C2323 = C3131 and C1212 = 0.5*(C1111-C1122)
489 : * @param input this is C1111, C1122, C1133, C3333, C2323.
490 : */
491 : void fillAxisymmetricRZFromInputVector(const std::vector<T> & input);
492 :
493 : /**
494 : * fillGeneralFromInputVector takes 81 inputs to fill the Rank-4 tensor
495 : * No symmetries are explicitly maintained
496 : * @param input C(i,j,k,l) = input[i*N*N*N + j*N*N + k*N + l]
497 : */
498 : void fillGeneralFromInputVector(const std::vector<T> & input);
499 :
500 : /**
501 : * fillPrincipalFromInputVector takes 9 inputs to fill a Rank-4 tensor
502 : * C1111 = input0
503 : * C1122 = input1
504 : * C1133 = input2
505 : * C2211 = input3
506 : * C2222 = input4
507 : * C2233 = input5
508 : * C3311 = input6
509 : * C3322 = input7
510 : * C3333 = input8
511 : * with all other components being zero
512 : */
513 :
514 : void fillPrincipalFromInputVector(const std::vector<T> & input);
515 :
516 : /**
517 : * fillGeneralOrhotropicFromInputVector takes 10 inputs to fill the Rank-4 tensor
518 : * It defines a general orthotropic tensor for which some constraints among
519 : * elastic parameters exist
520 : * @param input Ea, Eb, Ec, Gab, Gbc, Gca, nuba, nuca, nucb, nuab, nuac, nubc
521 : */
522 : void fillGeneralOrthotropicFromInputVector(const std::vector<T> & input);
523 :
524 : template <class T2>
525 : friend void dataStore(std::ostream &, RankFourTensorTempl<T2> &, void *);
526 :
527 : template <class T2>
528 : friend void dataLoad(std::istream &, RankFourTensorTempl<T2> &, void *);
529 :
530 : template <typename T2>
531 : friend class RankTwoTensorTempl;
532 : template <typename T2>
533 : friend class RankFourTensorTempl;
534 : template <typename T2>
535 : friend class RankThreeTensorTempl;
536 : };
537 :
538 : namespace MetaPhysicL
539 : {
540 : template <typename T>
541 : struct RawType<RankFourTensorTempl<T>>
542 : {
543 : typedef RankFourTensorTempl<typename RawType<T>::value_type> value_type;
544 :
545 6 : static value_type value(const RankFourTensorTempl<T> & in)
546 : {
547 6 : constexpr auto N = RankFourTensorTempl<T>::N;
548 6 : value_type ret;
549 24 : for (auto i : libMesh::make_range(N))
550 72 : for (auto j : libMesh::make_range(N))
551 216 : for (auto k : libMesh::make_range(N))
552 648 : for (auto l : libMesh::make_range(N))
553 486 : ret(i, j, k, l) = raw_value(in(i, j, k, l));
554 :
555 6 : return ret;
556 : }
557 : };
558 : }
559 :
560 : template <typename T1, typename T2>
561 : inline auto
562 38 : operator*(const T1 & a, const RankFourTensorTempl<T2> & b) ->
563 : typename std::enable_if<libMesh::ScalarTraits<T1>::value,
564 : RankFourTensorTempl<decltype(T1() * T2())>>::type
565 : {
566 38 : return b * a;
567 : }
568 :
569 : template <typename T>
570 : template <typename T2>
571 492 : RankFourTensorTempl<T>::RankFourTensorTempl(const RankFourTensorTempl<T2> & copy)
572 : {
573 492 : for (auto i : libMesh::make_range(N4))
574 486 : _vals[i] = copy._vals[i];
575 6 : }
576 :
577 : template <typename T>
578 : template <typename T2>
579 : auto
580 41 : RankFourTensorTempl<T>::operator*(const T2 & b) const ->
581 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
582 : RankFourTensorTempl<decltype(T() * T2())>>::type
583 : {
584 : typedef decltype(T() * T2()) ValueType;
585 41 : RankFourTensorTempl<ValueType> result;
586 :
587 3362 : for (auto i : libMesh::make_range(N4))
588 3321 : result._vals[i] = _vals[i] * b;
589 :
590 41 : return result;
591 0 : }
592 :
593 : template <typename T>
594 : template <typename T2>
595 : auto
596 5 : RankFourTensorTempl<T>::operator/(const T2 & b) const ->
597 : typename std::enable_if<libMesh::ScalarTraits<T2>::value,
598 : RankFourTensorTempl<decltype(T() / T2())>>::type
599 : {
600 5 : RankFourTensorTempl<decltype(T() / T2())> result;
601 410 : for (auto i : libMesh::make_range(N4))
602 405 : result._vals[i] = _vals[i] / b;
603 5 : return result;
604 0 : }
605 :
606 : template <typename T>
607 : template <typename T2>
608 : void
609 2 : RankFourTensorTempl<T>::fillSymmetric9FromInputVector(const T2 & input)
610 : {
611 : mooseAssert(input.size() == 9,
612 : "To use fillSymmetric9FromInputVector, your input must have size 9.");
613 2 : zero();
614 :
615 2 : (*this)(0, 0, 0, 0) = input[0]; // C1111
616 2 : (*this)(1, 1, 1, 1) = input[3]; // C2222
617 2 : (*this)(2, 2, 2, 2) = input[5]; // C3333
618 :
619 2 : (*this)(0, 0, 1, 1) = input[1]; // C1122
620 2 : (*this)(1, 1, 0, 0) = input[1];
621 :
622 2 : (*this)(0, 0, 2, 2) = input[2]; // C1133
623 2 : (*this)(2, 2, 0, 0) = input[2];
624 :
625 2 : (*this)(1, 1, 2, 2) = input[4]; // C2233
626 2 : (*this)(2, 2, 1, 1) = input[4];
627 :
628 2 : (*this)(1, 2, 1, 2) = input[6]; // C2323
629 2 : (*this)(2, 1, 2, 1) = input[6];
630 2 : (*this)(2, 1, 1, 2) = input[6];
631 2 : (*this)(1, 2, 2, 1) = input[6];
632 :
633 2 : (*this)(0, 2, 0, 2) = input[7]; // C1313
634 2 : (*this)(2, 0, 2, 0) = input[7];
635 2 : (*this)(2, 0, 0, 2) = input[7];
636 2 : (*this)(0, 2, 2, 0) = input[7];
637 :
638 2 : (*this)(0, 1, 0, 1) = input[8]; // C1212
639 2 : (*this)(1, 0, 1, 0) = input[8];
640 2 : (*this)(1, 0, 0, 1) = input[8];
641 2 : (*this)(0, 1, 1, 0) = input[8];
642 2 : }
643 : template <typename T>
644 : template <typename T2>
645 : void
646 6 : RankFourTensorTempl<T>::fillSymmetric21FromInputVector(const T2 & input)
647 : {
648 : mooseAssert(input.size() == 21,
649 : "To use fillSymmetric21FromInputVector, your input must have size 21.");
650 :
651 6 : (*this)(0, 0, 0, 0) = input[0]; // C1111
652 6 : (*this)(1, 1, 1, 1) = input[6]; // C2222
653 6 : (*this)(2, 2, 2, 2) = input[11]; // C3333
654 :
655 6 : (*this)(0, 0, 1, 1) = input[1]; // C1122
656 6 : (*this)(1, 1, 0, 0) = input[1];
657 :
658 6 : (*this)(0, 0, 2, 2) = input[2]; // C1133
659 6 : (*this)(2, 2, 0, 0) = input[2];
660 :
661 6 : (*this)(1, 1, 2, 2) = input[7]; // C2233
662 6 : (*this)(2, 2, 1, 1) = input[7];
663 :
664 6 : (*this)(0, 0, 0, 2) = input[4]; // C1113
665 6 : (*this)(0, 0, 2, 0) = input[4];
666 6 : (*this)(0, 2, 0, 0) = input[4];
667 6 : (*this)(2, 0, 0, 0) = input[4];
668 :
669 6 : (*this)(0, 0, 0, 1) = input[5]; // C1112
670 6 : (*this)(0, 0, 1, 0) = input[5];
671 6 : (*this)(0, 1, 0, 0) = input[5];
672 6 : (*this)(1, 0, 0, 0) = input[5];
673 :
674 6 : (*this)(1, 1, 1, 2) = input[8]; // C2223
675 6 : (*this)(1, 1, 2, 1) = input[8];
676 6 : (*this)(1, 2, 1, 1) = input[8];
677 6 : (*this)(2, 1, 1, 1) = input[8];
678 :
679 6 : (*this)(1, 1, 1, 0) = input[10];
680 6 : (*this)(1, 1, 0, 1) = input[10];
681 6 : (*this)(1, 0, 1, 1) = input[10];
682 6 : (*this)(0, 1, 1, 1) = input[10]; // C2212 //flipped for filling purposes
683 :
684 6 : (*this)(2, 2, 2, 1) = input[12];
685 6 : (*this)(2, 2, 1, 2) = input[12];
686 6 : (*this)(2, 1, 2, 2) = input[12];
687 6 : (*this)(1, 2, 2, 2) = input[12]; // C3323 //flipped for filling purposes
688 :
689 6 : (*this)(2, 2, 2, 0) = input[13];
690 6 : (*this)(2, 2, 0, 2) = input[13];
691 6 : (*this)(2, 0, 2, 2) = input[13];
692 6 : (*this)(0, 2, 2, 2) = input[13]; // C3313 //flipped for filling purposes
693 :
694 6 : (*this)(0, 0, 1, 2) = input[3]; // C1123
695 6 : (*this)(0, 0, 2, 1) = input[3];
696 6 : (*this)(1, 2, 0, 0) = input[3];
697 6 : (*this)(2, 1, 0, 0) = input[3];
698 :
699 6 : (*this)(1, 1, 0, 2) = input[9];
700 6 : (*this)(1, 1, 2, 0) = input[9];
701 6 : (*this)(0, 2, 1, 1) = input[9]; // C2213 //flipped for filling purposes
702 6 : (*this)(2, 0, 1, 1) = input[9];
703 :
704 6 : (*this)(2, 2, 0, 1) = input[14];
705 6 : (*this)(2, 2, 1, 0) = input[14];
706 6 : (*this)(0, 1, 2, 2) = input[14]; // C3312 //flipped for filling purposes
707 6 : (*this)(1, 0, 2, 2) = input[14];
708 :
709 6 : (*this)(1, 2, 1, 2) = input[15]; // C2323
710 6 : (*this)(2, 1, 2, 1) = input[15];
711 6 : (*this)(2, 1, 1, 2) = input[15];
712 6 : (*this)(1, 2, 2, 1) = input[15];
713 :
714 6 : (*this)(0, 2, 0, 2) = input[18]; // C1313
715 6 : (*this)(2, 0, 2, 0) = input[18];
716 6 : (*this)(2, 0, 0, 2) = input[18];
717 6 : (*this)(0, 2, 2, 0) = input[18];
718 :
719 6 : (*this)(0, 1, 0, 1) = input[20]; // C1212
720 6 : (*this)(1, 0, 1, 0) = input[20];
721 6 : (*this)(1, 0, 0, 1) = input[20];
722 6 : (*this)(0, 1, 1, 0) = input[20];
723 :
724 6 : (*this)(1, 2, 0, 2) = input[16];
725 6 : (*this)(0, 2, 1, 2) = input[16]; // C2313 //flipped for filling purposes
726 6 : (*this)(2, 1, 0, 2) = input[16];
727 6 : (*this)(1, 2, 2, 0) = input[16];
728 6 : (*this)(2, 0, 1, 2) = input[16];
729 6 : (*this)(0, 2, 2, 1) = input[16];
730 6 : (*this)(2, 1, 2, 0) = input[16];
731 6 : (*this)(2, 0, 2, 1) = input[16];
732 :
733 6 : (*this)(1, 2, 0, 1) = input[17];
734 6 : (*this)(0, 1, 1, 2) = input[17]; // C2312 //flipped for filling purposes
735 6 : (*this)(2, 1, 0, 1) = input[17];
736 6 : (*this)(1, 2, 1, 0) = input[17];
737 6 : (*this)(1, 0, 1, 2) = input[17];
738 6 : (*this)(0, 1, 2, 1) = input[17];
739 6 : (*this)(2, 1, 1, 0) = input[17];
740 6 : (*this)(1, 0, 2, 1) = input[17];
741 :
742 6 : (*this)(0, 2, 0, 1) = input[19];
743 6 : (*this)(0, 1, 0, 2) = input[19]; // C1312 //flipped for filling purposes
744 6 : (*this)(2, 0, 0, 1) = input[19];
745 6 : (*this)(0, 2, 1, 0) = input[19];
746 6 : (*this)(1, 0, 0, 2) = input[19];
747 6 : (*this)(0, 1, 2, 0) = input[19];
748 6 : (*this)(2, 0, 1, 0) = input[19];
749 6 : (*this)(1, 0, 2, 0) = input[19];
750 6 : }
751 :
752 : template <typename T>
753 : RankFourTensorTempl<T>
754 4 : RankFourTensorTempl<T>::inverse() const
755 : {
756 :
757 : // The inverse of a 3x3x3x3 in the C_ijkl*A_klmn = de_im de_jn sense is
758 : // simply the inverse of the 9x9 matrix of the tensor entries.
759 : // So all we need to do is inverse _vals (with the appropriate row-major
760 : // storage)
761 :
762 4 : RankFourTensorTempl<T> result;
763 :
764 : if constexpr (RankFourTensorTempl<T>::N4 * sizeof(T) > EIGEN_STACK_ALLOCATION_LIMIT)
765 : {
766 : // Allocate on the heap if you're going to exceed the stack size limit
767 1 : Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> mat(9, 9);
768 82 : for (auto i : libMesh::make_range(9 * 9))
769 81 : mat(i) = _vals[i];
770 :
771 1 : mat = mat.inverse();
772 :
773 82 : for (auto i : libMesh::make_range(9 * 9))
774 81 : result._vals[i] = mat(i);
775 1 : }
776 : else
777 : {
778 : // Allocate on the stack if small enough
779 3 : const Eigen::Map<const Eigen::Matrix<T, 9, 9, Eigen::RowMajor>> mat(&_vals[0]);
780 3 : Eigen::Map<Eigen::Matrix<T, 9, 9, Eigen::RowMajor>> res(&result._vals[0]);
781 3 : res = mat.inverse();
782 : }
783 :
784 4 : return result;
785 0 : }
786 :
787 : template <typename T>
788 : template <int m>
789 : RankThreeTensorTempl<T>
790 2 : RankFourTensorTempl<T>::contraction(const libMesh::VectorValue<T> & b) const
791 : {
792 2 : RankThreeTensorTempl<T> result;
793 : static constexpr std::size_t z[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
794 : std::size_t x[4];
795 8 : for (x[0] = 0; x[0] < N; ++x[0])
796 24 : for (x[1] = 0; x[1] < N; ++x[1])
797 72 : for (x[2] = 0; x[2] < N; ++x[2])
798 216 : for (x[3] = 0; x[3] < N; ++x[3])
799 162 : 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]);
800 :
801 4 : return result;
802 : }
|
