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 "RankFourTensor.h"
13 :
14 : // MOOSE includes
15 : #include "RankTwoTensor.h"
16 : #include "RankThreeTensor.h"
17 : #include "MooseEnum.h"
18 : #include "MooseException.h"
19 : #include "MooseUtils.h"
20 : #include "MatrixTools.h"
21 : #include "PermutationTensor.h"
22 :
23 : #include "libmesh/utility.h"
24 : #include "libmesh/tensor_value.h"
25 : #include "libmesh/vector_value.h"
26 :
27 : // Eigen needs LU
28 : #include "Eigen/LU"
29 :
30 : // C++ includes
31 : #include <iomanip>
32 : #include <ostream>
33 :
34 : namespace MathUtils
35 : {
36 : template <>
37 : void mooseSetToZero<RankFourTensorTempl<Real>>(RankFourTensorTempl<Real> & v);
38 : template <>
39 : void mooseSetToZero<RankFourTensorTempl<ADReal>>(RankFourTensorTempl<ADReal> & v);
40 : }
41 :
42 : template <typename T>
43 : MooseEnum
44 0 : RankFourTensorTempl<T>::fillMethodEnum()
45 : {
46 : return MooseEnum("antisymmetric symmetric9 symmetric21 general_isotropic symmetric_isotropic "
47 : "symmetric_isotropic_E_nu antisymmetric_isotropic axisymmetric_rz general "
48 0 : "principal orthotropic");
49 : }
50 :
51 : template <typename T>
52 2482004 : RankFourTensorTempl<T>::RankFourTensorTempl()
53 : {
54 : mooseAssert(N == 3, "RankFourTensorTempl<T> is currently only tested for 3 dimensions.");
55 :
56 196151708 : for (auto i : make_range(N4))
57 193759614 : _vals[i] = 0.0;
58 2392094 : }
59 :
60 : template <typename T>
61 31 : RankFourTensorTempl<T>::RankFourTensorTempl(const InitMethod init)
62 : {
63 31 : unsigned int index = 0;
64 31 : switch (init)
65 : {
66 0 : case initNone:
67 0 : break;
68 :
69 1 : case initIdentity:
70 1 : zero();
71 4 : for (auto i : make_range(N))
72 3 : (*this)(i, i, i, i) = 1.0;
73 1 : break;
74 :
75 27 : case initIdentityFour:
76 108 : for (auto i : make_range(N))
77 324 : for (auto j : make_range(N))
78 972 : for (auto k : make_range(N))
79 2916 : for (auto l : make_range(N))
80 2187 : _vals[index++] = Real(i == k && j == l);
81 27 : break;
82 :
83 3 : case initIdentitySymmetricFour:
84 12 : for (auto i : make_range(N))
85 36 : for (auto j : make_range(N))
86 108 : for (auto k : make_range(N))
87 324 : for (auto l : make_range(N))
88 243 : _vals[index++] = 0.5 * Real(i == k && j == l) + 0.5 * Real(i == l && j == k);
89 3 : break;
90 :
91 0 : case initIdentityDeviatoric:
92 0 : for (unsigned int i = 0; i < N; ++i)
93 0 : for (unsigned int j = 0; j < N; ++j)
94 0 : for (unsigned int k = 0; k < N; ++k)
95 0 : for (unsigned int l = 0; l < N; ++l)
96 : {
97 0 : _vals[index] = Real(i == k && j == l);
98 0 : if ((i == j) && (k == l))
99 0 : _vals[index] -= 1.0 / 3.0;
100 0 : index++;
101 : }
102 0 : break;
103 :
104 0 : default:
105 0 : mooseError("Unknown RankFourTensorTempl<T> initialization pattern.");
106 : }
107 31 : }
108 :
109 : template <typename T>
110 17 : RankFourTensorTempl<T>::RankFourTensorTempl(const std::vector<T> & input, FillMethod fill_method)
111 : {
112 17 : fillFromInputVector(input, fill_method);
113 17 : }
114 :
115 : template <typename T>
116 : void
117 4 : RankFourTensorTempl<T>::zero()
118 : {
119 328 : for (auto i : make_range(N4))
120 324 : _vals[i] = 0.0;
121 4 : }
122 :
123 : template <typename T>
124 : template <template <typename> class Tensor, typename T2>
125 : auto
126 2 : RankFourTensorTempl<T>::operator*(const Tensor<T2> & b) const ->
127 : typename std::enable_if<TwoTensorMultTraits<Tensor, T2>::value,
128 : RankTwoTensorTempl<decltype(T() * T2())>>::type
129 : {
130 : typedef decltype(T() * T2()) ValueType;
131 2 : RankTwoTensorTempl<ValueType> result;
132 :
133 2 : unsigned int index = 0;
134 20 : for (unsigned int ij = 0; ij < N2; ++ij)
135 : {
136 18 : ValueType tmp = 0;
137 180 : for (unsigned int kl = 0; kl < N2; ++kl)
138 162 : tmp += _vals[index++] * b(kl / LIBMESH_DIM, kl % LIBMESH_DIM);
139 18 : result._coords[ij] = tmp;
140 : }
141 :
142 2 : return result;
143 0 : }
144 :
145 : template <typename T>
146 : RankFourTensorTempl<T> &
147 3 : RankFourTensorTempl<T>::operator*=(const T & a)
148 : {
149 246 : for (auto i : make_range(N4))
150 243 : _vals[i] *= a;
151 3 : return *this;
152 : }
153 :
154 : template <typename T>
155 : RankFourTensorTempl<T> &
156 1 : RankFourTensorTempl<T>::operator/=(const T & a)
157 : {
158 82 : for (auto i : make_range(N4))
159 81 : _vals[i] /= a;
160 1 : return *this;
161 : }
162 :
163 : template <typename T>
164 : RankFourTensorTempl<T> &
165 37 : RankFourTensorTempl<T>::operator+=(const RankFourTensorTempl<T> & a)
166 : {
167 3034 : for (auto i : make_range(N4))
168 2997 : _vals[i] += a._vals[i];
169 37 : return *this;
170 : }
171 :
172 : template <typename T>
173 : template <typename T2>
174 : auto
175 55 : RankFourTensorTempl<T>::operator+(const RankFourTensorTempl<T2> & b) const
176 : -> RankFourTensorTempl<decltype(T() + T2())>
177 : {
178 55 : RankFourTensorTempl<decltype(T() + T2())> result;
179 4510 : for (auto i : make_range(N4))
180 4455 : result._vals[i] = _vals[i] + b._vals[i];
181 55 : return result;
182 0 : }
183 :
184 : template <typename T>
185 : RankFourTensorTempl<T> &
186 1 : RankFourTensorTempl<T>::operator-=(const RankFourTensorTempl<T> & a)
187 : {
188 82 : for (auto i : make_range(N4))
189 81 : _vals[i] -= a._vals[i];
190 1 : return *this;
191 : }
192 :
193 : template <typename T>
194 : template <typename T2>
195 : auto
196 792105 : RankFourTensorTempl<T>::operator-(const RankFourTensorTempl<T2> & b) const
197 : -> RankFourTensorTempl<decltype(T() - T2())>
198 : {
199 792105 : RankFourTensorTempl<decltype(T() - T2())> result;
200 64952610 : for (auto i : make_range(N4))
201 64160505 : result._vals[i] = _vals[i] - b._vals[i];
202 792105 : return result;
203 0 : }
204 :
205 : template <typename T>
206 : RankFourTensorTempl<T>
207 1 : RankFourTensorTempl<T>::operator-() const
208 : {
209 1 : RankFourTensorTempl<T> result;
210 82 : for (auto i : make_range(N4))
211 81 : result._vals[i] = -_vals[i];
212 1 : return result;
213 0 : }
214 :
215 : template <typename T>
216 : template <typename T2>
217 : auto
218 4 : RankFourTensorTempl<T>::operator*(const RankFourTensorTempl<T2> & b) const
219 : -> RankFourTensorTempl<decltype(T() * T2())>
220 : {
221 : typedef decltype(T() * T2()) ValueType;
222 4 : RankFourTensorTempl<ValueType> result;
223 :
224 16 : for (auto i : make_range(N))
225 48 : for (auto j : make_range(N))
226 144 : for (auto k : make_range(N))
227 432 : for (auto l : make_range(N))
228 1296 : for (auto p : make_range(N))
229 3888 : for (auto q : make_range(N))
230 2916 : result(i, j, k, l) += (*this)(i, j, p, q) * b(p, q, k, l);
231 :
232 4 : return result;
233 0 : }
234 :
235 : template <typename T>
236 : T
237 792105 : RankFourTensorTempl<T>::L2norm() const
238 : {
239 792105 : T l2 = 0;
240 :
241 64952610 : for (auto i : make_range(N4))
242 64160505 : l2 += Utility::pow<2>(_vals[i]);
243 :
244 792105 : return std::sqrt(l2);
245 0 : }
246 :
247 : template <typename T>
248 : RankFourTensorTempl<T>
249 0 : RankFourTensorTempl<T>::invSymm() const
250 : {
251 0 : mooseError("The invSymm operation calls to LAPACK and only supports plain Real type tensors.");
252 : }
253 :
254 : template <>
255 : RankFourTensorTempl<Real>
256 3 : RankFourTensorTempl<Real>::invSymm() const
257 : {
258 3 : unsigned int ntens = N * (N + 1) / 2;
259 3 : int nskip = N - 1;
260 :
261 3 : RankFourTensorTempl<Real> result;
262 3 : std::vector<PetscScalar> mat;
263 3 : mat.assign(ntens * ntens, 0);
264 :
265 : // We use the LAPACK matrix inversion routine here. Form the matrix
266 : //
267 : // mat[0] mat[1] mat[2] mat[3] mat[4] mat[5]
268 : // mat[6] mat[7] mat[8] mat[9] mat[10] mat[11]
269 : // mat[12] mat[13] mat[14] mat[15] mat[16] mat[17]
270 : // mat[18] mat[19] mat[20] mat[21] mat[22] mat[23]
271 : // mat[24] mat[25] mat[26] mat[27] mat[28] mat[29]
272 : // mat[30] mat[31] mat[32] mat[33] mat[34] mat[35]
273 : //
274 : // This is filled from the indpendent components of C assuming
275 : // the symmetry C_ijkl = C_ijlk = C_jikl.
276 : //
277 : // If there are two rank-four tensors X and Y then the reason for
278 : // this filling becomes apparent if we want to calculate
279 : // X_ijkl*Y_klmn = Z_ijmn
280 : // For denote the "mat" versions of X, Y and Z by x, y and z.
281 : // Then
282 : // z_ab = x_ac*y_cb
283 : // Eg
284 : // z_00 = Z_0000 = X_0000*Y_0000 + X_0011*Y_1111 + X_0022*Y_2200 + 2*X_0001*Y_0100 +
285 : // 2*X_0002*Y_0200 + 2*X_0012*Y_1200 (the factors of 2 come from the assumed symmetries)
286 : // z_03 = 2*Z_0001 = X_0000*2*Y_0001 + X_0011*2*Y_1101 + X_0022*2*Y_2201 + 2*X_0001*2*Y_0101 +
287 : // 2*X_0002*2*Y_0201 + 2*X_0012*2*Y_1201
288 : // z_22 = 2*Z_0102 = X_0100*2*Y_0002 + X_0111*2*X_1102 + X_0122*2*Y_2202 + 2*X_0101*2*Y_0102 +
289 : // 2*X_0102*2*Y_0202 + 2*X_0112*2*Y_1202
290 : // Finally, we use LAPACK to find x^-1, and put it back into rank-4 tensor form
291 : //
292 : // mat[0] = C(0,0,0,0)
293 : // mat[1] = C(0,0,1,1)
294 : // mat[2] = C(0,0,2,2)
295 : // mat[3] = C(0,0,0,1)*2
296 : // mat[4] = C(0,0,0,2)*2
297 : // mat[5] = C(0,0,1,2)*2
298 :
299 : // mat[6] = C(1,1,0,0)
300 : // mat[7] = C(1,1,1,1)
301 : // mat[8] = C(1,1,2,2)
302 : // mat[9] = C(1,1,0,1)*2
303 : // mat[10] = C(1,1,0,2)*2
304 : // mat[11] = C(1,1,1,2)*2
305 :
306 : // mat[12] = C(2,2,0,0)
307 : // mat[13] = C(2,2,1,1)
308 : // mat[14] = C(2,2,2,2)
309 : // mat[15] = C(2,2,0,1)*2
310 : // mat[16] = C(2,2,0,2)*2
311 : // mat[17] = C(2,2,1,2)*2
312 :
313 : // mat[18] = C(0,1,0,0)
314 : // mat[19] = C(0,1,1,1)
315 : // mat[20] = C(0,1,2,2)
316 : // mat[21] = C(0,1,0,1)*2
317 : // mat[22] = C(0,1,0,2)*2
318 : // mat[23] = C(0,1,1,2)*2
319 :
320 : // mat[24] = C(0,2,0,0)
321 : // mat[25] = C(0,2,1,1)
322 : // mat[26] = C(0,2,2,2)
323 : // mat[27] = C(0,2,0,1)*2
324 : // mat[28] = C(0,2,0,2)*2
325 : // mat[29] = C(0,2,1,2)*2
326 :
327 : // mat[30] = C(1,2,0,0)
328 : // mat[31] = C(1,2,1,1)
329 : // mat[32] = C(1,2,2,2)
330 : // mat[33] = C(1,2,0,1)*2
331 : // mat[34] = C(1,2,0,2)*2
332 : // mat[35] = C(1,2,1,2)*2
333 :
334 3 : unsigned int index = 0;
335 12 : for (auto i : make_range(N))
336 36 : for (auto j : make_range(N))
337 108 : for (auto k : make_range(N))
338 324 : for (auto l : make_range(N))
339 : {
340 243 : if (i == j)
341 81 : mat[k == l ? i * ntens + k : i * ntens + k + nskip + l] += _vals[index];
342 : else
343 : // i!=j
344 162 : mat[k == l ? (nskip + i + j) * ntens + k : (nskip + i + j) * ntens + k + nskip + l] +=
345 162 : _vals[index]; // note the +=, which results in double-counting and is rectified
346 : // below
347 243 : index++;
348 : }
349 :
350 12 : for (unsigned int i = 3; i < ntens; ++i)
351 63 : for (auto j : make_range(ntens))
352 54 : mat[i * ntens + j] /= 2.0; // because of double-counting above
353 :
354 : // use LAPACK to find the inverse
355 3 : MatrixTools::inverse(mat, ntens);
356 :
357 : // build the resulting rank-four tensor
358 : // using the inverse of the above algorithm
359 3 : index = 0;
360 12 : for (auto i : make_range(N))
361 36 : for (auto j : make_range(N))
362 108 : for (auto k : make_range(N))
363 324 : for (auto l : make_range(N))
364 : {
365 243 : if (i == j)
366 81 : result._vals[index] =
367 81 : k == l ? mat[i * ntens + k] : mat[i * ntens + k + nskip + l] / 2.0;
368 : else
369 : // i!=j
370 270 : result._vals[index] = k == l ? mat[(nskip + i + j) * ntens + k]
371 108 : : mat[(nskip + i + j) * ntens + k + nskip + l] / 2.0;
372 243 : index++;
373 : }
374 :
375 6 : return result;
376 3 : }
377 :
378 : template <typename T>
379 : void
380 1 : RankFourTensorTempl<T>::rotate(const TypeTensor<T> & R)
381 : {
382 1 : RankFourTensorTempl<T> old = *this;
383 :
384 1 : unsigned int index = 0;
385 4 : for (auto i : make_range(N))
386 12 : for (auto j : make_range(N))
387 36 : for (auto k : make_range(N))
388 108 : for (auto l : make_range(N))
389 : {
390 81 : unsigned int index2 = 0;
391 81 : T sum = 0.0;
392 324 : for (auto m : make_range(N))
393 : {
394 243 : const T & a = R(i, m);
395 972 : for (auto n : make_range(N))
396 : {
397 729 : const T & ab = a * R(j, n);
398 2916 : for (auto o : make_range(N))
399 : {
400 2187 : const T & abc = ab * R(k, o);
401 8748 : for (auto p : make_range(N))
402 6561 : sum += abc * R(l, p) * old._vals[index2++];
403 : }
404 : }
405 : }
406 81 : _vals[index++] = sum;
407 : }
408 1 : }
409 :
410 : template <typename T>
411 : void
412 0 : RankFourTensorTempl<T>::print(std::ostream & stm) const
413 : {
414 0 : for (auto i : make_range(N))
415 0 : for (auto j : make_range(N))
416 : {
417 0 : stm << "i = " << i << " j = " << j << '\n';
418 0 : for (auto k : make_range(N))
419 : {
420 0 : for (auto l : make_range(N))
421 0 : stm << std::setw(15) << (*this)(i, j, k, l) << " ";
422 :
423 0 : stm << '\n';
424 : }
425 : }
426 :
427 0 : stm << std::flush;
428 0 : }
429 :
430 : template <typename T>
431 : void
432 1 : RankFourTensorTempl<T>::printReal(std::ostream & stm) const
433 : {
434 4 : for (unsigned int i = 0; i < N; ++i)
435 12 : for (unsigned int j = 0; j < N; ++j)
436 : {
437 9 : stm << "i = " << i << " j = " << j << '\n';
438 36 : for (unsigned int k = 0; k < N; ++k)
439 : {
440 108 : for (unsigned int l = 0; l < N; ++l)
441 81 : stm << std::setw(15) << MetaPhysicL::raw_value((*this)(i, j, k, l)) << " ";
442 :
443 27 : stm << '\n';
444 : }
445 : }
446 :
447 1 : stm << std::flush;
448 1 : }
449 :
450 : template <typename T>
451 : RankFourTensorTempl<T>
452 1 : RankFourTensorTempl<T>::transposeMajor() const
453 : {
454 1 : RankFourTensorTempl<T> result;
455 :
456 1 : unsigned int index = 0;
457 4 : for (auto i : make_range(N))
458 12 : for (auto j : make_range(N))
459 36 : for (auto k : make_range(N))
460 108 : for (auto l : make_range(N))
461 81 : result._vals[index++] = _vals[k * N3 + i * N + j + l * N2];
462 :
463 1 : return result;
464 0 : }
465 :
466 : template <typename T>
467 : RankFourTensorTempl<T>
468 1 : RankFourTensorTempl<T>::transposeIj() const
469 : {
470 1 : RankFourTensorTempl<T> result;
471 :
472 4 : for (auto i : make_range(N))
473 12 : for (auto j : make_range(N))
474 36 : for (auto k : make_range(N))
475 108 : for (auto l : make_range(N))
476 81 : result(i, j, k, l) = (*this)(j, i, k, l);
477 :
478 1 : return result;
479 0 : }
480 :
481 : template <typename T>
482 : RankFourTensorTempl<T>
483 0 : RankFourTensorTempl<T>::transposeKl() const
484 : {
485 0 : RankFourTensorTempl<T> result;
486 :
487 0 : for (auto i : make_range(N))
488 0 : for (auto j : make_range(N))
489 0 : for (auto k : make_range(N))
490 0 : for (auto l : make_range(N))
491 0 : result(i, j, k, l) = (*this)(i, j, l, k);
492 :
493 0 : return result;
494 0 : }
495 :
496 : template <typename T>
497 : void
498 0 : RankFourTensorTempl<T>::surfaceFillFromInputVector(const std::vector<T> & input)
499 : {
500 0 : zero();
501 :
502 0 : if (input.size() == 9)
503 : {
504 : // then fill from vector C_1111, C_1112, C_1122, C_1212, C_1222, C_1211, C_2211, C_2212, C_2222
505 0 : (*this)(0, 0, 0, 0) = input[0];
506 0 : (*this)(0, 0, 0, 1) = input[1];
507 0 : (*this)(0, 0, 1, 1) = input[2];
508 0 : (*this)(0, 1, 0, 1) = input[3];
509 0 : (*this)(0, 1, 1, 1) = input[4];
510 0 : (*this)(0, 1, 0, 0) = input[5];
511 0 : (*this)(1, 1, 0, 0) = input[6];
512 0 : (*this)(1, 1, 0, 1) = input[7];
513 0 : (*this)(1, 1, 1, 1) = input[8];
514 :
515 : // fill in remainders from C_ijkl = C_ijlk = C_jikl
516 0 : (*this)(0, 0, 1, 0) = (*this)(0, 0, 0, 1);
517 0 : (*this)(0, 1, 1, 0) = (*this)(0, 1, 0, 1);
518 0 : (*this)(1, 0, 0, 0) = (*this)(0, 1, 0, 0);
519 0 : (*this)(1, 0, 0, 1) = (*this)(0, 1, 0, 1);
520 0 : (*this)(1, 0, 1, 1) = (*this)(0, 1, 1, 1);
521 0 : (*this)(1, 0, 0, 0) = (*this)(0, 1, 0, 0);
522 0 : (*this)(1, 1, 1, 0) = (*this)(1, 1, 0, 1);
523 : }
524 0 : else if (input.size() == 2)
525 : {
526 : // only two independent constants, C_1111 and C_1122
527 0 : (*this)(0, 0, 0, 0) = input[0];
528 0 : (*this)(0, 0, 1, 1) = input[1];
529 : // use symmetries
530 0 : (*this)(1, 1, 1, 1) = (*this)(0, 0, 0, 0);
531 0 : (*this)(1, 1, 0, 0) = (*this)(0, 0, 1, 1);
532 0 : (*this)(0, 1, 0, 1) = 0.5 * ((*this)(0, 0, 0, 0) - (*this)(0, 0, 1, 1));
533 0 : (*this)(1, 0, 0, 1) = (*this)(0, 1, 0, 1);
534 0 : (*this)(0, 1, 1, 0) = (*this)(0, 1, 0, 1);
535 0 : (*this)(1, 0, 1, 0) = (*this)(0, 1, 0, 1);
536 : }
537 : else
538 0 : mooseError("Please provide correct number of inputs for surface RankFourTensorTempl<T> "
539 : "initialization.");
540 0 : }
541 :
542 : template <typename T>
543 : void
544 69 : RankFourTensorTempl<T>::fillFromInputVector(const std::vector<T> & input, FillMethod fill_method)
545 : {
546 :
547 69 : switch (fill_method)
548 : {
549 0 : case antisymmetric:
550 0 : fillAntisymmetricFromInputVector(input);
551 0 : break;
552 2 : case symmetric9:
553 2 : fillSymmetric9FromInputVector(input);
554 2 : break;
555 2 : case symmetric21:
556 2 : fillSymmetric21FromInputVector(input);
557 2 : break;
558 0 : case general_isotropic:
559 0 : fillGeneralIsotropicFromInputVector(input);
560 0 : break;
561 2 : case symmetric_isotropic:
562 2 : fillSymmetricIsotropicFromInputVector(input);
563 2 : break;
564 1 : case symmetric_isotropic_E_nu:
565 1 : fillSymmetricIsotropicEandNuFromInputVector(input);
566 1 : break;
567 0 : case antisymmetric_isotropic:
568 0 : fillAntisymmetricIsotropicFromInputVector(input);
569 0 : break;
570 1 : case axisymmetric_rz:
571 1 : fillAxisymmetricRZFromInputVector(input);
572 1 : break;
573 59 : case general:
574 59 : fillGeneralFromInputVector(input);
575 59 : break;
576 1 : case principal:
577 1 : fillPrincipalFromInputVector(input);
578 1 : break;
579 1 : case orthotropic:
580 1 : fillGeneralOrthotropicFromInputVector(input);
581 1 : break;
582 0 : default:
583 0 : mooseError("fillFromInputVector called with unknown fill_method of ", fill_method);
584 : }
585 69 : }
586 :
587 : template <typename T>
588 : void
589 0 : RankFourTensorTempl<T>::fillAntisymmetricFromInputVector(const std::vector<T> & input)
590 : {
591 0 : if (input.size() != 6)
592 0 : mooseError(
593 : "To use fillAntisymmetricFromInputVector, your input must have size 6. Yours has size ",
594 0 : input.size());
595 :
596 0 : zero();
597 :
598 0 : (*this)(0, 1, 0, 1) = input[0]; // B1212
599 0 : (*this)(0, 1, 0, 2) = input[1]; // B1213
600 0 : (*this)(0, 1, 1, 2) = input[2]; // B1223
601 :
602 0 : (*this)(0, 2, 0, 2) = input[3]; // B1313
603 0 : (*this)(0, 2, 1, 2) = input[4]; // B1323
604 :
605 0 : (*this)(1, 2, 1, 2) = input[5]; // B2323
606 :
607 : // symmetry on the two pairs
608 0 : (*this)(0, 2, 0, 1) = (*this)(0, 1, 0, 2);
609 0 : (*this)(1, 2, 0, 1) = (*this)(0, 1, 1, 2);
610 0 : (*this)(1, 2, 0, 2) = (*this)(0, 2, 1, 2);
611 : // have now got the upper parts of vals[0][1], vals[0][2] and vals[1][2]
612 :
613 : // fill in from antisymmetry relations
614 0 : for (auto i : make_range(N))
615 0 : for (auto j : make_range(N))
616 : {
617 0 : (*this)(0, 1, j, i) = -(*this)(0, 1, i, j);
618 0 : (*this)(0, 2, j, i) = -(*this)(0, 2, i, j);
619 0 : (*this)(1, 2, j, i) = -(*this)(1, 2, i, j);
620 : }
621 : // have now got all of vals[0][1], vals[0][2] and vals[1][2]
622 :
623 : // fill in from antisymmetry relations
624 0 : for (auto i : make_range(N))
625 0 : for (auto j : make_range(N))
626 : {
627 0 : (*this)(1, 0, i, j) = -(*this)(0, 1, i, j);
628 0 : (*this)(2, 0, i, j) = -(*this)(0, 2, i, j);
629 0 : (*this)(2, 1, i, j) = -(*this)(1, 2, i, j);
630 : }
631 0 : }
632 :
633 : template <typename T>
634 : void
635 0 : RankFourTensorTempl<T>::fillGeneralIsotropicFromInputVector(const std::vector<T> & input)
636 : {
637 0 : if (input.size() != 3)
638 0 : mooseError("To use fillGeneralIsotropicFromInputVector, your input must have size 3. Yours "
639 : "has size ",
640 0 : input.size());
641 :
642 0 : fillGeneralIsotropic(input[0], input[1], input[2]);
643 0 : }
644 :
645 : template <typename T>
646 : void
647 0 : RankFourTensorTempl<T>::fillGeneralIsotropic(const T & i0, const T & i1, const T & i2)
648 : {
649 0 : for (auto i : make_range(N))
650 0 : for (auto j : make_range(N))
651 0 : for (auto k : make_range(N))
652 0 : for (auto l : make_range(N))
653 : {
654 0 : (*this)(i, j, k, l) = i0 * Real(i == j) * Real(k == l) +
655 0 : i1 * Real(i == k) * Real(j == l) + i1 * Real(i == l) * Real(j == k);
656 0 : for (auto m : make_range(N))
657 0 : (*this)(i, j, k, l) +=
658 0 : i2 * Real(PermutationTensor::eps(i, j, m)) * Real(PermutationTensor::eps(k, l, m));
659 : }
660 0 : }
661 :
662 : template <typename T>
663 : void
664 0 : RankFourTensorTempl<T>::fillAntisymmetricIsotropicFromInputVector(const std::vector<T> & input)
665 : {
666 0 : if (input.size() != 1)
667 0 : mooseError("To use fillAntisymmetricIsotropicFromInputVector, your input must have size 1. "
668 : "Yours has size ",
669 0 : input.size());
670 :
671 0 : fillGeneralIsotropic(0.0, 0.0, input[0]);
672 0 : }
673 :
674 : template <typename T>
675 : void
676 0 : RankFourTensorTempl<T>::fillAntisymmetricIsotropic(const T & i0)
677 : {
678 0 : fillGeneralIsotropic(0.0, 0.0, i0);
679 0 : }
680 :
681 : template <typename T>
682 : void
683 2 : RankFourTensorTempl<T>::fillSymmetricIsotropicFromInputVector(const std::vector<T> & input)
684 : {
685 : mooseAssert(input.size() == 2,
686 : "To use fillSymmetricIsotropicFromInputVector, your input must have size 2.");
687 2 : fillSymmetricIsotropic(input[0], input[1]);
688 2 : }
689 :
690 : template <typename T>
691 : void
692 3 : RankFourTensorTempl<T>::fillSymmetricIsotropic(const T & lambda, const T & G)
693 : {
694 : // clang-format off
695 6 : fillSymmetric21FromInputVector(std::array<T,21>
696 3 : {{lambda + 2.0 * G, lambda, lambda, 0.0, 0.0, 0.0,
697 3 : lambda + 2.0 * G, lambda, 0.0, 0.0, 0.0,
698 3 : lambda + 2.0 * G, 0.0, 0.0, 0.0,
699 0 : G, 0.0, 0.0,
700 0 : G, 0.0,
701 : G}});
702 : // clang-format on
703 3 : }
704 :
705 : template <typename T>
706 : void
707 1 : RankFourTensorTempl<T>::fillSymmetricIsotropicEandNuFromInputVector(const std::vector<T> & input)
708 : {
709 1 : if (input.size() != 2)
710 0 : mooseError(
711 : "To use fillSymmetricIsotropicEandNuFromInputVector, your input must have size 2. Yours "
712 : "has size ",
713 0 : input.size());
714 :
715 1 : fillSymmetricIsotropicEandNu(input[0], input[1]);
716 1 : }
717 :
718 : template <typename T>
719 : void
720 1 : RankFourTensorTempl<T>::fillSymmetricIsotropicEandNu(const T & E, const T & nu)
721 : {
722 : // Calculate lambda and the shear modulus from the given young's modulus and poisson's ratio
723 1 : const T & lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
724 1 : const T & G = E / (2.0 * (1.0 + nu));
725 :
726 1 : fillSymmetricIsotropic(lambda, G);
727 1 : }
728 :
729 : template <typename T>
730 : void
731 1 : RankFourTensorTempl<T>::fillAxisymmetricRZFromInputVector(const std::vector<T> & input)
732 : {
733 : mooseAssert(input.size() == 5,
734 : "To use fillAxisymmetricRZFromInputVector, your input must have size 5.");
735 :
736 : // C1111 C1122 C1133 0 0 0
737 : // C2222 C2233=C1133 0 0 0
738 : // C3333 0 0 0
739 : // C2323 0 0
740 : // C3131=C2323 0
741 : // C1212
742 : // clang-format off
743 10 : fillSymmetric21FromInputVector(std::array<T,21>
744 1 : {{input[0],input[1],input[2], 0.0, 0.0, 0.0,
745 1 : input[0],input[2], 0.0, 0.0, 0.0,
746 1 : input[3], 0.0, 0.0, 0.0,
747 1 : input[4], 0.0, 0.0,
748 1 : input[4], 0.0,
749 1 : (input[0] - input[1]) * 0.5}});
750 : // clang-format on
751 1 : }
752 :
753 : template <typename T>
754 : void
755 59 : RankFourTensorTempl<T>::fillGeneralFromInputVector(const std::vector<T> & input)
756 : {
757 59 : if (input.size() != 81)
758 0 : mooseError("To use fillGeneralFromInputVector, your input must have size 81. Yours has size ",
759 0 : input.size());
760 :
761 4838 : for (auto i : make_range(N4))
762 4779 : _vals[i] = input[i];
763 59 : }
764 :
765 : template <typename T>
766 : void
767 1 : RankFourTensorTempl<T>::fillPrincipalFromInputVector(const std::vector<T> & input)
768 : {
769 1 : if (input.size() != 9)
770 0 : mooseError("To use fillPrincipalFromInputVector, your input must have size 9. Yours has size ",
771 0 : input.size());
772 :
773 1 : zero();
774 :
775 1 : (*this)(0, 0, 0, 0) = input[0];
776 1 : (*this)(0, 0, 1, 1) = input[1];
777 1 : (*this)(0, 0, 2, 2) = input[2];
778 1 : (*this)(1, 1, 0, 0) = input[3];
779 1 : (*this)(1, 1, 1, 1) = input[4];
780 1 : (*this)(1, 1, 2, 2) = input[5];
781 1 : (*this)(2, 2, 0, 0) = input[6];
782 1 : (*this)(2, 2, 1, 1) = input[7];
783 1 : (*this)(2, 2, 2, 2) = input[8];
784 1 : }
785 :
786 : template <typename T>
787 : void
788 1 : RankFourTensorTempl<T>::fillGeneralOrthotropicFromInputVector(const std::vector<T> & input)
789 : {
790 1 : if (input.size() != 12)
791 0 : mooseError("To use fillGeneralOrhotropicFromInputVector, your input must have size 12. Yours "
792 : "has size ",
793 0 : input.size());
794 :
795 1 : const T & Ea = input[0];
796 1 : const T & Eb = input[1];
797 1 : const T & Ec = input[2];
798 1 : const T & Gab = input[3];
799 1 : const T & Gbc = input[4];
800 1 : const T & Gca = input[5];
801 1 : const T & nuba = input[6];
802 1 : const T & nuca = input[7];
803 1 : const T & nucb = input[8];
804 1 : const T & nuab = input[9];
805 1 : const T & nuac = input[10];
806 1 : const T & nubc = input[11];
807 :
808 : // Input must satisfy constraints.
809 1 : bool preserve_symmetry = MooseUtils::absoluteFuzzyEqual(nuab * Eb, nuba * Ea) &&
810 2 : MooseUtils::absoluteFuzzyEqual(nuca * Ea, nuac * Ec) &&
811 2 : MooseUtils::absoluteFuzzyEqual(nubc * Ec, nucb * Eb);
812 :
813 1 : if (!preserve_symmetry)
814 0 : mooseError("Orthotropic elasticity tensor input is not consistent with symmetry requirements. "
815 : "Check input for accuracy");
816 :
817 1 : unsigned int ntens = N * (N + 1) / 2;
818 :
819 1 : std::vector<T> mat;
820 1 : mat.assign(ntens * ntens, 0);
821 :
822 1 : T k = 1 - nuab * nuba - nubc * nucb - nuca * nuac - nuab * nubc * nuca - nuba * nucb * nuac;
823 :
824 1 : bool is_positive_definite =
825 1 : (k > 0) && (1 - nubc * nucb) > 0 && (1 - nuac * nuca) > 0 && (1 - nuab * nuba) > 0;
826 1 : if (!is_positive_definite)
827 0 : mooseError("Orthotropic elasticity tensor input is not positive definite. Check input for "
828 : "accuracy");
829 :
830 1 : mat[0] = Ea * (1 - nubc * nucb) / k;
831 1 : mat[1] = Ea * (nubc * nuca + nuba) / k;
832 1 : mat[2] = Ea * (nuba * nucb + nuca) / k;
833 :
834 1 : mat[6] = Eb * (nuac * nucb + nuab) / k;
835 1 : mat[7] = Eb * (1 - nuac * nuca) / k;
836 1 : mat[8] = Eb * (nuab * nuca + nucb) / k;
837 :
838 1 : mat[12] = Ec * (nuab * nubc + nuac) / k;
839 1 : mat[13] = Ec * (nuac * nuba + nubc) / k;
840 1 : mat[14] = Ec * (1 - nuab * nuba) / k;
841 :
842 1 : mat[21] = 2 * Gab;
843 1 : mat[28] = 2 * Gca;
844 1 : mat[35] = 2 * Gbc;
845 :
846 : // Switching from Voigt to fourth order tensor
847 : // Copied from existing code (invSymm)
848 1 : int nskip = N - 1;
849 :
850 1 : unsigned int index = 0;
851 4 : for (auto i : make_range(N))
852 12 : for (auto j : make_range(N))
853 36 : for (auto k : make_range(N))
854 108 : for (auto l : make_range(N))
855 : {
856 81 : if (i == j)
857 27 : (*this)._vals[index] =
858 27 : k == l ? mat[i * ntens + k] : mat[i * ntens + k + nskip + l] / 2.0;
859 : else
860 90 : (*this)._vals[index] = k == l ? mat[(nskip + i + j) * ntens + k]
861 36 : : mat[(nskip + i + j) * ntens + k + nskip + l] / 2.0;
862 81 : index++;
863 : }
864 1 : }
865 :
866 : template <typename T>
867 : RankTwoTensorTempl<T>
868 0 : RankFourTensorTempl<T>::innerProductTranspose(const RankTwoTensorTempl<T> & b) const
869 : {
870 0 : RankTwoTensorTempl<T> result;
871 :
872 0 : unsigned int index = 0;
873 0 : for (unsigned int ij = 0; ij < N2; ++ij)
874 : {
875 0 : T bb = b._coords[ij];
876 0 : for (unsigned int kl = 0; kl < N2; ++kl)
877 0 : result._coords[kl] += _vals[index++] * bb;
878 : }
879 :
880 0 : return result;
881 0 : }
882 :
883 : template <typename T>
884 : T
885 0 : RankFourTensorTempl<T>::contractionIj(unsigned int i,
886 : unsigned int j,
887 : const RankTwoTensorTempl<T> & M) const
888 : {
889 0 : T val = 0;
890 0 : for (unsigned int k = 0; k < N; k++)
891 0 : for (unsigned int l = 0; l < N; l++)
892 0 : val += (*this)(i, j, k, l) * M(k, l);
893 :
894 0 : return val;
895 0 : }
896 :
897 : template <typename T>
898 : T
899 0 : RankFourTensorTempl<T>::contractionKl(unsigned int k,
900 : unsigned int l,
901 : const RankTwoTensorTempl<T> & M) const
902 : {
903 0 : T val = 0;
904 0 : for (unsigned int i = 0; i < N; i++)
905 0 : for (unsigned int j = 0; j < N; j++)
906 0 : val += (*this)(i, j, k, l) * M(i, j);
907 :
908 0 : return val;
909 0 : }
910 :
911 : template <typename T>
912 : T
913 1 : RankFourTensorTempl<T>::sum3x3() const
914 : {
915 : // used in the volumetric locking correction
916 1 : T sum = 0;
917 4 : for (auto i : make_range(N))
918 12 : for (auto j : make_range(N))
919 9 : sum += (*this)(i, i, j, j);
920 1 : return sum;
921 0 : }
922 :
923 : template <typename T>
924 : libMesh::VectorValue<T>
925 1 : RankFourTensorTempl<T>::sum3x1() const
926 : {
927 : // used for volumetric locking correction
928 1 : libMesh::VectorValue<T> a(3);
929 1 : a(0) = (*this)(0, 0, 0, 0) + (*this)(0, 0, 1, 1) + (*this)(0, 0, 2, 2); // C0000 + C0011 + C0022
930 1 : a(1) = (*this)(1, 1, 0, 0) + (*this)(1, 1, 1, 1) + (*this)(1, 1, 2, 2); // C1100 + C1111 + C1122
931 1 : a(2) = (*this)(2, 2, 0, 0) + (*this)(2, 2, 1, 1) + (*this)(2, 2, 2, 2); // C2200 + C2211 + C2222
932 1 : return a;
933 0 : }
934 :
935 : template <typename T>
936 : RankFourTensorTempl<T>
937 0 : RankFourTensorTempl<T>::tripleProductJkl(const RankTwoTensorTempl<T> & A,
938 : const RankTwoTensorTempl<T> & B,
939 : const RankTwoTensorTempl<T> & C) const
940 : {
941 0 : RankFourTensorTempl<T> R;
942 0 : for (unsigned int i = 0; i < N; i++)
943 0 : for (unsigned int j = 0; j < N; j++)
944 0 : for (unsigned int k = 0; k < N; k++)
945 0 : for (unsigned int l = 0; l < N; l++)
946 0 : for (unsigned int m = 0; m < N; m++)
947 0 : for (unsigned int n = 0; n < N; n++)
948 0 : for (unsigned int t = 0; t < N; t++)
949 0 : R(i, j, k, l) += (*this)(i, m, n, t) * A(j, m) * B(k, n) * C(l, t);
950 :
951 0 : return R;
952 0 : }
953 :
954 : template <typename T>
955 : RankFourTensorTempl<T>
956 0 : RankFourTensorTempl<T>::tripleProductIkl(const RankTwoTensorTempl<T> & A,
957 : const RankTwoTensorTempl<T> & B,
958 : const RankTwoTensorTempl<T> & C) const
959 : {
960 0 : RankFourTensorTempl<T> R;
961 0 : for (unsigned int i = 0; i < N; i++)
962 0 : for (unsigned int j = 0; j < N; j++)
963 0 : for (unsigned int k = 0; k < N; k++)
964 0 : for (unsigned int l = 0; l < N; l++)
965 0 : for (unsigned int m = 0; m < N; m++)
966 0 : for (unsigned int n = 0; n < N; n++)
967 0 : for (unsigned int t = 0; t < N; t++)
968 0 : R(i, j, k, l) += (*this)(m, j, n, t) * A(i, m) * B(k, n) * C(l, t);
969 :
970 0 : return R;
971 0 : }
972 :
973 : template <typename T>
974 : RankFourTensorTempl<T>
975 0 : RankFourTensorTempl<T>::tripleProductIjl(const RankTwoTensorTempl<T> & A,
976 : const RankTwoTensorTempl<T> & B,
977 : const RankTwoTensorTempl<T> & C) const
978 : {
979 0 : RankFourTensorTempl<T> R;
980 0 : for (unsigned int i = 0; i < N; i++)
981 0 : for (unsigned int j = 0; j < N; j++)
982 0 : for (unsigned int k = 0; k < N; k++)
983 0 : for (unsigned int l = 0; l < N; l++)
984 0 : for (unsigned int m = 0; m < N; m++)
985 0 : for (unsigned int n = 0; n < N; n++)
986 0 : for (unsigned int t = 0; t < N; t++)
987 0 : R(i, j, k, l) += (*this)(m, n, k, t) * A(i, m) * B(j, n) * C(l, t);
988 :
989 0 : return R;
990 0 : }
991 :
992 : template <typename T>
993 : RankFourTensorTempl<T>
994 0 : RankFourTensorTempl<T>::tripleProductIjk(const RankTwoTensorTempl<T> & A,
995 : const RankTwoTensorTempl<T> & B,
996 : const RankTwoTensorTempl<T> & C) const
997 : {
998 0 : RankFourTensorTempl<T> R;
999 0 : for (unsigned int i = 0; i < N; i++)
1000 0 : for (unsigned int j = 0; j < N; j++)
1001 0 : for (unsigned int k = 0; k < N; k++)
1002 0 : for (unsigned int l = 0; l < N; l++)
1003 0 : for (unsigned int m = 0; m < N; m++)
1004 0 : for (unsigned int n = 0; n < N; n++)
1005 0 : for (unsigned int t = 0; t < N; t++)
1006 0 : R(i, j, k, l) += (*this)(m, n, t, l) * A(i, m) * B(j, n) * C(k, t);
1007 :
1008 0 : return R;
1009 0 : }
1010 :
1011 : template <typename T>
1012 : RankFourTensorTempl<T>
1013 0 : RankFourTensorTempl<T>::singleProductI(const RankTwoTensorTempl<T> & A) const
1014 : {
1015 0 : RankFourTensorTempl<T> R;
1016 :
1017 0 : for (unsigned int i = 0; i < N; i++)
1018 0 : for (unsigned int j = 0; j < N; j++)
1019 0 : for (unsigned int k = 0; k < N; k++)
1020 0 : for (unsigned int l = 0; l < N; l++)
1021 0 : for (unsigned int m = 0; m < N; m++)
1022 0 : R(i, j, k, l) += (*this)(m, j, k, l) * A(i, m);
1023 :
1024 0 : return R;
1025 0 : }
1026 :
1027 : template <typename T>
1028 : RankFourTensorTempl<T>
1029 0 : RankFourTensorTempl<T>::singleProductJ(const RankTwoTensorTempl<T> & A) const
1030 : {
1031 0 : RankFourTensorTempl<T> R;
1032 :
1033 0 : for (unsigned int i = 0; i < N; i++)
1034 0 : for (unsigned int j = 0; j < N; j++)
1035 0 : for (unsigned int k = 0; k < N; k++)
1036 0 : for (unsigned int l = 0; l < N; l++)
1037 0 : for (unsigned int m = 0; m < N; m++)
1038 0 : R(i, j, k, l) += (*this)(i, m, k, l) * A(j, m);
1039 :
1040 0 : return R;
1041 0 : }
1042 :
1043 : template <typename T>
1044 : RankFourTensorTempl<T>
1045 0 : RankFourTensorTempl<T>::singleProductK(const RankTwoTensorTempl<T> & A) const
1046 : {
1047 0 : RankFourTensorTempl<T> R;
1048 :
1049 0 : for (unsigned int i = 0; i < N; i++)
1050 0 : for (unsigned int j = 0; j < N; j++)
1051 0 : for (unsigned int k = 0; k < N; k++)
1052 0 : for (unsigned int l = 0; l < N; l++)
1053 0 : for (unsigned int m = 0; m < N; m++)
1054 0 : R(i, j, k, l) += (*this)(i, j, m, l) * A(k, m);
1055 :
1056 0 : return R;
1057 0 : }
1058 :
1059 : template <typename T>
1060 : RankFourTensorTempl<T>
1061 0 : RankFourTensorTempl<T>::singleProductL(const RankTwoTensorTempl<T> & A) const
1062 : {
1063 0 : RankFourTensorTempl<T> R;
1064 :
1065 0 : for (unsigned int i = 0; i < N; i++)
1066 0 : for (unsigned int j = 0; j < N; j++)
1067 0 : for (unsigned int k = 0; k < N; k++)
1068 0 : for (unsigned int l = 0; l < N; l++)
1069 0 : for (unsigned int m = 0; m < N; m++)
1070 0 : R(i, j, k, l) += (*this)(i, j, k, m) * A(l, m);
1071 :
1072 0 : return R;
1073 0 : }
1074 :
1075 : template <typename T>
1076 : bool
1077 1 : RankFourTensorTempl<T>::isSymmetric() const
1078 : {
1079 3 : for (auto i : make_range(1u, N))
1080 5 : for (auto j : make_range(i))
1081 9 : for (auto k : make_range(1u, N))
1082 15 : for (auto l : make_range(k))
1083 : {
1084 : // minor symmetries
1085 18 : if ((*this)(i, j, k, l) != (*this)(j, i, k, l) ||
1086 9 : (*this)(i, j, k, l) != (*this)(i, j, l, k))
1087 0 : return false;
1088 :
1089 : // major symmetry
1090 9 : if ((*this)(i, j, k, l) != (*this)(k, l, i, j))
1091 0 : return false;
1092 : }
1093 1 : return true;
1094 : }
1095 :
1096 : template <typename T>
1097 : bool
1098 1 : RankFourTensorTempl<T>::isIsotropic() const
1099 : {
1100 : // prerequisite is symmetry
1101 1 : if (!isSymmetric())
1102 0 : return false;
1103 :
1104 : // inspect shear components
1105 1 : const T & mu = (*this)(0, 1, 0, 1);
1106 : // ...diagonal
1107 1 : if ((*this)(1, 2, 1, 2) != mu || (*this)(2, 0, 2, 0) != mu)
1108 0 : return false;
1109 : // ...off-diagonal
1110 1 : if ((*this)(2, 0, 1, 2) != 0.0 || (*this)(0, 1, 1, 2) != 0.0 || (*this)(0, 1, 2, 0) != 0.0)
1111 0 : return false;
1112 :
1113 : // off diagonal blocks in Voigt
1114 4 : for (auto i : make_range(N))
1115 12 : for (auto j : make_range(N))
1116 9 : if (_vals[i * (N3 + N2) + ((j + 1) % N) * N + (j + 2) % N] != 0.0)
1117 0 : return false;
1118 :
1119 : // top left block
1120 1 : const T & K1 = (*this)(0, 0, 0, 0);
1121 1 : const T & K2 = (*this)(0, 0, 1, 1);
1122 1 : if (!MooseUtils::relativeFuzzyEqual(K1 - 4.0 * mu / 3.0, K2 + 2.0 * mu / 3.0))
1123 0 : return false;
1124 1 : if ((*this)(1, 1, 1, 1) != K1 || (*this)(2, 2, 2, 2) != K1)
1125 0 : return false;
1126 3 : for (auto i : make_range(1u, N))
1127 5 : for (auto j : make_range(i))
1128 3 : if ((*this)(i, i, j, j) != K2)
1129 0 : return false;
1130 :
1131 1 : return true;
1132 : }
|
