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 10 : return std::sqrt(l2);
257 0 : }
258 :
259 : template <typename T>
260 : void
261 2 : SymmetricRankFourTensorTempl<T>::print(std::ostream & stm) const
262 : {
263 14 : for (const auto i : make_range(N))
264 : {
265 84 : for (const auto j : make_range(N))
266 72 : stm << std::setw(15) << _vals[i * N + j] << " ";
267 12 : stm << '\n';
268 : }
269 2 : stm << std::flush;
270 2 : }
271 :
272 : template <typename T>
273 : void
274 2 : SymmetricRankFourTensorTempl<T>::printReal(std::ostream & stm) const
275 : {
276 14 : for (const auto i : make_range(N))
277 : {
278 84 : for (const auto j : make_range(N))
279 72 : stm << std::setw(15) << MetaPhysicL::raw_value(_vals[i * N + j]) << " ";
280 12 : stm << '\n';
281 : }
282 2 : stm << std::flush;
283 2 : }
284 :
285 : template <typename T>
286 : SymmetricRankFourTensorTempl<T>
287 2 : SymmetricRankFourTensorTempl<T>::transposeMajor() const
288 : {
289 2 : std::size_t index = 0;
290 2 : SymmetricRankFourTensorTempl<T> ret;
291 14 : for (const auto i : make_range(N))
292 84 : for (const auto j : make_range(N))
293 72 : ret._vals[index++] = _vals[i + N * j];
294 2 : return ret;
295 0 : }
296 :
297 : template <typename T>
298 : void
299 94 : SymmetricRankFourTensorTempl<T>::fillFromInputVector(const std::vector<T> & input,
300 : FillMethod fill_method)
301 : {
302 :
303 94 : switch (fill_method)
304 : {
305 42 : case symmetric9:
306 42 : fillSymmetric9FromInputVector(input);
307 42 : break;
308 42 : case symmetric21:
309 42 : fillSymmetric21FromInputVector(input);
310 42 : break;
311 2 : case symmetric_isotropic:
312 2 : fillSymmetricIsotropicFromInputVector(input);
313 2 : break;
314 2 : case symmetric_isotropic_E_nu:
315 2 : fillSymmetricIsotropicEandNuFromInputVector(input);
316 2 : break;
317 2 : case axisymmetric_rz:
318 2 : fillAxisymmetricRZFromInputVector(input);
319 2 : break;
320 2 : case principal:
321 2 : fillPrincipalFromInputVector(input);
322 2 : break;
323 2 : case orthotropic:
324 2 : fillGeneralOrthotropicFromInputVector(input);
325 2 : break;
326 0 : default:
327 0 : mooseError("fillFromInputVector called with unknown fill_method of ", fill_method);
328 : }
329 94 : }
330 :
331 : template <typename T>
332 : void
333 2 : SymmetricRankFourTensorTempl<T>::fillSymmetricIsotropicFromInputVector(const std::vector<T> & input)
334 : {
335 : mooseAssert(input.size() == 2,
336 : "To use fillSymmetricIsotropicFromInputVector, your input must have size 2.");
337 2 : fillSymmetricIsotropic(input[0], input[1]);
338 2 : }
339 :
340 : template <typename T>
341 : void
342 4 : SymmetricRankFourTensorTempl<T>::fillSymmetricIsotropic(const T & lambda, const T & G)
343 : {
344 : // clang-format off
345 8 : fillSymmetric21FromInputVector(std::array<T,21>
346 4 : {{lambda + 2.0 * G, lambda, lambda, 0.0, 0.0, 0.0,
347 4 : lambda + 2.0 * G, lambda, 0.0, 0.0, 0.0,
348 4 : lambda + 2.0 * G, 0.0, 0.0, 0.0,
349 0 : G, 0.0, 0.0,
350 0 : G, 0.0,
351 : G}});
352 : // clang-format on
353 4 : }
354 :
355 : template <typename T>
356 : void
357 2 : SymmetricRankFourTensorTempl<T>::fillSymmetricIsotropicEandNuFromInputVector(
358 : const std::vector<T> & input)
359 : {
360 2 : if (input.size() != 2)
361 0 : mooseError(
362 : "To use fillSymmetricIsotropicEandNuFromInputVector, your input must have size 2. Yours "
363 : "has size ",
364 0 : input.size());
365 :
366 2 : fillSymmetricIsotropicEandNu(input[0], input[1]);
367 2 : }
368 :
369 : template <typename T>
370 : void
371 2 : SymmetricRankFourTensorTempl<T>::fillSymmetricIsotropicEandNu(const T & E, const T & nu)
372 : {
373 : // Calculate lambda and the shear modulus from the given young's modulus and poisson's ratio
374 2 : const T & lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
375 2 : const T & G = E / (2.0 * (1.0 + nu));
376 :
377 2 : fillSymmetricIsotropic(lambda, G);
378 2 : }
379 :
380 : template <typename T>
381 : void
382 2 : SymmetricRankFourTensorTempl<T>::fillAxisymmetricRZFromInputVector(const std::vector<T> & input)
383 : {
384 : mooseAssert(input.size() == 5,
385 : "To use fillAxisymmetricRZFromInputVector, your input must have size 5.");
386 :
387 : // C1111 C1122 C1133 0 0 0
388 : // C2222 C2233=C1133 0 0 0
389 : // C3333 0 0 0
390 : // C2323 0 0
391 : // C3131=C2323 0
392 : // C1212
393 : // clang-format off
394 20 : fillSymmetric21FromInputVector(std::array<T,21>
395 2 : {{input[0],input[1],input[2], 0.0, 0.0, 0.0,
396 2 : input[0],input[2], 0.0, 0.0, 0.0,
397 2 : input[3], 0.0, 0.0, 0.0,
398 2 : input[4], 0.0, 0.0,
399 2 : input[4], 0.0,
400 2 : (input[0] - input[1]) * 0.5}});
401 : // clang-format on
402 2 : }
403 :
404 : template <typename T>
405 : void
406 2 : SymmetricRankFourTensorTempl<T>::fillPrincipalFromInputVector(const std::vector<T> & input)
407 : {
408 2 : if (input.size() != 9)
409 0 : mooseError("To use fillPrincipalFromInputVector, your input must have size 9. Yours has size ",
410 0 : input.size());
411 :
412 2 : zero();
413 :
414 : // top left block
415 2 : _vals[0] = input[0];
416 2 : _vals[1] = input[1];
417 2 : _vals[2] = input[2];
418 2 : _vals[6] = input[3];
419 2 : _vals[7] = input[4];
420 2 : _vals[8] = input[5];
421 2 : _vals[12] = input[6];
422 2 : _vals[13] = input[7];
423 2 : _vals[14] = input[8];
424 2 : }
425 :
426 : template <typename T>
427 : void
428 2 : SymmetricRankFourTensorTempl<T>::fillGeneralOrthotropicFromInputVector(const std::vector<T> & input)
429 : {
430 : mooseAssert(LIBMESH_DIM == 3, "This method assumes LIBMESH_DIM == 3");
431 2 : if (input.size() != 12)
432 0 : mooseError("To use fillGeneralOrhotropicFromInputVector, your input must have size 12. Yours "
433 : "has size ",
434 0 : input.size());
435 :
436 2 : const T & Ea = input[0];
437 2 : const T & Eb = input[1];
438 2 : const T & Ec = input[2];
439 2 : const T & Gab = input[3];
440 2 : const T & Gbc = input[4];
441 2 : const T & Gca = input[5];
442 2 : const T & nuba = input[6];
443 2 : const T & nuca = input[7];
444 2 : const T & nucb = input[8];
445 2 : const T & nuab = input[9];
446 2 : const T & nuac = input[10];
447 2 : const T & nubc = input[11];
448 :
449 : // Input must satisfy constraints.
450 2 : bool preserve_symmetry = MooseUtils::relativeFuzzyEqual(nuab * Eb, nuba * Ea) &&
451 4 : MooseUtils::relativeFuzzyEqual(nuca * Ea, nuac * Ec) &&
452 4 : MooseUtils::relativeFuzzyEqual(nubc * Ec, nucb * Eb);
453 :
454 2 : if (!preserve_symmetry)
455 0 : mooseError("Orthotropic elasticity tensor input is not consistent with symmetry requirements. "
456 : "Check input for accuracy");
457 :
458 2 : zero();
459 2 : T k = 1 - nuab * nuba - nubc * nucb - nuca * nuac - nuab * nubc * nuca - nuba * nucb * nuac;
460 :
461 2 : bool is_positive_definite =
462 2 : (k > 0) && (1 - nubc * nucb) > 0 && (1 - nuac * nuca) > 0 && (1 - nuab * nuba) > 0;
463 2 : if (!is_positive_definite)
464 0 : mooseError("Orthotropic elasticity tensor input is not positive definite. Check input for "
465 : "accuracy");
466 :
467 2 : _vals[0] = Ea * (1 - nubc * nucb) / k;
468 2 : _vals[1] = Ea * (nubc * nuca + nuba) / k;
469 2 : _vals[2] = Ea * (nuba * nucb + nuca) / k;
470 :
471 2 : _vals[6] = Eb * (nuac * nucb + nuab) / k;
472 2 : _vals[7] = Eb * (1 - nuac * nuca) / k;
473 2 : _vals[8] = Eb * (nuab * nuca + nucb) / k;
474 :
475 2 : _vals[12] = Ec * (nuab * nubc + nuac) / k;
476 2 : _vals[13] = Ec * (nuac * nuba + nubc) / k;
477 2 : _vals[14] = Ec * (1 - nuab * nuba) / k;
478 :
479 2 : _vals[21] = 2 * Gbc;
480 2 : _vals[28] = 2 * Gca;
481 2 : _vals[35] = 2 * Gab;
482 2 : }
483 :
484 : template <typename T>
485 : T
486 2 : SymmetricRankFourTensorTempl<T>::sum3x3() const
487 : {
488 : mooseAssert(LIBMESH_DIM == 3, "This method assumes LIBMESH_DIM == 3");
489 : // summation of Ciijj used in the volumetric locking correction
490 2 : T sum = 0;
491 8 : for (const auto i : make_range(3))
492 24 : for (const auto j : make_range(3))
493 18 : sum += (*this)(i, j);
494 2 : return sum;
495 0 : }
496 :
497 : template <typename T>
498 : libMesh::VectorValue<T>
499 2 : SymmetricRankFourTensorTempl<T>::sum3x1() const
500 : {
501 : mooseAssert(LIBMESH_DIM == 3, "This method assumes LIBMESH_DIM == 3");
502 : // used for volumetric locking correction
503 2 : return libMesh::VectorValue<T>(_vals[0] + _vals[1] + _vals[2],
504 2 : _vals[6] + _vals[7] + _vals[8],
505 4 : _vals[12] + _vals[13] + _vals[14]);
506 : }
507 :
508 : template <typename T>
509 : bool
510 8 : SymmetricRankFourTensorTempl<T>::isSymmetric() const
511 : {
512 44 : for (unsigned int i = 0; i < N; ++i)
513 256 : for (unsigned int j = 0; j < N; ++j)
514 : // major symmetry
515 220 : if (_vals[i + N * j] != _vals[N * i + j])
516 2 : return false;
517 6 : return true;
518 : }
519 :
520 : template <typename T>
521 : bool
522 2 : SymmetricRankFourTensorTempl<T>::isIsotropic() const
523 : {
524 : // prerequisite is symmetry
525 2 : if (!isSymmetric())
526 0 : return false;
527 :
528 : // inspect shear components
529 2 : const T & mu = _vals[35];
530 :
531 : // ...diagonal
532 2 : if (_vals[28] != mu || _vals[21] != mu)
533 0 : return false;
534 :
535 : // ...off-diagonal
536 2 : if (_vals[22] != 0.0 || _vals[23] != 0.0 || _vals[29] != 0.0)
537 0 : return false;
538 :
539 : // off diagonal blocks in Voigt
540 8 : for (const auto i : make_range(3))
541 24 : for (const auto j : make_range(3))
542 18 : if (_vals[3 + i + N * j] != 0.0)
543 0 : return false;
544 :
545 : // top left block
546 2 : const T & K1 = _vals[0];
547 2 : const T & K2 = _vals[1];
548 2 : if (!MooseUtils::relativeFuzzyEqual(K1 - 2.0 * mu / 3.0, K2 + mu / 3.0))
549 0 : return false;
550 2 : if (_vals[7] != K1 || _vals[14] != K1)
551 0 : return false;
552 :
553 6 : for (const auto i : make_range(1, 3))
554 10 : for (const auto j : make_range(i))
555 6 : if (_vals[i + N * j] != K2)
556 0 : return false;
557 :
558 2 : return true;
559 : }
|
