AnisotropicElasticityTensor.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
10  : ElasticityTensor(false), _euler_angle(3), _c11(0), _c12(0), _c44(0)
11 {
12 }
13
14 void
16 {
17  _euler_angle[0] = a1;
18 }
19
20 void
22 {
23  _euler_angle[1] = a2;
24 }
25
26 void
28 {
29  _euler_angle[2] = a3;
30 }
31
32 void
34 {
35  _c11 = c11;
36 }
37
38 void
40 {
41  _c12 = c12;
42 }
43
44 void
46 {
47  _c44 = c44;
48 }
49
50 void
52 {
53
54  // form_r_matrix();
55  // Form rotation matrix from Euler angles
56  const Real phi1 = _euler_angle[0] * (M_PI / 180.0);
57  const Real phi = _euler_angle[1] * (M_PI / 180.0);
58  const Real phi2 = _euler_angle[2] * (M_PI / 180.0);
59
60  Real cp1 = cos(phi1);
61  Real cp2 = cos(phi2);
62  Real cp = cos(phi);
63
64  Real sp1 = sin(phi1);
65  Real sp2 = sin(phi2);
66  Real sp = sin(phi);
67
68  DenseMatrix<Real> R; // Rotational Matrix
69  R(0, 0) = cp1 * cp2 - sp1 * sp2 * cp;
70  R(0, 1) = sp1 * cp2 + cp1 * sp2 * cp;
71  R(0, 2) = sp2 * sp;
72  R(1, 0) = -cp1 * sp2 - sp1 * cp2 * cp;
73  R(1, 1) = -sp1 * sp2 + cp1 * cp2 * cp;
74  R(1, 2) = cp2 * sp;
75  R(2, 0) = sp1 * sp;
76  R(2, 1) = -cp1 * sp;
77  R(2, 2) = cp;
78
79  // Initialize material Dt matrix
80  DenseMatrix<Real> Dt; // 6 x 6 Material Matrix
81  Dt(0, 0) = Dt(1, 1) = Dt(2, 2) = _c11;
82  Dt(0, 1) = Dt(0, 2) = Dt(1, 0) = Dt(2, 0) = Dt(1, 2) = Dt(2, 1) = _c12;
83  Dt(3, 3) = Dt(4, 4) = Dt(5, 5) = 2 * _c44;
84
85  // Form Q = R dyadic R and Q transpose
86  DenseMatrix<Real> Q, Qt; // Q = R (dyadic) R and Q transpose
87  for (unsigned int i = 0; i < 3; i++)
88  for (unsigned int j = 0; j < 3; j++)
89  for (unsigned int k = 0; k < 3; k++)
90  for (unsigned int l = 0; l < 3; l++)
91  Q(((i * 3) + k), ((j * 3) + l)) = R(i, j) * R(k, l);
92
93  for (unsigned int p = 0; p < 9; p++)
94  for (unsigned int q = 0; q < 9; q++)
95  Qt(q, p) = Q(p, q);
96
97  // Form two kinds of transformation matrix
98  // TransD6toD9 transfroms Dt[6][6] to Dmat[9][9]
99  // TransD9toD6 transforms Dmat[9][9] to Dt[6][6]
100  DenseMatrix<Real> trans_d6_to_d9, transpose_trans_d6_to_d9;
101  DenseMatrix<Real> trans_d9_to_d6, transpose_trans_d9_to_d6;
102  Real sqrt2 = std::sqrt(2.0);
103  Real a = 1 / sqrt2;
104
105  trans_d6_to_d9(0, 0) = trans_d6_to_d9(4, 1) = trans_d6_to_d9(8, 2) = 1.0;
106  trans_d6_to_d9(1, 3) = trans_d6_to_d9(3, 3) = a;
107  trans_d6_to_d9(2, 4) = trans_d6_to_d9(6, 4) = a;
108  trans_d6_to_d9(5, 5) = trans_d6_to_d9(7, 5) = a;
109
110  for (unsigned int i = 0; i < 9; i++)
111  for (unsigned int j = 0; j < 6; j++)
112  transpose_trans_d6_to_d9(j, i) = trans_d6_to_d9(i, j);
113
114  trans_d9_to_d6(0, 0) = trans_d9_to_d6(1, 4) = trans_d9_to_d6(2, 8) = 1.0;
115  trans_d9_to_d6(3, 3) = trans_d9_to_d6(4, 6) = trans_d9_to_d6(5, 7) = sqrt2;
116
117  for (unsigned int i = 0; i < 6; i++)
118  for (unsigned int j = 0; j < 9; j++)
119  transpose_trans_d9_to_d6(j, i) = trans_d9_to_d6(i, j);
120
121  // The function makes use of TransD6toD9 matrix to transfrom Dt[6][6] to Dmat[9][9]
122  // Dmat = T * Dt * TT
123
124  DenseMatrix<Real> outputMatrix(9, 6);
125
126  outputMatrix = trans_d6_to_d9;
127  outputMatrix.right_multiply(Dt);
128
129  DenseMatrix<Real> Dmat; // 9 x 9 Material Matrix
130  Dmat = outputMatrix;
131  Dmat.right_multiply(transpose_trans_d6_to_d9);
132
133  // The function makes use of Q matrix to rotate Dmat[9][9] to QDmat[9][9]
134  // QDmat = QT * Dmat * Q
135
136  DenseMatrix<Real> outputMatrix99(9, 9);
137
138  outputMatrix99 = Qt;
139  outputMatrix99.right_multiply(Dmat);
140
141  DenseMatrix<Real> QDmat; // Rotated Material Matrix
142  QDmat = outputMatrix99;
143  QDmat.right_multiply(Q);
144
145  // Convert 9X9 matrix QDmat to 81 vector
146
147  unsigned int count = 0;
148
149  for (unsigned int j = 0; j < 9; j++)
150  for (unsigned int i = 0; i < 9; i++)
151  {
152  _values[count] = QDmat(i, j);
153  count = count + 1;
154  }
155 }
This class defines a basic set of capabilities any elasticity tensor should have. ...
virtual void calculateEntries(unsigned int qp)
Fill in the matrix.
void setFirstEulerAngle(const Real a1)
Set the first euler angle.
void setMaterialConstantc11(const Real c11)
Set the material constant c11.
void setMaterialConstantc12(const Real c12)
Set the material constant c22.
void setThirdEulerAngle(const Real a3)
Set the third euler angle.
void setMaterialConstantc44(const Real c44)
Set the material constant c44.
void setSecondEulerAngle(const Real a2)
Set the second euler angle.