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