www.mooseframework.org
RankFourTensor.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 #include "RankFourTensor.h"
8 
9 // MOOSE includes
10 #include "RankTwoTensor.h"
11 #include "MooseEnum.h"
12 #include "MooseException.h"
13 #include "MooseUtils.h"
14 #include "MatrixTools.h"
15 #include "MaterialProperty.h"
16 #include "PermutationTensor.h"
17 
18 #include "libmesh/utility.h"
19 
20 // C++ includes
21 #include <iomanip>
22 #include <ostream>
23 
24 template <>
25 void
27 {
28  v.zero();
29 }
30 
31 template <>
32 void
33 dataStore(std::ostream & stream, RankFourTensor & rft, void * context)
34 {
35  dataStore(stream, rft._vals, context);
36 }
37 
38 template <>
39 void
40 dataLoad(std::istream & stream, RankFourTensor & rft, void * context)
41 {
42  dataLoad(stream, rft._vals, context);
43 }
44 
47 {
48  return MooseEnum("antisymmetric symmetric9 symmetric21 general_isotropic symmetric_isotropic "
49  "symmetric_isotropic_E_nu antisymmetric_isotropic axisymmetric_rz general "
50  "principal");
51 }
52 
54 {
55  mooseAssert(N == 3, "RankFourTensor is currently only tested for 3 dimensions.");
56 
57  unsigned int index = 0;
58  for (unsigned int i = 0; i < N4; ++i)
59  _vals[index++] = 0.0;
60 }
61 
63 {
64  unsigned int index = 0;
65  switch (init)
66  {
67  case initNone:
68  break;
69 
70  case initIdentity:
71  zero();
72  for (unsigned int i = 0; i < N; ++i)
73  (*this)(i, i, i, i) = 1.0;
74  break;
75 
76  case initIdentityFour:
77  for (unsigned int i = 0; i < N; ++i)
78  for (unsigned int j = 0; j < N; ++j)
79  for (unsigned int k = 0; k < N; ++k)
80  for (unsigned int l = 0; l < N; ++l)
81  _vals[index++] = (i == k) && (j == l);
82  break;
83 
85  for (unsigned int i = 0; i < N; ++i)
86  for (unsigned int j = 0; j < N; ++j)
87  for (unsigned int k = 0; k < N; ++k)
88  for (unsigned int l = 0; l < N; ++l)
89  _vals[index++] = 0.5 * ((i == k) && (j == l)) + 0.5 * ((i == l) && (j == k));
90  break;
91 
92  default:
93  mooseError("Unknown RankFourTensor initialization pattern.");
94  }
95 }
96 
97 RankFourTensor::RankFourTensor(const std::vector<Real> & input, FillMethod fill_method)
98 {
99  fillFromInputVector(input, fill_method);
100 }
101 
102 void
104 {
105  for (unsigned int i = 0; i < N4; ++i)
106  _vals[i] = 0.0;
107 }
108 
111 {
112  for (unsigned int i = 0; i < N4; ++i)
113  _vals[i] = a._vals[i];
114  return *this;
115 }
116 
118 {
119  RankTwoTensor result;
120 
121  unsigned int index = 0;
122  for (unsigned int ij = 0; ij < N2; ++ij)
123  {
124  Real tmp = 0;
125  for (unsigned int kl = 0; kl < N2; ++kl)
126  tmp += _vals[index++] * b._vals[kl];
127  result._vals[ij] = tmp;
128  }
129 
130  return result;
131 }
132 
134 {
135  RealTensorValue result;
136 
137  unsigned int index = 0;
138  for (unsigned int i = 0; i < N; ++i)
139  for (unsigned int j = 0; j < N; ++j)
140  for (unsigned int k = 0; k < N; ++k)
141  for (unsigned int l = 0; l < N; ++l)
142  result(i, j) += _vals[index++] * b(k, l);
143 
144  return result;
145 }
146 
148 {
149  RankFourTensor result;
150 
151  for (unsigned int i = 0; i < N4; ++i)
152  result._vals[i] = _vals[i] * b;
153 
154  return result;
155 }
156 
159 {
160  for (unsigned int i = 0; i < N4; ++i)
161  _vals[i] *= a;
162  return *this;
163 }
164 
166 RankFourTensor::operator/(const Real b) const
167 {
168  RankFourTensor result;
169  for (unsigned int i = 0; i < N4; ++i)
170  result._vals[i] = _vals[i] / b;
171  return result;
172 }
173 
176 {
177  for (unsigned int i = 0; i < N4; ++i)
178  _vals[i] /= a;
179  return *this;
180 }
181 
184 {
185  for (unsigned int i = 0; i < N4; ++i)
186  _vals[i] += a._vals[i];
187  return *this;
188 }
189 
192 {
193  RankFourTensor result;
194  for (unsigned int i = 0; i < N4; ++i)
195  result._vals[i] = _vals[i] + b._vals[i];
196  return result;
197 }
198 
201 {
202  for (unsigned int i = 0; i < N4; ++i)
203  _vals[i] -= a._vals[i];
204  return *this;
205 }
206 
209 {
210  RankFourTensor result;
211  for (unsigned int i = 0; i < N4; ++i)
212  result._vals[i] = _vals[i] - b._vals[i];
213  return result;
214 }
215 
218 {
219  RankFourTensor result;
220  for (unsigned int i = 0; i < N4; ++i)
221  result._vals[i] = -_vals[i];
222  return result;
223 }
224 
226 {
227  RankFourTensor result;
228 
229  unsigned int index = 0;
230  unsigned int ij1 = 0;
231  for (unsigned int i = 0; i < N; ++i)
232  {
233  for (unsigned int j = 0; j < N; ++j)
234  {
235  for (unsigned int k = 0; k < N; ++k)
236  {
237  for (unsigned int l = 0; l < N; ++l)
238  {
239  Real sum = 0;
240  for (unsigned int p = 0; p < N; ++p)
241  {
242  unsigned int p1 = N * p;
243  for (unsigned int q = 0; q < N; ++q)
244  sum += _vals[ij1 + p1 + q] * b(p, q, k, l);
245  }
246  result._vals[index++] = sum;
247  }
248  }
249  ij1 += N2;
250  }
251  }
252 
253  return result;
254 }
255 
256 Real
258 {
259  Real l2 = 0;
260 
261  for (unsigned int i = 0; i < N4; ++i)
262  l2 += Utility::pow<2>(_vals[i]);
263 
264  return std::sqrt(l2);
265 }
266 
269 {
270  unsigned int ntens = N * (N + 1) / 2;
271  int nskip = N - 1;
272 
273  RankFourTensor result;
274  std::vector<PetscScalar> mat;
275  mat.assign(ntens * ntens, 0);
276 
277  // We use the LAPACK matrix inversion routine here. Form the matrix
278  //
279  // mat[0] mat[1] mat[2] mat[3] mat[4] mat[5]
280  // mat[6] mat[7] mat[8] mat[9] mat[10] mat[11]
281  // mat[12] mat[13] mat[14] mat[15] mat[16] mat[17]
282  // mat[18] mat[19] mat[20] mat[21] mat[22] mat[23]
283  // mat[24] mat[25] mat[26] mat[27] mat[28] mat[29]
284  // mat[30] mat[31] mat[32] mat[33] mat[34] mat[35]
285  //
286  // This is filled from the indpendent components of C assuming
287  // the symmetry C_ijkl = C_ijlk = C_jikl.
288  //
289  // If there are two rank-four tensors X and Y then the reason for
290  // this filling becomes apparent if we want to calculate
291  // X_ijkl*Y_klmn = Z_ijmn
292  // For denote the "mat" versions of X, Y and Z by x, y and z.
293  // Then
294  // z_ab = x_ac*y_cb
295  // Eg
296  // z_00 = Z_0000 = X_0000*Y_0000 + X_0011*Y_1111 + X_0022*Y_2200 + 2*X_0001*Y_0100 +
297  // 2*X_0002*Y_0200 + 2*X_0012*Y_1200 (the factors of 2 come from the assumed symmetries)
298  // 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 +
299  // 2*X_0002*2*Y_0201 + 2*X_0012*2*Y_1201
300  // 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 +
301  // 2*X_0102*2*Y_0202 + 2*X_0112*2*Y_1202
302  // Finally, we use LAPACK to find x^-1, and put it back into rank-4 tensor form
303  //
304  // mat[0] = C(0,0,0,0)
305  // mat[1] = C(0,0,1,1)
306  // mat[2] = C(0,0,2,2)
307  // mat[3] = C(0,0,0,1)*2
308  // mat[4] = C(0,0,0,2)*2
309  // mat[5] = C(0,0,1,2)*2
310 
311  // mat[6] = C(1,1,0,0)
312  // mat[7] = C(1,1,1,1)
313  // mat[8] = C(1,1,2,2)
314  // mat[9] = C(1,1,0,1)*2
315  // mat[10] = C(1,1,0,2)*2
316  // mat[11] = C(1,1,1,2)*2
317 
318  // mat[12] = C(2,2,0,0)
319  // mat[13] = C(2,2,1,1)
320  // mat[14] = C(2,2,2,2)
321  // mat[15] = C(2,2,0,1)*2
322  // mat[16] = C(2,2,0,2)*2
323  // mat[17] = C(2,2,1,2)*2
324 
325  // mat[18] = C(0,1,0,0)
326  // mat[19] = C(0,1,1,1)
327  // mat[20] = C(0,1,2,2)
328  // mat[21] = C(0,1,0,1)*2
329  // mat[22] = C(0,1,0,2)*2
330  // mat[23] = C(0,1,1,2)*2
331 
332  // mat[24] = C(0,2,0,0)
333  // mat[25] = C(0,2,1,1)
334  // mat[26] = C(0,2,2,2)
335  // mat[27] = C(0,2,0,1)*2
336  // mat[28] = C(0,2,0,2)*2
337  // mat[29] = C(0,2,1,2)*2
338 
339  // mat[30] = C(1,2,0,0)
340  // mat[31] = C(1,2,1,1)
341  // mat[32] = C(1,2,2,2)
342  // mat[33] = C(1,2,0,1)*2
343  // mat[34] = C(1,2,0,2)*2
344  // mat[35] = C(1,2,1,2)*2
345 
346  unsigned int index = 0;
347  for (unsigned int i = 0; i < N; ++i)
348  for (unsigned int j = 0; j < N; ++j)
349  for (unsigned int k = 0; k < N; ++k)
350  for (unsigned int l = 0; l < N; ++l)
351  {
352  if (i == j)
353  mat[k == l ? i * ntens + k : i * ntens + k + nskip + l] += _vals[index];
354  else // i!=j
355  mat[k == l ? (nskip + i + j) * ntens + k : (nskip + i + j) * ntens + k + nskip + l] +=
356  _vals[index]; // note the +=, which results in double-counting and is rectified
357  // below
358  index++;
359  }
360 
361  for (unsigned int i = 3; i < ntens; ++i)
362  for (unsigned int j = 0; j < ntens; ++j)
363  mat[i * ntens + j] /= 2.0; // because of double-counting above
364 
365  // use LAPACK to find the inverse
366  MatrixTools::inverse(mat, ntens);
367 
368  // build the resulting rank-four tensor
369  // using the inverse of the above algorithm
370  index = 0;
371  for (unsigned int i = 0; i < N; ++i)
372  for (unsigned int j = 0; j < N; ++j)
373  for (unsigned int k = 0; k < N; ++k)
374  for (unsigned int l = 0; l < N; ++l)
375  {
376  if (i == j)
377  result._vals[index] =
378  k == l ? mat[i * ntens + k] : mat[i * ntens + k + nskip + l] / 2.0;
379  else // i!=j
380  result._vals[index] = k == l ? mat[(nskip + i + j) * ntens + k]
381  : mat[(nskip + i + j) * ntens + k + nskip + l] / 2.0;
382  index++;
383  }
384 
385  return result;
386 }
387 
388 void
390 {
391  RankFourTensor old = *this;
392 
393  unsigned int index = 0;
394  for (unsigned int i = 0; i < N; ++i)
395  {
396  for (unsigned int j = 0; j < N; ++j)
397  {
398  for (unsigned int k = 0; k < N; ++k)
399  {
400  for (unsigned int l = 0; l < N; ++l)
401  {
402  unsigned int index2 = 0;
403  Real sum = 0.0;
404  for (unsigned int m = 0; m < N; ++m)
405  {
406  const Real a = R(i, m);
407  for (unsigned int n = 0; n < N; ++n)
408  {
409  const Real ab = a * R(j, n);
410  for (unsigned int o = 0; o < N; ++o)
411  {
412  const Real abc = ab * R(k, o);
413  for (unsigned int p = 0; p < N; ++p)
414  sum += abc * R(l, p) * old._vals[index2++];
415  }
416  }
417  }
418  _vals[index++] = sum;
419  }
420  }
421  }
422  }
423 }
424 
425 void
427 {
428  RankFourTensor old = *this;
429 
430  unsigned int index = 0;
431  unsigned int i1 = 0;
432  for (unsigned int i = 0; i < N; ++i)
433  {
434  unsigned int j1 = 0;
435  for (unsigned int j = 0; j < N; ++j)
436  {
437  unsigned int k1 = 0;
438  for (unsigned int k = 0; k < N; ++k)
439  {
440  unsigned int l1 = 0;
441  for (unsigned int l = 0; l < N; ++l)
442  {
443  unsigned int index2 = 0;
444  Real sum = 0.0;
445  for (unsigned int m = 0; m < N; ++m)
446  {
447  const Real a = R._vals[i1 + m];
448  for (unsigned int n = 0; n < N; ++n)
449  {
450  const Real ab = a * R._vals[j1 + n];
451  for (unsigned int o = 0; o < N; ++o)
452  {
453  const Real abc = ab * R._vals[k1 + o];
454  for (unsigned int p = 0; p < N; ++p)
455  sum += abc * R._vals[l1 + p] * old._vals[index2++];
456  }
457  }
458  }
459  _vals[index++] = sum;
460  l1 += N;
461  }
462  k1 += N;
463  }
464  j1 += N;
465  }
466  i1 += N;
467  }
468 }
469 
470 void
471 RankFourTensor::print(std::ostream & stm) const
472 {
473  for (unsigned int i = 0; i < N; ++i)
474  for (unsigned int j = 0; j < N; ++j)
475  {
476  stm << "i = " << i << " j = " << j << '\n';
477  for (unsigned int k = 0; k < N; ++k)
478  {
479  for (unsigned int l = 0; l < N; ++l)
480  stm << std::setw(15) << (*this)(i, j, k, l) << " ";
481 
482  stm << '\n';
483  }
484  }
485 }
486 
489 {
490  RankFourTensor result;
491 
492  unsigned int index = 0;
493  unsigned int i1 = 0;
494  for (unsigned int i = 0; i < N; ++i)
495  {
496  for (unsigned int j = 0; j < N; ++j)
497  {
498  for (unsigned int k = 0; k < N; ++k)
499  {
500  unsigned int ijk1 = k * N3 + i1 + j;
501  for (unsigned int l = 0; l < N; ++l)
502  result._vals[index++] = _vals[ijk1 + l * N2];
503  }
504  }
505  i1 += N;
506  }
507 
508  return result;
509 }
510 
511 void
512 RankFourTensor::surfaceFillFromInputVector(const std::vector<Real> & input)
513 {
514  zero();
515 
516  if (input.size() == 9)
517  {
518  // then fill from vector C_1111, C_1112, C_1122, C_1212, C_1222, C_1211, C_2211, C_2212, C_2222
519  (*this)(0, 0, 0, 0) = input[0];
520  (*this)(0, 0, 0, 1) = input[1];
521  (*this)(0, 0, 1, 1) = input[2];
522  (*this)(0, 1, 0, 1) = input[3];
523  (*this)(0, 1, 1, 1) = input[4];
524  (*this)(0, 1, 0, 0) = input[5];
525  (*this)(1, 1, 0, 0) = input[6];
526  (*this)(1, 1, 0, 1) = input[7];
527  (*this)(1, 1, 1, 1) = input[8];
528 
529  // fill in remainders from C_ijkl = C_ijlk = C_jikl
530  (*this)(0, 0, 1, 0) = (*this)(0, 0, 0, 1);
531  (*this)(0, 1, 1, 0) = (*this)(0, 1, 0, 1);
532  (*this)(1, 0, 0, 0) = (*this)(0, 1, 0, 0);
533  (*this)(1, 0, 0, 1) = (*this)(0, 1, 0, 1);
534  (*this)(1, 0, 1, 1) = (*this)(0, 1, 1, 1);
535  (*this)(1, 0, 0, 0) = (*this)(0, 1, 0, 0);
536  (*this)(1, 1, 1, 0) = (*this)(1, 1, 0, 1);
537  }
538  else if (input.size() == 2)
539  {
540  // only two independent constants, C_1111 and C_1122
541  (*this)(0, 0, 0, 0) = input[0];
542  (*this)(0, 0, 1, 1) = input[1];
543  // use symmetries
544  (*this)(1, 1, 1, 1) = (*this)(0, 0, 0, 0);
545  (*this)(1, 1, 0, 0) = (*this)(0, 0, 1, 1);
546  (*this)(0, 1, 0, 1) = 0.5 * ((*this)(0, 0, 0, 0) - (*this)(0, 0, 1, 1));
547  (*this)(1, 0, 0, 1) = (*this)(0, 1, 0, 1);
548  (*this)(0, 1, 1, 0) = (*this)(0, 1, 0, 1);
549  (*this)(1, 0, 1, 0) = (*this)(0, 1, 0, 1);
550  }
551  else
552  mooseError(
553  "Please provide correct number of inputs for surface RankFourTensor initialization.");
554 }
555 
556 void
557 RankFourTensor::fillFromInputVector(const std::vector<Real> & input, FillMethod fill_method)
558 {
559  zero();
560 
561  switch (fill_method)
562  {
563  case antisymmetric:
565  break;
566  case symmetric9:
567  fillSymmetricFromInputVector(input, false);
568  break;
569  case symmetric21:
570  fillSymmetricFromInputVector(input, true);
571  break;
572  case general_isotropic:
574  break;
575  case symmetric_isotropic:
577  break;
580  break;
583  break;
584  case axisymmetric_rz:
586  break;
587  case general:
589  break;
590  case principal:
592  break;
593  default:
594  mooseError("fillFromInputVector called with unknown fill_method of ", fill_method);
595  }
596 }
597 
598 void
599 RankFourTensor::fillSymmetricFromInputVector(const std::vector<Real> & input, bool all)
600 {
601  if ((all == true && input.size() != 21) || (all == false && input.size() != 9))
602  mooseError("Please check the number of entries in the stiffness input vector.");
603 
604  zero();
605 
606  if (all == true)
607  {
608  (*this)(0, 0, 0, 0) = input[0]; // C1111
609  (*this)(0, 0, 1, 1) = input[1]; // C1122
610  (*this)(0, 0, 2, 2) = input[2]; // C1133
611  (*this)(0, 0, 1, 2) = input[3]; // C1123
612  (*this)(0, 0, 0, 2) = input[4]; // C1113
613  (*this)(0, 0, 0, 1) = input[5]; // C1112
614 
615  (*this)(1, 1, 1, 1) = input[6]; // C2222
616  (*this)(1, 1, 2, 2) = input[7]; // C2233
617  (*this)(1, 1, 1, 2) = input[8]; // C2223
618  (*this)(0, 2, 1, 1) = input[9]; // C2213 //flipped for filling purposes
619  (*this)(0, 1, 1, 1) = input[10]; // C2212 //flipped for filling purposes
620 
621  (*this)(2, 2, 2, 2) = input[11]; // C3333
622  (*this)(1, 2, 2, 2) = input[12]; // C3323 //flipped for filling purposes
623  (*this)(0, 2, 2, 2) = input[13]; // C3313 //flipped for filling purposes
624  (*this)(0, 1, 2, 2) = input[14]; // C3312 //flipped for filling purposes
625 
626  (*this)(1, 2, 1, 2) = input[15]; // C2323
627  (*this)(0, 2, 1, 2) = input[16]; // C2313 //flipped for filling purposes
628  (*this)(0, 1, 1, 2) = input[17]; // C2312 //flipped for filling purposes
629 
630  (*this)(0, 2, 0, 2) = input[18]; // C1313
631  (*this)(0, 1, 0, 2) = input[19]; // C1312 //flipped for filling purposes
632 
633  (*this)(0, 1, 0, 1) = input[20]; // C1212
634  }
635  else
636  {
637  (*this)(0, 0, 0, 0) = input[0]; // C1111
638  (*this)(0, 0, 1, 1) = input[1]; // C1122
639  (*this)(0, 0, 2, 2) = input[2]; // C1133
640  (*this)(1, 1, 1, 1) = input[3]; // C2222
641  (*this)(1, 1, 2, 2) = input[4]; // C2233
642  (*this)(2, 2, 2, 2) = input[5]; // C3333
643  (*this)(1, 2, 1, 2) = input[6]; // C2323
644  (*this)(0, 2, 0, 2) = input[7]; // C1313
645  (*this)(0, 1, 0, 1) = input[8]; // C1212
646  }
647 
648  // fill in from symmetry relations
649  for (unsigned int i = 0; i < N; ++i)
650  for (unsigned int j = 0; j < N; ++j)
651  for (unsigned int k = 0; k < N; ++k)
652  for (unsigned int l = 0; l < N; ++l)
653  (*this)(i, j, l, k) = (*this)(j, i, k, l) = (*this)(j, i, l, k) = (*this)(k, l, i, j) =
654  (*this)(l, k, j, i) = (*this)(k, l, j, i) = (*this)(l, k, i, j) = (*this)(i, j, k, l);
655 }
656 
657 void
658 RankFourTensor::fillAntisymmetricFromInputVector(const std::vector<Real> & input)
659 {
660  if (input.size() != 6)
661  mooseError(
662  "To use fillAntisymmetricFromInputVector, your input must have size 6. Yours has size ",
663  input.size());
664 
665  zero();
666 
667  (*this)(0, 1, 0, 1) = input[0]; // B1212
668  (*this)(0, 1, 0, 2) = input[1]; // B1213
669  (*this)(0, 1, 1, 2) = input[2]; // B1223
670 
671  (*this)(0, 2, 0, 2) = input[3]; // B1313
672  (*this)(0, 2, 1, 2) = input[4]; // B1323
673 
674  (*this)(1, 2, 1, 2) = input[5]; // B2323
675 
676  // symmetry on the two pairs
677  (*this)(0, 2, 0, 1) = (*this)(0, 1, 0, 2);
678  (*this)(1, 2, 0, 1) = (*this)(0, 1, 1, 2);
679  (*this)(1, 2, 0, 2) = (*this)(0, 2, 1, 2);
680  // have now got the upper parts of vals[0][1], vals[0][2] and vals[1][2]
681 
682  // fill in from antisymmetry relations
683  for (unsigned int i = 0; i < N; ++i)
684  for (unsigned int j = 0; j < N; ++j)
685  {
686  (*this)(0, 1, j, i) = -(*this)(0, 1, i, j);
687  (*this)(0, 2, j, i) = -(*this)(0, 2, i, j);
688  (*this)(1, 2, j, i) = -(*this)(1, 2, i, j);
689  }
690  // have now got all of vals[0][1], vals[0][2] and vals[1][2]
691 
692  // fill in from antisymmetry relations
693  for (unsigned int i = 0; i < N; ++i)
694  for (unsigned int j = 0; j < N; ++j)
695  {
696  (*this)(1, 0, i, j) = -(*this)(0, 1, i, j);
697  (*this)(2, 0, i, j) = -(*this)(0, 2, i, j);
698  (*this)(2, 1, i, j) = -(*this)(1, 2, i, j);
699  }
700 }
701 
702 void
704 {
705  if (input.size() != 3)
706  mooseError(
707  "To use fillGeneralIsotropicFromInputVector, your input must have size 3. Yours has size ",
708  input.size());
709 
710  fillGeneralIsotropic(input[0], input[1], input[2]);
711 }
712 
713 void
714 RankFourTensor::fillGeneralIsotropic(Real i0, Real i1, Real i2)
715 {
716  for (unsigned int i = 0; i < N; ++i)
717  for (unsigned int j = 0; j < N; ++j)
718  for (unsigned int k = 0; k < N; ++k)
719  for (unsigned int l = 0; l < N; ++l)
720  {
721  (*this)(i, j, k, l) =
722  i0 * (i == j) * (k == l) + i1 * (i == k) * (j == l) + i1 * (i == l) * (j == k);
723  for (unsigned int m = 0; m < N; ++m)
724  (*this)(i, j, k, l) +=
725  i2 * PermutationTensor::eps(i, j, m) * PermutationTensor::eps(k, l, m);
726  }
727 }
728 
729 void
731 {
732  if (input.size() != 1)
733  mooseError("To use fillAntisymmetricIsotropicFromInputVector, your input must have size 1. "
734  "Yours has size ",
735  input.size());
736  fillGeneralIsotropic(0.0, 0.0, input[0]);
737 }
738 
739 void
741 {
742  fillGeneralIsotropic(0.0, 0.0, i0);
743 }
744 
745 void
747 {
748  if (input.size() != 2)
749  mooseError("To use fillSymmetricIsotropicFromInputVector, your input must have size 2. Yours "
750  "has size ",
751  input.size());
752  fillGeneralIsotropic(input[0], input[1], 0.0);
753 }
754 
755 void
757 {
758  fillGeneralIsotropic(i0, i1, 0.0);
759 }
760 
761 void
763 {
764  if (input.size() != 2)
765  mooseError(
766  "To use fillSymmetricIsotropicEandNuFromInputVector, your input must have size 2. Yours "
767  "has size ",
768  input.size());
769 
770  fillSymmetricIsotropicEandNu(input[0], input[1]);
771 }
772 
773 void
775 {
776  // Calculate lambda and the shear modulus from the given young's modulus and poisson's ratio
777  const Real lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
778  const Real G = E / (2.0 * (1.0 + nu));
779 
780  fillGeneralIsotropic(lambda, G, 0.0);
781 }
782 
783 void
785 {
786  if (input.size() != 5)
787  mooseError("To use fillAxisymmetricRZFromInputVector, your input must have size 5. Your "
788  "vector has size ",
789  input.size());
790 
791  // C1111 C1122 C1133 C2222 C2233=C1133
793  input[1],
794  input[2],
795  input[0],
796  input[2],
797  // C3333 C2323 C3131=C2323 C1212
798  input[3],
799  input[4],
800  input[4],
801  (input[0] - input[1]) * 0.5},
802  false);
803 }
804 
805 void
806 RankFourTensor::fillGeneralFromInputVector(const std::vector<Real> & input)
807 {
808  if (input.size() != 81)
809  mooseError("To use fillGeneralFromInputVector, your input must have size 81. Yours has size ",
810  input.size());
811 
812  for (unsigned int i = 0; i < N4; ++i)
813  _vals[i] = input[i];
814 }
815 
816 void
817 RankFourTensor::fillPrincipalFromInputVector(const std::vector<Real> & input)
818 {
819  if (input.size() != 9)
820  mooseError("To use fillPrincipalFromInputVector, your input must have size 9. Yours has size ",
821  input.size());
822 
823  zero();
824 
825  (*this)(0, 0, 0, 0) = input[0];
826  (*this)(0, 0, 1, 1) = input[1];
827  (*this)(0, 0, 2, 2) = input[2];
828  (*this)(1, 1, 0, 0) = input[3];
829  (*this)(1, 1, 1, 1) = input[4];
830  (*this)(1, 1, 2, 2) = input[5];
831  (*this)(2, 2, 0, 0) = input[6];
832  (*this)(2, 2, 1, 1) = input[7];
833  (*this)(2, 2, 2, 2) = input[8];
834 }
835 
838 {
839  RankTwoTensor result;
840 
841  unsigned int index = 0;
842  for (unsigned int ij = 0; ij < N2; ++ij)
843  {
844  Real bb = b._vals[ij];
845  for (unsigned int kl = 0; kl < N2; ++kl)
846  result._vals[kl] += _vals[index++] * bb;
847  }
848 
849  return result;
850 }
851 
852 Real
854 {
855  // summation of Ciijj for i and j ranging from 0 to 2 - used in the volumetric locking correction
856  return (*this)(0, 0, 0, 0) + (*this)(0, 0, 1, 1) + (*this)(0, 0, 2, 2) + (*this)(1, 1, 0, 0) +
857  (*this)(1, 1, 1, 1) + (*this)(1, 1, 2, 2) + (*this)(2, 2, 0, 0) + (*this)(2, 2, 1, 1) +
858  (*this)(2, 2, 2, 2);
859 }
860 
863 {
864  // used for volumetric locking correction
865  RealGradient a(3);
866  a(0) = (*this)(0, 0, 0, 0) + (*this)(0, 0, 1, 1) + (*this)(0, 0, 2, 2); // C0000 + C0011 + C0022
867  a(1) = (*this)(1, 1, 0, 0) + (*this)(1, 1, 1, 1) + (*this)(1, 1, 2, 2); // C1100 + C1111 + C1122
868  a(2) = (*this)(2, 2, 0, 0) + (*this)(2, 2, 1, 1) + (*this)(2, 2, 2, 2); // C2200 + C2211 + C2222
869  return a;
870 }
871 
872 bool
874 {
875  for (unsigned int i = 1; i < N; ++i)
876  for (unsigned int j = 0; j < i; ++j)
877  for (unsigned int k = 1; k < N; ++k)
878  for (unsigned int l = 0; l < k; ++l)
879  {
880  // minor symmetries
881  if ((*this)(i, j, k, l) != (*this)(j, i, k, l) ||
882  (*this)(i, j, k, l) != (*this)(i, j, l, k))
883  return false;
884 
885  // major symmetry
886  if ((*this)(i, j, k, l) != (*this)(k, l, i, j))
887  return false;
888  }
889  return true;
890 }
891 
892 bool
894 {
895  // prerequisite is symmetry
896  if (!isSymmetric())
897  return false;
898 
899  // inspect shear components
900  const Real mu = (*this)(0, 1, 0, 1);
901  // ...diagonal
902  if ((*this)(1, 2, 1, 2) != mu || (*this)(2, 0, 2, 0) != mu)
903  return false;
904  // ...off-diagonal
905  if ((*this)(2, 0, 1, 2) != 0.0 || (*this)(0, 1, 1, 2) != 0.0 || (*this)(0, 1, 2, 0) != 0.0)
906  return false;
907 
908  // off diagonal blocks in Voigt
909  unsigned int i1 = 0;
910  for (unsigned int i = 0; i < N; ++i)
911  {
912  for (unsigned int j = 0; j < N; ++j)
913  if (_vals[i1 + ((j + 1) % N) * N + (j + 2) % N] != 0.0)
914  return false;
915  i1 += N3 + N2;
916  }
917 
918  // top left block
919  const Real K1 = (*this)(0, 0, 0, 0);
920  const Real K2 = (*this)(0, 0, 1, 1);
921  if (!MooseUtils::relativeFuzzyEqual(K1 - 4.0 * mu / 3.0, K2 + 2.0 * mu / 3.0))
922  return false;
923  if ((*this)(1, 1, 1, 1) != K1 || (*this)(2, 2, 2, 2) != K1)
924  return false;
925  for (unsigned int i = 1; i < N; ++i)
926  for (unsigned int j = 0; j < i; ++j)
927  if ((*this)(i, i, j, j) != K2)
928  return false;
929 
930  return true;
931 }
RankFourTensor operator/(const Real a) const
C_ijkl/a.
bool relativeFuzzyEqual(const libMesh::Real &var1, const libMesh::Real &var2, const libMesh::Real &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
Function to check whether two variables are equal within a relative tolerance.
void fillAntisymmetricIsotropicFromInputVector(const std::vector< Real > &input)
fillAntisymmetricIsotropicFromInputVector takes 1 input to fill the the antisymmetric Rank-4 tensor w...
void dataLoad(std::istream &stream, RankFourTensor &rft, void *context)
RealVectorValue RealGradient
Definition: Assembly.h:43
RankFourTensor operator-() const
-C_ijkl
void fillGeneralFromInputVector(const std::vector< Real > &input)
fillGeneralFromInputVector takes 81 inputs to fill the Rank-4 tensor No symmetries are explicitly mai...
int eps(unsigned int i, unsigned int j)
2D version
FillMethod
To fill up the 81 entries in the 4th-order tensor, fillFromInputVector is called with one of the foll...
void fillSymmetricIsotropicEandNu(Real E, Real nu)
Real _vals[N2]
The values of the rank-two tensor stored by index=(i * LIBMESH_DIM + j)
RankTwoTensor innerProductTranspose(const RankTwoTensor &) const
Inner product of the major transposed tensor with a rank two tensor.
void zero()
Zeros out the tensor.
void mooseError(Args &&...args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:182
void fillGeneralIsotropic(Real i0, Real i1, Real i2)
Vector-less fill API functions. See docs of the corresponding ...FromInputVector methods.
void rotate(const T &R)
Rotate the tensor using C_ijkl = R_im R_in R_ko R_lp C_mnop.
RankFourTensor & operator=(const RankFourTensor &a)
copies values from a into this tensor
void surfaceFillFromInputVector(const std::vector< Real > &input)
Fills the tensor entries ignoring the last dimension (ie, C_ijkl=0 if any of i, j, k, or l = 3).
bool isIsotropic() const
checks if the tensor is isotropic
void fillSymmetricIsotropicEandNuFromInputVector(const std::vector< Real > &input)
fillSymmetricIsotropicEandNuFromInputVector is a variation of the fillSymmetricIsotropicFromInputVect...
RankFourTensor()
Default constructor; fills to zero.
RankFourTensor invSymm() const
This returns A_ijkl such that C_ijkl*A_klmn = 0.5*(de_im de_jn + de_in de_jm) This routine assumes th...
Real L2norm() const
sqrt(C_ijkl*C_ijkl)
void inverse(const std::vector< std::vector< Real >> &m, std::vector< std::vector< Real >> &m_inv)
Inverse the dense square matrix m using LAPACK routines.
Definition: MatrixTools.C:25
void fillFromInputVector(const std::vector< Real > &input, FillMethod fill_method)
fillFromInputVector takes some number of inputs to fill the Rank-4 tensor.
InitMethod
Initialization method.
RankTwoTensor operator*(const RankTwoTensor &a) const
C_ijkl*a_kl.
bool isSymmetric() const
checks if the tensor is symmetric
RealGradient sum3x1() const
Calculates the vector a[i] = sum over j Ciijj for i and j varying from 0 to 2.
RankFourTensor transposeMajor() const
Transpose the tensor by swapping the first pair with the second pair of indices.
static constexpr unsigned int N2
PetscInt m
RankFourTensor & operator*=(const Real a)
C_ijkl *= a.
void fillAntisymmetricFromInputVector(const std::vector< Real > &input)
fillAntisymmetricFromInputVector takes 6 inputs to fill the the antisymmetric Rank-4 tensor with the ...
void fillAntisymmetricIsotropic(Real i0)
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:37
RankTwoTensor is designed to handle the Stress or Strain Tensor for a fully anisotropic material...
Definition: RankTwoTensor.h:45
static constexpr unsigned int N4
Real _vals[N4]
The values of the rank-four tensor stored by index=(((i * LIBMESH_DIM + j) * LIBMESH_DIM + k) * LIBME...
RankFourTensor & operator+=(const RankFourTensor &a)
C_ijkl += a_ijkl for all i, j, k, l.
RankFourTensor & operator-=(const RankFourTensor &a)
C_ijkl -= a_ijkl.
void fillSymmetricIsotropicFromInputVector(const std::vector< Real > &input)
fillSymmetricIsotropicFromInputVector takes 2 inputs to fill the the symmetric Rank-4 tensor with the...
void fillSymmetricFromInputVector(const std::vector< Real > &input, bool all)
fillSymmetricFromInputVector takes either 21 (all=true) or 9 (all=false) inputs to fill in the Rank-4...
void fillGeneralIsotropicFromInputVector(const std::vector< Real > &input)
fillGeneralIsotropicFromInputVector takes 3 inputs to fill the Rank-4 tensor with symmetries C_ijkl =...
void fillSymmetricIsotropic(Real i0, Real i1)
Real sum3x3() const
Calculates the sum of Ciijj for i and j varying from 0 to 2.
PetscInt n
RankFourTensor operator+(const RankFourTensor &a) const
C_ijkl + a_ijkl.
void fillPrincipalFromInputVector(const std::vector< Real > &input)
fillPrincipalFromInputVector takes 9 inputs to fill a Rank-4 tensor C1111 = input0 C1122 = input1 C11...
void dataStore(std::ostream &stream, RankFourTensor &rft, void *context)
static MooseEnum fillMethodEnum()
Static method for use in validParams for getting the "fill_method".
RankFourTensor & operator/=(const Real a)
C_ijkl /= a for all i, j, k, l.
static constexpr unsigned int N3
TensorValue< Real > RealTensorValue
Definition: Assembly.h:45
void print(std::ostream &stm=Moose::out) const
Print the rank four tensor.
void fillAxisymmetricRZFromInputVector(const std::vector< Real > &input)
fillAxisymmetricRZFromInputVector takes 5 inputs to fill the axisymmetric Rank-4 tensor with the appr...
static constexpr unsigned int N
Dimensionality of rank-four tensor.
RankFourTensor is designed to handle any N-dimensional fourth order tensor, C.
void mooseSetToZero< RankFourTensor >(RankFourTensor &v)
Helper function template specialization to set an object to zero.