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