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