Alamo
Matrix4_MajorMinor.H
Go to the documentation of this file.
1#ifndef SET_MATRIX4_MAJORMINOR_H
2#define SET_MATRIX4_MAJORMINOR_H
3
4#include "Base.H"
5#include <math.h>
6
7namespace Set
8{
9/// \brief Data structure for a 4th order 3D tensor with major and minor symmetry
10///
11/// 2D version. See full explanation below.
12///
13template<>
15{
16 Scalar data[6] = {NAN,NAN,NAN,NAN,NAN,NAN};
17
18public:
20 //Matrix4(const Matrix4<2,Sym::MajorMinor> &in)
21 //{
22 // for (int i = 0; i < 6; i++) data[i] = in.data[i];
23 //}
25 const Scalar & operator () (const int i, const int j, const int k, const int l) const
26 {
27 int uid = i + 2*j + 4*k + 8*l;
28 if (uid==0 ) return data[0]; // [0000]
29 else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1]; // [0001] [0010] [0100] [1000]
30 else if (uid==12 || uid==3 ) return data[2]; // [0011] [1100]
31 else if (uid==10 || uid==6 || uid==9 || uid==5 ) return data[3]; // [0101] [1001] [0110] [1010]
32 else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[4]; // [0111] [1011] [1101] [1110]
33 else if (uid==15 ) return data[5]; // [1111]
34 else Util::Abort(INFO,"Index out of range");
35 return Set::Garbage;
36 }
38 Scalar & operator () (const int i, const int j, const int k, const int l)
39 {
40 int uid = i + 2*j + 4*k + 8*l;
41 if (uid==0 ) return data[0]; // [0000]
42 else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1]; // [0001] [0010] [0100] [1000]
43 else if (uid==12 || uid==3 ) return data[2]; // [0011] [1100]
44 else if (uid==10 || uid==6 || uid==9 || uid==5 ) return data[3]; // [0101] [1001] [0110] [1010]
45 else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[4]; // [0111] [1011] [1101] [1110]
46 else if (uid==15 ) return data[5]; // [1111]
47 else Util::Abort(INFO,"Index out of range");
48 return Set::Garbage;
49 }
50 void Print (std::ostream& os)
51 {
52 os.precision(4);
53 for (int k = 0; k < 2; k++)
54 {
55 for (int i = -1; i < 3; i++)
56 {
57 for (int l = 0; l < 2; l++)
58 {
59 if (i==-1) os << "┌ ┐ ";
60 else if (i == 2) os << "└ ┘ ";
61 else
62 {
63 os << "│ ";
64 for (int j = 0; j < 2; j++)
65 {
66 const Set::Scalar &val = (*this)(i,j,k,l);
67 os << std::scientific << std::setw(11) << val ; //(fabs(val)>1E-10 ? val : 0);
68 os << " ";
69 }
70 os << "│ ";
71 }
72 }
73 os<<std::endl;
74 }
75 }
76 }
77
79 {
80 return std::sqrt(data[0]*data[0] + data[1]*data[1] + data[2]*data[2] + data[3]*data[3] + data[4]*data[4] + data[5]*data[5]);
81 }
82 bool contains_nan() const
83 {
84 if (std::isnan(data[ 0])) return true;
85 if (std::isnan(data[ 1])) return true;
86 if (std::isnan(data[ 2])) return true;
87 if (std::isnan(data[ 3])) return true;
88 if (std::isnan(data[ 4])) return true;
89 if (std::isnan(data[ 5])) return true;
90 return false;
91 }
92
93 AMREX_GPU_HOST_DEVICE void operator += (Matrix4<2,Sym::MajorMinor> a) {for (int i = 0; i < 6; i++) data[i] += a.data[i];}
94 AMREX_GPU_HOST_DEVICE void operator -= (Matrix4<2,Sym::MajorMinor> a) {for (int i = 0; i < 6; i++) data[i] -= a.data[i];}
95 AMREX_GPU_HOST_DEVICE void operator *= (Matrix4<2,Sym::MajorMinor> a) {for (int i = 0; i < 6; i++) data[i] *= a.data[i];}
96 AMREX_GPU_HOST_DEVICE void operator /= (Matrix4<2,Sym::MajorMinor> a) {for (int i = 0; i < 6; i++) data[i] /= a.data[i];}
97 AMREX_GPU_HOST_DEVICE void operator *= (Set::Scalar alpha) {for (int i = 0; i < 6; i++) data[i] *= alpha;}
98 AMREX_GPU_HOST_DEVICE void operator /= (Set::Scalar alpha) {for (int i = 0; i < 6; i++) data[i] /= alpha;}
99
101 {
103 for (int i = 0 ; i < 6; i++) ret.data[i] = (Set::Scalar)i;
104 return ret;
105 }
107 {
109 for (int i = 0 ; i < 6; i++) ret.data[i] = (Util::Random());
110 return ret;
111 }
113 {
115 for (int i = 0 ; i < 6; i++) ret.data[i] = 0.0;
116 return ret;
117 }
119 Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
120 {
122 Set::Scalar Ctmp[3][3][3][3];
123
124 for(int i = 0; i < 3; i++)
125 for(int j = 0; j < 3; j++)
126 for(int k = 0; k < 3; k++)
127 for(int l = 0; l < 3; l++)
128 {
129 if(i == j && j == k && k == l) Ctmp[i][j][k][l] = C11;
130 else if (i==k && j==l) Ctmp[i][j][k][l] = C44;
131 else if (i==j && k==l) Ctmp[i][j][k][l] = C12;
132 else Ctmp[i][j][k][l] = 0.0;
133 }
134 for(int p = 0; p < 2; p++)
135 for(int q = 0; q < 2; q++)
136 for(int s = 0; s < 2; s++)
137 for(int t = 0; t < 2; t++)
138 {
139 ret(p,q,s,t) = 0.0;
140 for(int i = 0; i < 3; i++)
141 for(int j = 0; j < 3; j++)
142 for(int k = 0; k < 3; k++)
143 for(int l = 0; l < 3; l++)
144 ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
145 }
146 return ret;
147 }
150 {
151 Eigen::Matrix3d R;
152 R = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
153 Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
154 Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
155 return Matrix4<2,Sym::MajorMinor>::Cubic(C11,C12,C44,R);
156 }
159 {
160 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
161 return Cubic(C11,C12,C44,R);
162 }
163
164
165 static Matrix4<2,Sym::MajorMinor> Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
166 {
168 Set::Scalar Ctmp[3][3][3][3] = {{{{0.0}}}};
169
170 for(int i = 0; i < 3; i++)
171 for(int j = 0; j < 3; j++)
172 for(int k = 0; k < 3; k++)
173 for(int l = 0; l < 3; l++)
174 {
175 if(i==j && j==k && k==l)
176 {
177 if (i < 2) Ctmp[i][j][k][l] = C11;
178 else Ctmp[i][j][k][l] = C33;
179 }
180 else if(i==j && k==l)
181 {
182 if (i<2 && k<2 && i!=k) Ctmp[i][j][k][l] = C12;
183 else if ((i<2 && k==2) || (i==2 && k<2)) Ctmp[i][j][k][l] = C13;
184 }
185 else if((i==k && j==l) && (i!=j))
186 {
187 Ctmp[i][j][k][l] = C44;
188 }
189 else
190 {
191 Ctmp[i][j][k][l] = 0.0;
192 }
193 }
194
195 // C66 = 0.5*(C11 - C12)
196 Ctmp[0][1][0][1] = 0.5 * (C11 - C12);
197 Ctmp[1][0][1][0] = 0.5 * (C11 - C12);
198
199 for(int p = 0; p < 2; p++)
200 for(int q = 0; q < 2; q++)
201 for(int s = 0; s < 2; s++)
202 for(int t = 0; t < 2; t++)
203 {
204 ret(p,q,s,t) = 0.0;
205 for(int i = 0; i < 3; i++)
206 for(int j = 0; j < 3; j++)
207 for(int k = 0; k < 3; k++)
208 for(int l = 0; l < 3; l++)
209 ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
210 }
211 return ret;
212 }
214 {
215 Eigen::Quaterniond q =
216 Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
217 Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
218 Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
219 Eigen::Matrix3d R = q.toRotationMatrix();
220 return Matrix4<2,Sym::MajorMinor>::Transverse(C11,C12,C13,C33,C44,R);
221 }
223 {
224 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
225 return Transverse(C11,C12,C13,C33,C44,R);
226 }
227
228
229 friend Eigen::Matrix<Set::Scalar,2,2> operator * (Matrix4<2,Sym::MajorMinor> a, Eigen::Matrix<Set::Scalar,2,2> b);
237};
240{
241 for (int i = 0 ; i < 6 ; i++) if (a.data[i] != b.data[i]) return false;
242 return true;
243}
246{
248 for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] + b.data[i];
249 return ret;
250}
253{
255 for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] - b.data[i];
256 return ret;
257}
260{
262 for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] * b;
263 return ret;
264}
272{
274 for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] / b;
275 return ret;
276}
278Eigen::Matrix<Set::Scalar,2,2> operator * (Matrix4<2,Sym::MajorMinor> a, Eigen::Matrix<Set::Scalar,2,2> b)
279{
280 Eigen::Matrix<Set::Scalar,2,2> ret = Eigen::Matrix<Set::Scalar,2,2>::Zero();
281 ret(0,0) = a.data[0]*b(0,0) + a.data[1]*(b(0,1) + b(1,0)) + a.data[2]*b(1,1);
282 ret(0,1) = a.data[1]*b(0,0) + a.data[3]*(b(0,1) + b(1,0)) + a.data[4]*b(1,1);
283 ret(1,1) = a.data[2]*b(0,0) + a.data[4]*(b(0,1) + b(1,0)) + a.data[5]*b(1,1);
284 ret(1,0) = ret(0,1);
285 return ret;
286}
288Eigen::Matrix<amrex::Real,2,1> operator * (Matrix4<2,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,2,2>,2> b);
289
292{
293 Set::Vector ret = Set::Vector::Zero();
294 ret(0) = a.data[0]*b(0,0,0) + a.data[1]*(b(0,1,0) + b(1,0,0) + b(0,0,1)) + a.data[2]*b(1,1,0) + a.data[3]*(b(0,1,1) + b(1,0,1)) + a.data[4]*b(1,1,1);
295 ret(1) = a.data[1]*b(0,0,0) + a.data[2]*b(0,0,1) + a.data[3]*(b(0,1,0) + b(1,0,0)) + a.data[4]*(b(1,1,0) + b(0,1,1) + b(1,0,1)) + a.data[5]*b(1,1,1);
296 return ret;
297}
298
299/// \brief Data structure for a 4th order 3D tensor with major and minor symmetry
300///
301/// Let the tensor
302/// \f$\mathbb{C}\in\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{R}^3\f$
303/// be have major symmetry (\f$\mathbb{C}_{ijkl}=\mathbb{C}_{klij}\f$)
304/// and minor symmetry (\f$\mathbb{C}_{ijkl}=\mathbb{C}_{jikl}=\mathbb{C}_{ijlk}\f$)
305/// Then there are only 21 unique elements (rather than 81).
306///
307/// This object acts like a 4D array such that `C(i,j,k,l)` returns the corresponding
308/// element, but symmetry is always obeyed. This allows the user code to be much prettier,
309/// while maintaining a relatively small object size.
310///
311template<>
313{
314 Scalar data[21] = { NAN,NAN,NAN,NAN,NAN,NAN,NAN,
316 NAN,NAN,NAN,NAN,NAN,NAN,NAN };
317
318public:
321 const Scalar & operator () (const int i, const int j, const int k, const int l) const
322 {
323 int uid = i + 3*j + 9*k + 27*l;
324 // [0000]
325 if (uid==0 ) return data[0];
326 // [0001] [0010] [0100] [1000]
327 else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
328 // [0002] [0020] [0200] [2000]
329 else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
330 // [0011] [1100]
331 else if (uid==36 || uid==4 ) return data[3];
332 // [0012] [0021] [1200] [2100]
333 else if (uid==63 || uid==45 || uid==7 || uid==5 ) return data[4];
334 // [0022] [2200]
335 else if (uid==72 || uid==8 ) return data[5];
336 // [0101] [0110] [1001] [1010]
337 else if (uid==30 || uid==12 || uid==28 || uid==10 ) return data[6];
338 // [0102] [0120] [1002] [1020] [0201] [2001] [0210] [2010]
339 else if (uid==57 || uid==21 || uid==33 || uid==15 || uid==55 || uid==19 || uid==29 || uid==11 ) return data[7];
340 // [0111] [1011] [1101] [1110]
341 else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[8];
342 // [0112] [1012] [0121] [1021] [1201] [1210] [2101] [2110]
343 else if (uid==66 || uid==48 || uid==64 || uid==46 || uid==34 || uid==16 || uid==32 || uid==14 ) return data[9];
344 // [0122] [1022] [2201] [2210]
345 else if (uid==75 || uid==73 || uid==35 || uid==17 ) return data[10];
346 // [0202] [2002] [0220] [2020]
347 else if (uid==60 || uid==24 || uid==56 || uid==20 ) return data[11];
348 // [0211] [2011] [1102] [1120]
349 else if (uid==42 || uid==58 || uid==22 || uid==38 ) return data[12];
350 // [0212] [2012] [0221] [2021] [1202] [1220] [2102] [2120]
351 else if (uid==69 || uid==51 || uid==61 || uid==25 || uid==65 || uid==47 || uid==59 || uid==23 ) return data[13];
352 // [0222] [2022] [2202] [2220]
353 else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[14];
354 // [1111]
355 else if (uid==40 ) return data[15];
356 // [1112] [1121] [1211] [2111]
357 else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[16];
358 // [1122] [2211]
359 else if (uid==76 || uid==44 ) return data[17];
360 // [1212] [2112] [1221] [2121]
361 else if (uid==70 || uid==52 || uid==68 || uid==50 ) return data[18];
362 // [1222] [2122] [2212] [2221]
363 else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[19];
364 // [2222]
365 else if (uid==80 ) return data[20];
366 else Util::Abort(INFO,"Index out of range");
367 return Set::Garbage;
368 }
370 Scalar & operator () (const int i, const int j, const int k, const int l)
371 {
372 int uid = i + 3*j + 9*k + 27*l;
373 // [0000]
374 if (uid==0 ) return data[0];
375 // [0001] [0010] [0100] [1000]
376 else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
377 // [0002] [0020] [0200] [2000]
378 else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
379 // [0011] [1100]
380 else if (uid==36 || uid==4 ) return data[3];
381 // [0012] [0021] [1200] [2100]
382 else if (uid==63 || uid==45 || uid==7 || uid==5 ) return data[4];
383 // [0022] [2200]
384 else if (uid==72 || uid==8 ) return data[5];
385 // [0101] [0110] [1001] [1010]
386 else if (uid==30 || uid==12 || uid==28 || uid==10 ) return data[6];
387 // [0102] [0120] [1002] [1020] [0201] [2001] [0210] [2010]
388 else if (uid==57 || uid==21 || uid==33 || uid==15 || uid==55 || uid==19 || uid==29 || uid==11 ) return data[7];
389 // [0111] [1011] [1101] [1110]
390 else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[8];
391 // [0112] [1012] [0121] [1021] [1201] [1210] [2101] [2110]
392 else if (uid==66 || uid==48 || uid==64 || uid==46 || uid==34 || uid==16 || uid==32 || uid==14 ) return data[9];
393 // [0122] [1022] [2201] [2210]
394 else if (uid==75 || uid==73 || uid==35 || uid==17 ) return data[10];
395 // [0202] [2002] [0220] [2020]
396 else if (uid==60 || uid==24 || uid==56 || uid==20 ) return data[11];
397 // [0211] [2011] [1102] [1120]
398 else if (uid==42 || uid==58 || uid==22 || uid==38 ) return data[12];
399 // [0212] [2012] [0221] [2021] [1202] [1220] [2102] [2120]
400 else if (uid==69 || uid==51 || uid==61 || uid==25 || uid==65 || uid==47 || uid==59 || uid==23 ) return data[13];
401 // [0222] [2022] [2202] [2220]
402 else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[14];
403 // [1111]
404 else if (uid==40 ) return data[15];
405 // [1112] [1121] [1211] [2111]
406 else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[16];
407 // [1122] [2211]
408 else if (uid==76 || uid==44 ) return data[17];
409 // [1212] [2112] [1221] [2121]
410 else if (uid==70 || uid==52 || uid==68 || uid==50 ) return data[18];
411 // [1222] [2122] [2212] [2221]
412 else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[19];
413 // [2222]
414 else if (uid==80 ) return data[20];
415 else Util::Abort(INFO,"Index out of range");
416 return Set::Garbage;
417 }
418 void Print (std::ostream& os)
419 {
420 for (int i = 0; i < 21; i++)
421 os << "i = " << i << " " << data[i] << std::endl;
422
423 os.precision(4);
424 for (int k = 0; k < 3; k++)
425 {
426 for (int i = -1; i < 4; i++)
427 {
428 for (int l = 0; l < 3; l++)
429 {
430 if (i==-1) os << "┌ ┐ ";
431 else if (i == 3) os << "└ ┘ ";
432 else
433 {
434 os << "│ ";
435 for (int j = 0; j < 3; j++)
436 {
437 const Set::Scalar &val = (*this)(i,j,k,l);
438 os << std::scientific << std::setw(11) << val ; //(fabs(val)>1E-10 ? val : 0);
439 os << " ";
440 }
441 os << "│ ";
442 }
443 }
444 os<<std::endl;
445 }
446 }
447 }
448
450 {
451 Set::Scalar retsq = 0;
452 for (int i = 0; i < 21; i++) retsq += data[i];
453 return std::sqrt(retsq);
454 }
455 bool contains_nan() const
456 {
457 if (std::isnan(data[ 0])) return true;
458 if (std::isnan(data[ 1])) return true;
459 if (std::isnan(data[ 2])) return true;
460 if (std::isnan(data[ 3])) return true;
461 if (std::isnan(data[ 4])) return true;
462 if (std::isnan(data[ 5])) return true;
463 if (std::isnan(data[ 6])) return true;
464 if (std::isnan(data[ 7])) return true;
465 if (std::isnan(data[ 8])) return true;
466 if (std::isnan(data[ 9])) return true;
467 if (std::isnan(data[10])) return true;
468 if (std::isnan(data[11])) return true;
469 if (std::isnan(data[12])) return true;
470 if (std::isnan(data[13])) return true;
471 if (std::isnan(data[14])) return true;
472 if (std::isnan(data[15])) return true;
473 if (std::isnan(data[16])) return true;
474 if (std::isnan(data[17])) return true;
475 if (std::isnan(data[18])) return true;
476 if (std::isnan(data[19])) return true;
477 if (std::isnan(data[20])) return true;
478 return false;
479 }
480
481 //AMREX_GPU_HOST_DEVICE void operator = (Matrix4<3,Sym::MajorMinor> a) {for (int i = 0; i < 21; i++) data[i] = a.data[i];}
482 AMREX_GPU_HOST_DEVICE void operator += (Matrix4<3,Sym::MajorMinor> a) {for (int i = 0; i < 21; i++) data[i] += a.data[i];}
483 AMREX_GPU_HOST_DEVICE void operator -= (Matrix4<3,Sym::MajorMinor> a) {for (int i = 0; i < 21; i++) data[i] -= a.data[i];}
484 AMREX_GPU_HOST_DEVICE void operator *= (Matrix4<3,Sym::MajorMinor> a) {for (int i = 0; i < 21; i++) data[i] *= a.data[i];}
485 AMREX_GPU_HOST_DEVICE void operator /= (Matrix4<3,Sym::MajorMinor> a) {for (int i = 0; i < 21; i++) data[i] /= a.data[i];}
486 AMREX_GPU_HOST_DEVICE void operator *= (Set::Scalar alpha) {for (int i = 0; i < 21; i++) data[i] *= alpha;}
487 AMREX_GPU_HOST_DEVICE void operator /= (Set::Scalar alpha) {for (int i = 0; i < 21; i++) data[i] /= alpha;}
488
490 {
492 for (int i = 0 ; i < 21; i++) ret.data[i] = (Set::Scalar)i;
493 return ret;
494 }
496 {
498 for (int i = 0 ; i < 21; i++) ret.data[i] = (Util::Random());
499 return ret;
500 }
502 {
504 for (int i = 0 ; i < 21; i++) ret.data[i] = 0.0;
505 return ret;
506 }
508 Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
509 {
511 Set::Scalar Ctmp[3][3][3][3];
512
513 for(int i = 0; i < 3; i++)
514 for(int j = 0; j < 3; j++)
515 for(int k = 0; k < 3; k++)
516 for(int l = 0; l < 3; l++)
517 {
518 if(i == j && j == k && k == l) Ctmp[i][j][k][l] = C11;
519 else if (i==k && j==l) Ctmp[i][j][k][l] = C44;
520 else if (i==j && k==l) Ctmp[i][j][k][l] = C12;
521 else Ctmp[i][j][k][l] = 0.0;
522 }
523 for(int p = 0; p < 3; p++)
524 for(int q = 0; q < 3; q++)
525 for(int s = 0; s < 3; s++)
526 for(int t = 0; t < 3; t++)
527 {
528 ret(p,q,s,t) = 0.0;
529 for(int i = 0; i < 3; i++)
530 for(int j = 0; j < 3; j++)
531 for(int k = 0; k < 3; k++)
532 for(int l = 0; l < 3; l++)
533 ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
534 }
535 return ret;
536 }
539 {
540 Eigen::Matrix3d R;
541 R = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
542 Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
543 Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
544 return Matrix4<3,Sym::MajorMinor>::Cubic(C11,C12,C44,R);
545 }
548 {
549 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
550 return Cubic(C11,C12,C44,R);
551 }
552
553 static Matrix4<3,Sym::MajorMinor> Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
554 {
556 Set::Scalar Ctmp[3][3][3][3] = {{{{0.0}}}};
557
558 for(int i = 0; i < 3; i++)
559 for(int j = 0; j < 3; j++)
560 for(int k = 0; k < 3; k++)
561 for(int l = 0; l < 3; l++)
562 {
563 if(i==j && j==k && k==l)
564 {
565 if (i < 2) Ctmp[i][j][k][l] = C11;
566 else Ctmp[i][j][k][l] = C33;
567 }
568 else if(i==j && k==l)
569 {
570 if (i<2 && k<2 && i!=k) Ctmp[i][j][k][l] = C12;
571 else if ((i<2 && k==2) || (i==2 && k<2)) Ctmp[i][j][k][l] = C13;
572 }
573 else if((i==k && j==l) && (i!=j))
574 {
575 Ctmp[i][j][k][l] = C44;
576 }
577 else
578 {
579 Ctmp[i][j][k][l] = 0.0;
580 }
581 }
582
583 // C66 = 0.5*(C11 - C12)
584 Ctmp[0][1][0][1] = 0.5 * (C11 - C12);
585 Ctmp[1][0][1][0] = 0.5 * (C11 - C12);
586
587 for(int p = 0; p < 3; p++)
588 for(int q = 0; q < 3; q++)
589 for(int s = 0; s < 3; s++)
590 for(int t = 0; t < 3; t++)
591 {
592 ret(p,q,s,t) = 0.0;
593 for(int i = 0; i < 3; i++)
594 for(int j = 0; j < 3; j++)
595 for(int k = 0; k < 3; k++)
596 for(int l = 0; l < 3; l++)
597 ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
598 }
599 return ret;
600 }
602 {
603 Eigen::Quaterniond q =
604 Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
605 Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
606 Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
607 Eigen::Matrix3d R = q.toRotationMatrix();
608 return Matrix4<3,Sym::MajorMinor>::Transverse(C11,C12,C13,C33,C44,R);
609 }
611 {
612 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
613 return Transverse(C11,C12,C13,C33,C44,R);
614 }
615
623 friend Eigen::Matrix<amrex::Real,3,3> operator * (Matrix4<3,Sym::MajorMinor> a, Eigen::Matrix<amrex::Real,3,3> b);
624};
627{
628 for (int i = 0 ; i < 21 ; i++) if (a.data[i] != b.data[i]) return false;
629 return true;
630}
633{
635 for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] + b.data[i];
636 return ret;
637}
640{
642 for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] - b.data[i];
643 return ret;
644}
647{
649 for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] * b;
650 return ret;
651}
659{
661 for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] / b;
662 return ret;
663}
665Eigen::Matrix<amrex::Real,3,3> operator * (Matrix4<3,Sym::MajorMinor> a, Eigen::Matrix<amrex::Real,3,3> b)
666{
667 Eigen::Matrix<amrex::Real,3,3> ret = Eigen::Matrix<amrex::Real,3,3>::Zero();
668 ret(0,0) += a.data[ 0]*b(0,0) + a.data[ 1]*(b(0,1) + b(1,0)) + a.data[ 2]*(b(0,2) + b(2,0)) + a.data[ 3]*b(1,1) + a.data[ 4]*(b(1,2) + b(2,1)) + a.data[ 5]*b(2,2);
669 ret(1,1) += a.data[ 3]*b(0,0) + a.data[ 8]*(b(0,1) + b(1,0)) + a.data[ 12]*(b(0,2) + b(2,0)) + a.data[ 15]*b(1,1) + a.data[ 16]*(b(1,2) + b(2,1)) + a.data[ 17]*b(2,2);
670 ret(2,2) += a.data[ 5]*b(0,0) + a.data[ 10]*(b(0,1) + b(1,0)) + a.data[ 14]*(b(0,2) + b(2,0)) + a.data[ 17]*b(1,1) + a.data[ 19]*(b(1,2) + b(2,1)) + a.data[ 20]*b(2,2);
671 ret(1,2) += a.data[ 4]*b(0,0) + a.data[ 9]*(b(0,1) + b(1,0)) + a.data[ 13]*(b(0,2) + b(2,0)) + a.data[ 16]*b(1,1) + a.data[ 18]*(b(1,2) + b(2,1)) + a.data[ 19]*b(2,2);
672 ret(2,0) += a.data[ 2]*b(0,0) + a.data[ 7]*(b(0,1) + b(1,0)) + a.data[ 11]*(b(0,2) + b(2,0)) + a.data[ 12]*b(1,1) + a.data[ 13]*(b(1,2) + b(2,1)) + a.data[ 14]*b(2,2);
673 ret(0,1) += a.data[ 1]*b(0,0) + a.data[ 6]*(b(0,1) + b(1,0)) + a.data[ 7]*(b(0,2) + b(2,0)) + a.data[ 8]*b(1,1) + a.data[ 9]*(b(1,2) + b(2,1)) + a.data[ 10]*b(2,2);
674
675 ret(2,1) = ret(1,2);
676 ret(0,2) = ret(2,0);
677 ret(1,0) = ret(0,1);
678
679 return ret;
680}
682Eigen::Matrix<amrex::Real,2,2> operator * (Matrix4<3,Sym::MajorMinor> a, Eigen::Matrix<amrex::Real,2,2> b)
683{
684 Eigen::Matrix<amrex::Real,2,2> ret = Eigen::Matrix<amrex::Real,2,2>::Zero();
685 for (int i = 0; i < 2; i++)
686 for (int j = 0; j < 2; j++)
687 for (int k = 0; k < 2; k++)
688 for (int l = 0; l < 2; l++)
689 ret(i,j) += a(i,j,k,l)*b(k,l);
690 return ret;
691}
693Eigen::Matrix<amrex::Real,3,1> operator * (Matrix4<3,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,3,3>,3> b);
694
696Eigen::Matrix<amrex::Real,2,1> operator * (Matrix4<3,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,2,2>,2> b);
697
700{
701 // TODO: improve efficiency of this method
702 Set::Vector ret = Set::Vector::Zero();
703
704 ret(0) +=
705 a.data[ 0]*b(0,0,0) +
706 a.data[ 1]*(b(0,1,0) + b(1,0,0) + b(0,0,1)) +
707 a.data[ 2]*(b(0,2,0) + b(2,0,0) + b(0,0,2)) +
708 a.data[ 3]*b(1,1,0) +
709 a.data[ 4]*(b(1,2,0) + b(2,1,0)) +
710 a.data[ 5]*b(2,2,0) +
711 a.data[ 6]*(b(0,1,1) + b(1,0,1)) +
712 a.data[ 7]*(b(0,2,1) + b(2,0,1) + b(0,1,2) + b(1,0,2)) +
713 a.data[ 8]*b(1,1,1) +
714 a.data[ 9]*(b(1,2,1) + b(2,1,1)) +
715 a.data[ 10]*b(2,2,1) +
716 a.data[ 11]*(b(0,2,2) + b(2,0,2)) +
717 a.data[ 12]*b(1,1,2) +
718 a.data[ 13]*(b(1,2,2) + b(2,1,2)) +
719 a.data[ 14]*b(2,2,2);
720
721 ret(1) +=
722 a.data[ 1]*b(0,0,0) +
723 a.data[ 4]*b(0,0,2) +
724 a.data[ 3]*b(0,0,1) +
725 a.data[ 6]*(b(0,1,0) + b(1,0,0)) +
726 a.data[ 7]*(b(0,2,0) + b(2,0,0)) +
727 a.data[ 8]*(b(1,1,0) + b(0,1,1) + b(1,0,1)) +
728 a.data[ 9]*(b(1,2,0) + b(2,1,0) + b(0,1,2) + b(1,0,2)) +
729 a.data[ 10]*b(2,2,0) +
730 a.data[ 12]*(b(0,2,1) + b(2,0,1)) +
731 a.data[ 13]*(b(0,2,2) + b(2,0,2)) +
732 a.data[ 15]*b(1,1,1) +
733 a.data[ 16]*(b(1,2,1) + b(2,1,1) + b(1,1,2)) +
734 a.data[ 17]*b(2,2,1) +
735 a.data[ 18]*(b(1,2,2) + b(2,1,2)) +
736 a.data[ 19]*b(2,2,2);
737
738 ret(2) +=
739 a.data[ 2]*b(0,0,0) +
740 a.data[ 4]*b(0,0,1) +
741 a.data[ 5]*b(0,0,2) +
742 a.data[ 7]*(b(0,1,0) + b(1,0,0)) +
743 a.data[ 9]*(b(0,1,1) + b(1,0,1)) +
744 a.data[ 10]*(b(0,1,2) + b(1,0,2)) +
745 a.data[ 11]*(b(0,2,0) + b(2,0,0)) +
746 a.data[ 12]*b(1,1,0) +
747 a.data[ 13]*(b(1,2,0) + b(2,1,0) + b(0,2,1) + b(2,0,1)) +
748 a.data[ 14]*(b(2,2,0) + b(0,2,2) + b(2,0,2)) +
749 a.data[ 16]*b(1,1,1) +
750 a.data[ 17]*b(1,1,2) +
751 a.data[ 18]*(b(1,2,1) + b(2,1,1)) +
752 a.data[ 19]*(b(2,2,1) + b(1,2,2) + b(2,1,2)) +
753 a.data[ 20]*b(2,2,2);
754
755 return ret;
756}
757
758}
759#endif
Matrix< _Scalar, _Rows, _Cols > & operator*=(const amrex::Vector< amrex::Real > &x)
Definition Eigen_Amrex.H:18
std::time_t t
AMREX_FORCE_INLINE void operator+=(const OP_CLASS &rhs)
#define INFO
Definition Util.H:20
static Matrix4< 2, Sym::MajorMinor > Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Quaternion q)
static Matrix4< 2, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Quaternion q)
static Matrix4< 2, Sym::MajorMinor > Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Eigen::Matrix3d R=Eigen::Matrix3d::Identity())
static Matrix4< 2, Sym::MajorMinor > Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
static Matrix4< 2, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Eigen::Matrix3d R=Eigen::Matrix3d::Identity())
static Matrix4< 2, Sym::MajorMinor > Zero()
static Matrix4< 2, Sym::MajorMinor > Randomize()
static Matrix4< 2, Sym::MajorMinor > Increment()
static Matrix4< 2, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
static Matrix4< 3, Sym::MajorMinor > Increment()
static Matrix4< 3, Sym::MajorMinor > Randomize()
static Matrix4< 3, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Quaternion q)
static Matrix4< 3, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
static Matrix4< 3, Sym::MajorMinor > Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Eigen::Matrix3d R=Eigen::Matrix3d::Identity())
static Matrix4< 3, Sym::MajorMinor > Zero()
static Matrix4< 3, Sym::MajorMinor > Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
static Matrix4< 3, Sym::MajorMinor > Trans(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Quaternion q)
static Matrix4< 3, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Eigen::Matrix3d R=Eigen::Matrix3d::Identity())
A collection of data types and symmetry-reduced data structures.
Definition Base.H:18
AMREX_FORCE_INLINE Quaternion operator-(const Quaternion a, const Quaternion b)
Definition Base.H:100
static Scalar Garbage
Definition Base.H:118
Sym
Definition Base.H:197
@ MajorMinor
Definition Base.H:197
amrex::Real Scalar
Definition Base.H:19
AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE Matrix4< AMREX_SPACEDIM, Sym::Diagonal > operator/(const Matrix4< AMREX_SPACEDIM, Sym::Diagonal > &a, const Set::Scalar &b)
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:20
AMREX_FORCE_INLINE bool operator==(const Quaternion a, const Quaternion b)
Definition Base.H:108
AMREX_FORCE_INLINE Quaternion operator+(const Quaternion a, const Quaternion b)
Definition Base.H:92
AMREX_FORCE_INLINE Quaternion operator*(const Set::Scalar alpha, const Quaternion b)
Definition Base.H:78
void Abort(const char *msg)
Definition Util.cpp:170
Set::Scalar Random()
Definition Set.cpp:9