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