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 "SymmetricRankFourTensor.h"
13 :
14 : // MOOSE includes
15 : #include "SymmetricRankTwoTensor.h"
16 : #include "RankFourTensor.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 : // C++ includes
28 : #include <iomanip>
29 : #include <ostream>
30 :
31 : namespace MathUtils
32 : {
33 : template <>
34 : void mooseSetToZero<SymmetricRankFourTensorTempl<Real>>(SymmetricRankFourTensorTempl<Real> & v);
35 : template <>
36 : void mooseSetToZero<SymmetricRankFourTensorTempl<ADReal>>(SymmetricRankFourTensorTempl<ADReal> & v);
37 : }
38 :
39 : template <typename T>
40 : MooseEnum
41 0 : SymmetricRankFourTensorTempl<T>::fillMethodEnum()
42 : {
43 : return MooseEnum("symmetric9 symmetric21 symmetric_isotropic symmetric_isotropic_E_nu "
44 0 : "axisymmetric_rz principal orthotropic");
45 : }
46 :
47 : template <typename T>
48 7841 : SymmetricRankFourTensorTempl<T>::SymmetricRankFourTensorTempl()
49 : {
50 : mooseAssert(Ndim == 3,
51 : "SymmetricRankFourTensorTempl<T> is designed to only work in 3 dimensions.");
52 7841 : zero();
53 7841 : }
54 :
55 : template <typename T>
56 22 : SymmetricRankFourTensorTempl<T>::SymmetricRankFourTensorTempl(const InitMethod init)
57 : {
58 22 : switch (init)
59 : {
60 18 : case initNone:
61 18 : break;
62 :
63 2 : case initIdentity:
64 2 : zero();
65 8 : for (const auto i : make_range(Ndim))
66 6 : (*this)(i, i) = 1.0;
67 2 : break;
68 :
69 2 : case initIdentitySymmetricFour:
70 2 : zero();
71 14 : for (const auto i : make_range(N))
72 12 : (*this)(i, i) = 1.0;
73 2 : break;
74 :
75 0 : default:
76 0 : mooseError("Unknown SymmetricRankFourTensorTempl<T> initialization pattern.");
77 : }
78 22 : }
79 :
80 : template <typename T>
81 0 : SymmetricRankFourTensorTempl<T>::SymmetricRankFourTensorTempl(const RankFourTensorTempl<T> & t)
82 : {
83 0 : for (const auto a : make_range(N))
84 0 : for (const auto b : make_range(N))
85 : {
86 0 : const auto & idx = full_index[a][b];
87 0 : auto i = idx[0];
88 0 : auto j = idx[1];
89 0 : auto k = idx[2];
90 0 : auto l = idx[3];
91 0 : (*this)(a, b) =
92 0 : (t(i, j, k, l) + t(j, i, l, k) + t(j, i, k, l) + t(i, j, l, k)) / 4 * mandelFactor(a, b);
93 : }
94 0 : }
95 :
96 : template <typename T>
97 70 : SymmetricRankFourTensorTempl<T>::operator RankFourTensorTempl<T>()
98 : {
99 70 : auto & q = *this;
100 70 : RankFourTensorTempl<T> r;
101 490 : for (const auto a : make_range(N))
102 2940 : for (const auto b : make_range(N))
103 : {
104 2520 : const auto i = full_index[a][b][0];
105 2520 : const auto j = full_index[a][b][1];
106 2520 : const auto k = full_index[a][b][2];
107 2520 : const auto l = full_index[a][b][3];
108 :
109 : // Rijkl = Rjikl = Rijlk = Rjilk
110 2520 : r(i, j, k, l) = q(a, b) / mandelFactor(a, b);
111 2520 : r(j, i, k, l) = q(a, b) / mandelFactor(a, b);
112 2520 : r(i, j, l, k) = q(a, b) / mandelFactor(a, b);
113 2520 : r(j, i, l, k) = q(a, b) / mandelFactor(a, b);
114 : }
115 :
116 70 : return r;
117 0 : }
118 :
119 : template <typename T>
120 14 : SymmetricRankFourTensorTempl<T>::SymmetricRankFourTensorTempl(const std::vector<T> & input,
121 0 : FillMethod fill_method)
122 : {
123 14 : fillFromInputVector(input, fill_method);
124 14 : }
125 :
126 : template <typename T>
127 : void
128 7891 : SymmetricRankFourTensorTempl<T>::zero()
129 : {
130 7891 : std::fill(_vals.begin(), _vals.end(), 0.0);
131 7891 : }
132 :
133 : template <typename T>
134 : SymmetricRankFourTensorTempl<T>
135 4 : SymmetricRankFourTensorTempl<T>::rotationMatrix(const TypeTensor<T> & R)
136 : {
137 4 : SymmetricRankFourTensorTempl<T> M(initNone);
138 : const static std::array<std::size_t, 3> a = {{1, 0, 0}};
139 : const static std::array<std::size_t, 3> b = {{2, 2, 1}};
140 16 : for (std::size_t i = 0; i < 3; ++i)
141 48 : for (std::size_t j = 0; j < 3; ++j)
142 : {
143 36 : M(i, j) = R(i, j) * R(i, j);
144 36 : M(i + 3, j) = MathUtils::sqrt2 * R((i + 1) % 3, j) * R((i + 2) % 3, j);
145 36 : M(j, i + 3) = MathUtils::sqrt2 * R(j, (i + 1) % 3) * R(j, (i + 2) % 3);
146 36 : M(i + 3, j + 3) = R(a[i], a[j]) * R(b[i], b[j]) + R(a[i], b[j]) * R(b[i], a[j]);
147 : }
148 4 : return M;
149 0 : }
150 :
151 : template <typename T>
152 : void
153 2 : SymmetricRankFourTensorTempl<T>::rotate(const TypeTensor<T> & R)
154 : {
155 : // build 6x6 rotation matrix
156 2 : auto M = SymmetricRankFourTensorTempl<T>::rotationMatrix(R);
157 :
158 : // rotate tensor
159 2 : (*this) = M * (*this) * M.transposeMajor();
160 2 : }
161 :
162 : template <typename T>
163 : SymmetricRankFourTensorTempl<T> &
164 2 : SymmetricRankFourTensorTempl<T>::operator*=(const T & a)
165 : {
166 74 : for (const auto i : make_range(N2))
167 72 : _vals[i] *= a;
168 2 : return *this;
169 : }
170 :
171 : template <typename T>
172 : SymmetricRankFourTensorTempl<T> &
173 2 : SymmetricRankFourTensorTempl<T>::operator/=(const T & a)
174 : {
175 74 : for (const auto i : make_range(N2))
176 72 : _vals[i] /= a;
177 2 : return *this;
178 : }
179 :
180 : template <typename T>
181 : SymmetricRankFourTensorTempl<T> &
182 14 : SymmetricRankFourTensorTempl<T>::operator+=(const SymmetricRankFourTensorTempl<T> & a)
183 : {
184 518 : for (const auto i : make_range(N2))
185 504 : _vals[i] += a._vals[i];
186 14 : return *this;
187 : }
188 :
189 : template <typename T>
190 : template <typename T2>
191 : auto
192 2 : SymmetricRankFourTensorTempl<T>::operator+(const SymmetricRankFourTensorTempl<T2> & b) const
193 : -> SymmetricRankFourTensorTempl<decltype(T() + T2())>
194 : {
195 2 : SymmetricRankFourTensorTempl<decltype(T() + T2())> result;
196 74 : for (const auto i : make_range(N2))
197 72 : result._vals[i] = _vals[i] + b._vals[i];
198 2 : return result;
199 0 : }
200 :
201 : template <typename T>
202 : SymmetricRankFourTensorTempl<T> &
203 2 : SymmetricRankFourTensorTempl<T>::operator-=(const SymmetricRankFourTensorTempl<T> & a)
204 : {
205 74 : for (const auto i : make_range(N2))
206 72 : _vals[i] -= a._vals[i];
207 2 : return *this;
208 : }
209 :
210 : template <typename T>
211 : template <typename T2>
212 : auto
213 6 : SymmetricRankFourTensorTempl<T>::operator-(const SymmetricRankFourTensorTempl<T2> & b) const
214 : -> SymmetricRankFourTensorTempl<decltype(T() - T2())>
215 : {
216 6 : SymmetricRankFourTensorTempl<decltype(T() - T2())> result;
217 222 : for (const auto i : make_range(N2))
218 216 : result._vals[i] = _vals[i] - b._vals[i];
219 6 : return result;
220 0 : }
221 :
222 : template <typename T>
223 : SymmetricRankFourTensorTempl<T>
224 2 : SymmetricRankFourTensorTempl<T>::operator-() const
225 : {
226 2 : SymmetricRankFourTensorTempl<T> result;
227 74 : for (const auto i : make_range(N2))
228 72 : result._vals[i] = -_vals[i];
229 2 : return result;
230 0 : }
231 :
232 : template <typename T>
233 : template <typename T2>
234 : auto
235 6 : SymmetricRankFourTensorTempl<T>::operator*(const SymmetricRankFourTensorTempl<T2> & b) const
236 : -> SymmetricRankFourTensorTempl<decltype(T() * T2())>
237 : {
238 : typedef decltype(T() * T2()) ValueType;
239 6 : SymmetricRankFourTensorTempl<ValueType> result;
240 :
241 42 : for (const auto i : make_range(N))
242 252 : for (const auto j : make_range(N))
243 1512 : for (const auto p : make_range(N))
244 1296 : result(i, j) += (*this)(i, p) * b(p, j);
245 :
246 6 : return result;
247 0 : }
248 :
249 : template <typename T>
250 : T
251 10 : SymmetricRankFourTensorTempl<T>::L2norm() const
252 : {
253 10 : T l2 = Utility::pow<2>(_vals[0]);
254 360 : for (const auto i : make_range(1u, N2))
255 350 : l2 += Utility::pow<2>(_vals[i]);
256 : using std::sqrt;
257 10 : return sqrt(l2);
258 0 : }
259 :
260 : template <typename T>
261 : void
262 2 : SymmetricRankFourTensorTempl<T>::print(std::ostream & stm) const
263 : {
264 14 : for (const auto i : make_range(N))
265 : {
266 84 : for (const auto j : make_range(N))
267 72 : stm << std::setw(15) << _vals[i * N + j] << " ";
268 12 : stm << '\n';
269 : }
270 2 : stm << std::flush;
271 2 : }
272 :
273 : template <typename T>
274 : void
275 2 : SymmetricRankFourTensorTempl<T>::printReal(std::ostream & stm) const
276 : {
277 14 : for (const auto i : make_range(N))
278 : {
279 84 : for (const auto j : make_range(N))
280 72 : stm << std::setw(15) << MetaPhysicL::raw_value(_vals[i * N + j]) << " ";
281 12 : stm << '\n';
282 : }
283 2 : stm << std::flush;
284 2 : }
285 :
286 : template <typename T>
287 : SymmetricRankFourTensorTempl<T>
288 2 : SymmetricRankFourTensorTempl<T>::transposeMajor() const
289 : {
290 2 : std::size_t index = 0;
291 2 : SymmetricRankFourTensorTempl<T> ret;
292 14 : for (const auto i : make_range(N))
293 84 : for (const auto j : make_range(N))
294 72 : ret._vals[index++] = _vals[i + N * j];
295 2 : return ret;
296 0 : }
297 :
298 : template <typename T>
299 : void
300 94 : SymmetricRankFourTensorTempl<T>::fillFromInputVector(const std::vector<T> & input,
301 : FillMethod fill_method)
302 : {
303 :
304 94 : switch (fill_method)
305 : {
306 42 : case symmetric9:
307 42 : fillSymmetric9FromInputVector(input);
308 42 : break;
309 42 : case symmetric21:
310 42 : fillSymmetric21FromInputVector(input);
311 42 : break;
312 2 : case symmetric_isotropic:
313 2 : fillSymmetricIsotropicFromInputVector(input);
314 2 : break;
315 2 : case symmetric_isotropic_E_nu:
316 2 : fillSymmetricIsotropicEandNuFromInputVector(input);
317 2 : break;
318 2 : case axisymmetric_rz:
319 2 : fillAxisymmetricRZFromInputVector(input);
320 2 : break;
321 2 : case principal:
322 2 : fillPrincipalFromInputVector(input);
323 2 : break;
324 2 : case orthotropic:
325 2 : fillGeneralOrthotropicFromInputVector(input);
326 2 : break;
327 0 : default:
328 0 : mooseError("fillFromInputVector called with unknown fill_method of ", fill_method);
329 : }
330 94 : }
331 :
332 : template <typename T>
333 : void
334 2 : SymmetricRankFourTensorTempl<T>::fillSymmetricIsotropicFromInputVector(const std::vector<T> & input)
335 : {
336 : mooseAssert(input.size() == 2,
337 : "To use fillSymmetricIsotropicFromInputVector, your input must have size 2.");
338 2 : fillSymmetricIsotropic(input[0], input[1]);
339 2 : }
340 :
341 : template <typename T>
342 : void
343 4 : SymmetricRankFourTensorTempl<T>::fillSymmetricIsotropic(const T & lambda, const T & G)
344 : {
345 : // clang-format off
346 8 : fillSymmetric21FromInputVector(std::array<T,21>
347 4 : {{lambda + 2.0 * G, lambda, lambda, 0.0, 0.0, 0.0,
348 4 : lambda + 2.0 * G, lambda, 0.0, 0.0, 0.0,
349 4 : lambda + 2.0 * G, 0.0, 0.0, 0.0,
350 0 : G, 0.0, 0.0,
351 0 : G, 0.0,
352 : G}});
353 : // clang-format on
354 4 : }
355 :
356 : template <typename T>
357 : void
358 2 : SymmetricRankFourTensorTempl<T>::fillSymmetricIsotropicEandNuFromInputVector(
359 : const std::vector<T> & input)
360 : {
361 2 : if (input.size() != 2)
362 0 : mooseError(
363 : "To use fillSymmetricIsotropicEandNuFromInputVector, your input must have size 2. Yours "
364 : "has size ",
365 0 : input.size());
366 :
367 2 : fillSymmetricIsotropicEandNu(input[0], input[1]);
368 2 : }
369 :
370 : template <typename T>
371 : void
372 2 : SymmetricRankFourTensorTempl<T>::fillSymmetricIsotropicEandNu(const T & E, const T & nu)
373 : {
374 : // Calculate lambda and the shear modulus from the given young's modulus and poisson's ratio
375 2 : const T & lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
376 2 : const T & G = E / (2.0 * (1.0 + nu));
377 :
378 2 : fillSymmetricIsotropic(lambda, G);
379 2 : }
380 :
381 : template <typename T>
382 : void
383 2 : SymmetricRankFourTensorTempl<T>::fillAxisymmetricRZFromInputVector(const std::vector<T> & input)
384 : {
385 : mooseAssert(input.size() == 5,
386 : "To use fillAxisymmetricRZFromInputVector, your input must have size 5.");
387 :
388 : // C1111 C1122 C1133 0 0 0
389 : // C2222 C2233=C1133 0 0 0
390 : // C3333 0 0 0
391 : // C2323 0 0
392 : // C3131=C2323 0
393 : // C1212
394 : // clang-format off
395 20 : fillSymmetric21FromInputVector(std::array<T,21>
396 2 : {{input[0],input[1],input[2], 0.0, 0.0, 0.0,
397 2 : input[0],input[2], 0.0, 0.0, 0.0,
398 2 : input[3], 0.0, 0.0, 0.0,
399 2 : input[4], 0.0, 0.0,
400 2 : input[4], 0.0,
401 2 : (input[0] - input[1]) * 0.5}});
402 : // clang-format on
403 2 : }
404 :
405 : template <typename T>
406 : void
407 2 : SymmetricRankFourTensorTempl<T>::fillPrincipalFromInputVector(const std::vector<T> & input)
408 : {
409 2 : if (input.size() != 9)
410 0 : mooseError("To use fillPrincipalFromInputVector, your input must have size 9. Yours has size ",
411 0 : input.size());
412 :
413 2 : zero();
414 :
415 : // top left block
416 2 : _vals[0] = input[0];
417 2 : _vals[1] = input[1];
418 2 : _vals[2] = input[2];
419 2 : _vals[6] = input[3];
420 2 : _vals[7] = input[4];
421 2 : _vals[8] = input[5];
422 2 : _vals[12] = input[6];
423 2 : _vals[13] = input[7];
424 2 : _vals[14] = input[8];
425 2 : }
426 :
427 : template <typename T>
428 : void
429 2 : SymmetricRankFourTensorTempl<T>::fillGeneralOrthotropicFromInputVector(const std::vector<T> & input)
430 : {
431 : mooseAssert(LIBMESH_DIM == 3, "This method assumes LIBMESH_DIM == 3");
432 2 : if (input.size() != 12)
433 0 : mooseError("To use fillGeneralOrhotropicFromInputVector, your input must have size 12. Yours "
434 : "has size ",
435 0 : input.size());
436 :
437 2 : const T & Ea = input[0];
438 2 : const T & Eb = input[1];
439 2 : const T & Ec = input[2];
440 2 : const T & Gab = input[3];
441 2 : const T & Gbc = input[4];
442 2 : const T & Gca = input[5];
443 2 : const T & nuba = input[6];
444 2 : const T & nuca = input[7];
445 2 : const T & nucb = input[8];
446 2 : const T & nuab = input[9];
447 2 : const T & nuac = input[10];
448 2 : const T & nubc = input[11];
449 :
450 : // Input must satisfy constraints.
451 2 : bool preserve_symmetry = MooseUtils::relativeFuzzyEqual(nuab * Eb, nuba * Ea) &&
452 4 : MooseUtils::relativeFuzzyEqual(nuca * Ea, nuac * Ec) &&
453 4 : MooseUtils::relativeFuzzyEqual(nubc * Ec, nucb * Eb);
454 :
455 2 : if (!preserve_symmetry)
456 0 : mooseError("Orthotropic elasticity tensor input is not consistent with symmetry requirements. "
457 : "Check input for accuracy");
458 :
459 2 : zero();
460 2 : T k = 1 - nuab * nuba - nubc * nucb - nuca * nuac - nuab * nubc * nuca - nuba * nucb * nuac;
461 :
462 2 : bool is_positive_definite =
463 2 : (k > 0) && (1 - nubc * nucb) > 0 && (1 - nuac * nuca) > 0 && (1 - nuab * nuba) > 0;
464 2 : if (!is_positive_definite)
465 0 : mooseError("Orthotropic elasticity tensor input is not positive definite. Check input for "
466 : "accuracy");
467 :
468 2 : _vals[0] = Ea * (1 - nubc * nucb) / k;
469 2 : _vals[1] = Ea * (nubc * nuca + nuba) / k;
470 2 : _vals[2] = Ea * (nuba * nucb + nuca) / k;
471 :
472 2 : _vals[6] = Eb * (nuac * nucb + nuab) / k;
473 2 : _vals[7] = Eb * (1 - nuac * nuca) / k;
474 2 : _vals[8] = Eb * (nuab * nuca + nucb) / k;
475 :
476 2 : _vals[12] = Ec * (nuab * nubc + nuac) / k;
477 2 : _vals[13] = Ec * (nuac * nuba + nubc) / k;
478 2 : _vals[14] = Ec * (1 - nuab * nuba) / k;
479 :
480 2 : _vals[21] = 2 * Gbc;
481 2 : _vals[28] = 2 * Gca;
482 2 : _vals[35] = 2 * Gab;
483 2 : }
484 :
485 : template <typename T>
486 : T
487 2 : SymmetricRankFourTensorTempl<T>::sum3x3() const
488 : {
489 : mooseAssert(LIBMESH_DIM == 3, "This method assumes LIBMESH_DIM == 3");
490 : // summation of Ciijj used in the volumetric locking correction
491 2 : T sum = 0;
492 8 : for (const auto i : make_range(3))
493 24 : for (const auto j : make_range(3))
494 18 : sum += (*this)(i, j);
495 2 : return sum;
496 0 : }
497 :
498 : template <typename T>
499 : libMesh::VectorValue<T>
500 2 : SymmetricRankFourTensorTempl<T>::sum3x1() const
501 : {
502 : mooseAssert(LIBMESH_DIM == 3, "This method assumes LIBMESH_DIM == 3");
503 : // used for volumetric locking correction
504 2 : return libMesh::VectorValue<T>(_vals[0] + _vals[1] + _vals[2],
505 2 : _vals[6] + _vals[7] + _vals[8],
506 4 : _vals[12] + _vals[13] + _vals[14]);
507 : }
508 :
509 : template <typename T>
510 : bool
511 8 : SymmetricRankFourTensorTempl<T>::isSymmetric() const
512 : {
513 44 : for (unsigned int i = 0; i < N; ++i)
514 256 : for (unsigned int j = 0; j < N; ++j)
515 : // major symmetry
516 220 : if (_vals[i + N * j] != _vals[N * i + j])
517 2 : return false;
518 6 : return true;
519 : }
520 :
521 : template <typename T>
522 : bool
523 2 : SymmetricRankFourTensorTempl<T>::isIsotropic() const
524 : {
525 : // prerequisite is symmetry
526 2 : if (!isSymmetric())
527 0 : return false;
528 :
529 : // inspect shear components
530 2 : const T & mu = _vals[35];
531 :
532 : // ...diagonal
533 2 : if (_vals[28] != mu || _vals[21] != mu)
534 0 : return false;
535 :
536 : // ...off-diagonal
537 2 : if (_vals[22] != 0.0 || _vals[23] != 0.0 || _vals[29] != 0.0)
538 0 : return false;
539 :
540 : // off diagonal blocks in Voigt
541 8 : for (const auto i : make_range(3))
542 24 : for (const auto j : make_range(3))
543 18 : if (_vals[3 + i + N * j] != 0.0)
544 0 : return false;
545 :
546 : // top left block
547 2 : const T & K1 = _vals[0];
548 2 : const T & K2 = _vals[1];
549 2 : if (!MooseUtils::relativeFuzzyEqual(K1 - 2.0 * mu / 3.0, K2 + mu / 3.0))
550 0 : return false;
551 2 : if (_vals[7] != K1 || _vals[14] != K1)
552 0 : return false;
553 :
554 6 : for (const auto i : make_range(1, 3))
555 10 : for (const auto j : make_range(i))
556 6 : if (_vals[i + N * j] != K2)
557 0 : return false;
558 :
559 2 : return true;
560 : }
|
