www.mooseframework.org
SymmAnisotropicElasticityTensor.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 /****************************************************************/
8 #include <cmath>
9 #include <ostream>
10 
12  : SymmElasticityTensor(false),
13  _dmat(9, 9),
14  _qdmat(9, 9),
15  _dt(6, 6),
16  _qdt(6, 6),
17  _r(3, 3),
18  _q(9, 9),
19  _qt(9, 9),
20  _euler_angle(3),
21  _trans_d6_to_d9(9, 6),
22  //_transpose_trans_d6_to_d9(6,9),
23  _trans_d9_to_d6(6, 9),
24  //_transpose_trans_d9_to_d6(9,6),
25  _c11(0),
26  _c12(0),
27  _c44(0)
28 {
29  // form_transformation_t_matrix();
30 }
31 
33  bool all_21)
34  : SymmElasticityTensor(false),
35  _dmat(9, 9),
36  _qdmat(9, 9),
37  _dt(6, 6),
38  _qdt(6, 6),
39  _r(3, 3),
40  _q(9, 9),
41  _qt(9, 9),
42  _euler_angle(3),
43  _trans_d6_to_d9(9, 6),
44  //_transpose_trans_d6_to_d9(6,9),
45  _trans_d9_to_d6(6, 9),
46  //_transpose_trans_d9_to_d6(9,6),
47  _c11(0),
48  _c12(0),
49  _c44(0)
50 {
51  // test the input vector length to make sure it's correct
52  if ((all_21 == true && init_list.size() != 21) || (all_21 == false && init_list.size() != 9))
53  mooseError("Please correct the number of entries in the stiffness input.");
54 
55  if (all_21 == true)
56  {
57  for (int i = 0; i < 21; i++)
58  _val[i] = init_list[i];
59  }
60  else
61  {
62  _val[0] = init_list[0]; // C1111
63  _val[1] = init_list[1]; // C1122
64  _val[2] = init_list[2]; // C1133
65  _val[6] = init_list[3]; // C2222
66  _val[7] = init_list[4]; // C2233
67  _val[11] = init_list[5]; // C3333
68  _val[15] = init_list[6]; // C2323
69  _val[18] = init_list[7]; // C1313
70  _val[20] = init_list[8]; // C1212
71  }
72 
73  // form_transformation_t_matrix();
74 }
75 
78  : SymmElasticityTensor(false)
79 {
80  *this = a;
81 }
82 
83 void
85 {
86  _euler_angle[0] = a1;
87 }
88 
89 void
91 {
92  _euler_angle[1] = a2;
93 }
94 
95 void
97 {
98  _euler_angle[2] = a3;
99 }
100 
101 Real
103 {
104  return _euler_angle[0];
105 }
106 
107 Real
109 {
110  return _euler_angle[1];
111 }
112 
113 Real
115 {
116  return _euler_angle[2];
117 }
118 
119 void
121 {
122  _c11 = c11;
123  _val[0] = _val[6] = _val[11] = _c11;
124 }
125 
126 void
128 {
129  _c12 = c12;
130  _val[1] = _val[2] = _val[7] = _c12;
131 }
132 
133 void
135 {
136  _c44 = c44;
137  _val[15] = _val[18] = _val[20] = _c44;
138 }
139 
140 void
141 SymmAnisotropicElasticityTensor::rotate(const Real a1, const Real a2, const Real a3)
142 {
143  setFirstEulerAngle(a1);
145  setThirdEulerAngle(a3);
146 
147  // pulled from calculateEntries to sub in the initialize_anisotropic_material_dt_matrix() call
148  // calculateEntries(0);
149 
150  form_r_matrix();
155 
156  unsigned count(0);
157 
158  for (int j(0); j < 6; ++j)
159  {
160  for (int i(j); i < 6; ++i)
161  {
162  _val[count++] = _dt(i, j);
163  }
164  }
165 }
166 
167 void
169 {
170  Real phi1 = _euler_angle[0] * (libMesh::pi / 180.0);
171  Real phi = _euler_angle[1] * (libMesh::pi / 180.0);
172  Real phi2 = _euler_angle[2] * (libMesh::pi / 180.0);
173 
174  Real cp1 = std::cos(phi1);
175  Real cp2 = std::cos(phi2);
176  Real cp = std::cos(phi);
177 
178  Real sp1 = std::sin(phi1);
179  Real sp2 = std::sin(phi2);
180  Real sp = std::sin(phi);
181 
182  _r(0, 0) = cp1 * cp2 - sp1 * sp2 * cp;
183  _r(0, 1) = sp1 * cp2 + cp1 * sp2 * cp;
184  _r(0, 2) = sp2 * sp;
185  _r(1, 0) = -cp1 * sp2 - sp1 * cp2 * cp;
186  _r(1, 1) = -sp1 * sp2 + cp1 * cp2 * cp;
187  _r(1, 2) = cp2 * sp;
188  _r(2, 0) = sp1 * sp;
189  _r(2, 1) = -cp1 * sp;
190  _r(2, 2) = cp;
191 }
192 
193 // this is now obsolete
194 void
196 {
197  // This function initializes the 6 x 6 material Dt matrix for a cubic material
198 
199  _dt(0, 0) = _dt(1, 1) = _dt(2, 2) = _c11;
200  _dt(0, 1) = _dt(0, 2) = _dt(1, 0) = _dt(2, 0) = _dt(1, 2) = _dt(2, 1) = _c12;
201  // beware the factor of two here
202  _dt(3, 3) = _dt(4, 4) = _dt(5, 5) = 2 * _c44;
203 }
204 
205 // this is now obsolete
206 void
208 {
209  // This function initializes the 6 x 6 material Dt matrix for an anisotropic material
210 
211  int k = 0;
212  for (int i = 0; i < 6; i++)
213  {
214  for (int j = i; j < 6; j++)
215  {
216  _dt(i, j) = _dt(j, i) = _val[k++];
217  }
218  }
219 }
220 
221 void
223 {
224 
225  for (int i = 0; i < 3; ++i)
226  for (int j = 0; j < 3; ++j)
227  for (int k = 0; k < 3; ++k)
228  for (int l = 0; l < 3; ++l)
229  _q(((i * 3) + k), ((j * 3) + l)) = _r(i, j) * _r(k, l);
230 
231  /*for (int p = 0; p < 9; ++p)
232  for (int q = 0; q < 9; ++q)
233  _qt(q,p) = _q(p,q);*/
234 }
235 
236 void
238 {
239  // Forms to two kinds of transformation matrix
240  // TransD6toD9 transfroms Dt[6][6] to Dmat[9][9]
241  // TransD9toD6 transforms Dmat[9][9] to Dt[6][6]
242 
243  // Real sqrt2 = std::sqrt(2.0);
244  // Real a = 1/sqrt2;
245  Real a = 1.0;
246 
247  _trans_d6_to_d9(0, 0) = _trans_d6_to_d9(4, 1) = _trans_d6_to_d9(8, 2) = 1.0;
248  _trans_d6_to_d9(1, 3) = _trans_d6_to_d9(3, 3) = a;
249  _trans_d6_to_d9(5, 4) = _trans_d6_to_d9(7, 4) = a;
250  _trans_d6_to_d9(2, 5) = _trans_d6_to_d9(6, 5) = a;
251 
252  /*for (int i = 0; i < 9; ++i)
253  {
254  for (int j = 0; j < 6; ++j)
255  {
256  _transpose_trans_d6_to_d9(j,i) = _trans_d6_to_d9(i,j);
257  }
258  }*/
259 
260  _trans_d9_to_d6(0, 0) = _trans_d9_to_d6(1, 4) = _trans_d9_to_d6(2, 8) = 1.0;
261  // _trans_d9_to_d6(3,3) = _trans_d9_to_d6(4,7) = _trans_d9_to_d6(5,6) = sqrt2;
262  _trans_d9_to_d6(3, 3) = _trans_d9_to_d6(4, 7) = _trans_d9_to_d6(5, 6) = 1.0;
263 
264  /*for (int i = 0; i < 6; ++i)
265  {
266  for (int j = 0; j < 9; ++j)
267  {
268  _transpose_trans_d9_to_d6(j,i) = _trans_d9_to_d6(i,j);
269  }
270  }*/
271 }
272 
273 void
275 {
276 
277  // THIS TRANSFORMATION IS VALID ONLY WHEN THE INCOMING MATRIX HAS NOT BEEN ROTATED
278 
279  // The function makes use of TransD6toD9 matrix to transfrom
280  // Dt[6][6] to Dmat[9][9]
281  // Dmat = T * Dt * TT
282 
283  // DenseMatrix<Real> outputMatrix(9,6);
284 
285  // outputMatrix = _trans_d6_to_d9;
286  // outputMatrix.right_multiply(_dt);
287 
288  //_dmat = outputMatrix;
289  //_dmat.right_multiply_transpose(_trans_d6_to_d9);
290 
291  // Use the plug-and-chug functions given in SymmElasticityTensor
292  // to take the existing _val[] and put them into the 9x9 _dmat matrix.
293  SymmElasticityTensor temp_dt;
294  copyValues(temp_dt);
295 
296  ColumnMajorMatrix temp_dmat = temp_dt.columnMajorMatrix9x9();
297 
298  for (unsigned j(0); j < 9; ++j)
299  {
300  for (unsigned i(0); i < 9; ++i)
301  {
302  _dmat(i, j) = temp_dmat(i, j);
303  }
304  }
305 }
306 
307 void
309 {
310  // The function makes use of TransD6toD9 matrix to transfrom
311  // QDmat[9][9] to Dt[6][6]
312  // Dt = TT * QDmat * T
313 
314  // DenseMatrix<Real> outputMatrix(6,9);
315 
316  // outputMatrix = _trans_d9_to_d6;
317  // outputMatrix.right_multiply(_qdmat);
318  // _dt = outputMatrix;
319  // _dt.right_multiply(_transpose_trans_d9_to_d6);
320 
321  // The transformation below is general and should work whether or not the
322  // incoming tensor has been rotated.
323 
324  ColumnMajorMatrix tmp(9, 9);
325  for (unsigned j(0); j < 9; ++j)
326  {
327  for (unsigned i(0); i < 9; ++i)
328  {
329  tmp(i, j) = _qdmat(i, j);
330  }
331  }
332 
334  fred.convertFrom9x9(tmp);
335  ColumnMajorMatrix wilma = fred.columnMajorMatrix6x6();
336  for (unsigned j(0); j < 6; ++j)
337  {
338  for (unsigned i(0); i < 6; ++i)
339  {
340  _dt(i, j) = wilma(i, j);
341  }
342  }
343 }
344 
345 void
347 {
348  // The function makes use of Q matrix to rotate
349  // Dmat[9][9] to QDmat[9][9]
350  // QDmat = QT * Dmat * Q
351 
352  DenseMatrix<Real> outputMatrix(9, 9);
353 
354  _q.get_transpose(outputMatrix);
355  outputMatrix.right_multiply(_dmat);
356 
357  _qdmat = outputMatrix;
358  _qdmat.right_multiply(_q);
359 }
360 
361 void
363 {
364 
365  // The following four lines of code force the calculateEntries function to be useful
366  // only for CUBIC ANISOTROPIC materials.
367  zero();
371 
372  form_r_matrix();
373  // initialize_material_anisotropic_dt_matrix();
375  // form_transformation_t_matrix();
379 
380  unsigned count(0);
381 
382  for (int j(0); j < 6; ++j)
383  {
384  for (int i(j); i < 6; ++i)
385  {
386  _val[count++] = _dt(i, j);
387  }
388  }
389 }
390 
391 void
393 {
394  printf("\nSymmAnisotropicElasticityTensor::show_dt_matrix()\n");
395 
396  for (int j = 0; j < 6; ++j)
397  {
398  printf(" ");
399  for (int i = 0; i < 6; ++i)
400  {
401  printf("%12.4f ", _dt(i, j));
402  }
403  printf("\n");
404  }
405 }
406 
407 void
409 {
410  printf("\nSymmAnisotropicElasticityTensor::show_r_matrix() Euler angles are (%f, %f, %f)\n",
411  _euler_angle[0],
412  _euler_angle[1],
413  _euler_angle[2]);
414 
415  for (int j = 0; j < 3; ++j)
416  {
417  printf(" ");
418  for (int i = 0; i < 3; ++i)
419  {
420  printf("%8.4f ", _r(i, j));
421  }
422  printf("\n");
423  }
424 }
virtual void rotate(const Real a1, const Real a2, const Real a3)
Perform rotation around three axes.
This class defines a basic set of capabilities any elasticity tensor should have. ...
virtual void calculateEntries(unsigned int qp)
Fill in the matrix.
void setMaterialConstantc11(const Real c11)
Set the material constant c11; assumes cubic material.
void setMaterialConstantc44(const Real c44)
Set the material constant c44; assumes cubic material.
void setThirdEulerAngle(const Real a3)
Set the third euler angle.
ColumnMajorMatrix columnMajorMatrix6x6() const
void convertFrom9x9(const ColumnMajorMatrix &cmm)
void setSecondEulerAngle(const Real a2)
Set the second euler angle.
void setFirstEulerAngle(const Real a1)
Set the first euler angle.
ColumnMajorMatrix columnMajorMatrix9x9() const
void setMaterialConstantc12(const Real c12)
Set the material constant c22; assumes cubic material.
void copyValues(SymmElasticityTensor &rhs) const