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<>
14class Matrix4<2,Sym::MajorMinor>
15{
16 Scalar data[6] = {NAN,NAN,NAN,NAN,NAN,NAN};
17
18public:
19 AMREX_GPU_HOST_DEVICE Matrix4() {};
20 //Matrix4(const Matrix4<2,Sym::MajorMinor> &in)
21 //{
22 // for (int i = 0; i < 6; i++) data[i] = in.data[i];
23 //}
24 AMREX_FORCE_INLINE
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 }
37 AMREX_FORCE_INLINE
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
78 Set::Scalar Norm()
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
100 static Matrix4<2,Sym::MajorMinor> Increment()
101 {
102 Matrix4<2,Sym::MajorMinor> ret;
103 for (int i = 0 ; i < 6; i++) ret.data[i] = (Set::Scalar)i;
104 return ret;
105 }
106 static Matrix4<2,Sym::MajorMinor> Randomize()
107 {
108 Matrix4<2,Sym::MajorMinor> ret;
109 for (int i = 0 ; i < 6; i++) ret.data[i] = (Util::Random());
110 return ret;
111 }
112 static Matrix4<2,Sym::MajorMinor> Zero()
113 {
114 Matrix4<2,Sym::MajorMinor> ret;
115 for (int i = 0 ; i < 6; i++) ret.data[i] = 0.0;
116 return ret;
117 }
118 static Matrix4<2,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
119 Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
120 {
121 Set::Matrix4<2,Sym::MajorMinor> ret;
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 }
148 static Matrix4<2,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
149 Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
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 }
157 static Matrix4<2,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
158 Set::Quaternion q)
159 {
160 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
161 return Cubic(C11,C12,C44,R);
162 }
163
164
165 friend Eigen::Matrix<Set::Scalar,2,2> operator * (Matrix4<2,Sym::MajorMinor> a, Eigen::Matrix<Set::Scalar,2,2> b);
166 friend bool operator == (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b);
167 friend Matrix4<2,Sym::MajorMinor> operator + (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b);
168 friend Matrix4<2,Sym::MajorMinor> operator - (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b);
169 friend Matrix4<2,Sym::MajorMinor> operator * (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b);
170 friend Matrix4<2,Sym::MajorMinor> operator * (Set::Scalar b, Matrix4<2,Sym::MajorMinor> a);
171 friend Matrix4<2,Sym::MajorMinor> operator / (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b);
172 friend Set::Vector operator * (Matrix4<2,Sym::MajorMinor> a, Set::Matrix3 b);
173};
174AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
175bool operator == (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b)
176{
177 for (int i = 0 ; i < 6 ; i++) if (a.data[i] != b.data[i]) return false;
178 return true;
179}
180AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
181Matrix4<2,Sym::MajorMinor> operator + (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b)
182{
183 Matrix4<2,Sym::MajorMinor> ret;
184 for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] + b.data[i];
185 return ret;
186}
187AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
188Matrix4<2,Sym::MajorMinor> operator - (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b)
189{
190 Matrix4<2,Sym::MajorMinor> ret;
191 for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] - b.data[i];
192 return ret;
193}
194AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
195Matrix4<2,Sym::MajorMinor> operator * (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b)
196{
197 Matrix4<2,Sym::MajorMinor> ret;
198 for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] * b;
199 return ret;
200}
201AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
206AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
207Matrix4<2,Sym::MajorMinor> operator / (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b)
208{
209 Matrix4<2,Sym::MajorMinor> ret;
210 for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] / b;
211 return ret;
212}
213AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
214Eigen::Matrix<Set::Scalar,2,2> operator * (Matrix4<2,Sym::MajorMinor> a, Eigen::Matrix<Set::Scalar,2,2> b)
215{
216 Eigen::Matrix<Set::Scalar,2,2> ret = Eigen::Matrix<Set::Scalar,2,2>::Zero();
217 ret(0,0) = a.data[0]*b(0,0) + a.data[1]*(b(0,1) + b(1,0)) + a.data[2]*b(1,1);
218 ret(0,1) = a.data[1]*b(0,0) + a.data[3]*(b(0,1) + b(1,0)) + a.data[4]*b(1,1);
219 ret(1,1) = a.data[2]*b(0,0) + a.data[4]*(b(0,1) + b(1,0)) + a.data[5]*b(1,1);
220 ret(1,0) = ret(0,1);
221 return ret;
222}
223AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
224Eigen::Matrix<amrex::Real,2,1> operator * (Matrix4<2,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,2,2>,2> b);
225
226AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
227Set::Vector operator * (Matrix4<2,Sym::MajorMinor> a, Set::Matrix3 b)
228{
229 Set::Vector ret = Set::Vector::Zero();
230 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);
231 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);
232 return ret;
233}
234
235/// \brief Data structure for a 4th order 3D tensor with major and minor symmetry
236///
237/// Let the tensor
238/// \f$\mathbb{C}\in\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{R}^3\f$
239/// be have major symmetry (\f$\mathbb{C}_{ijkl}=\mathbb{C}_{klij}\f$)
240/// and minor symmetry (\f$\mathbb{C}_{ijkl}=\mathbb{C}_{jikl}=\mathbb{C}_{ijlk}\f$)
241/// Then there are only 21 unique elements (rather than 81).
242///
243/// This object acts like a 4D array such that `C(i,j,k,l)` returns the corresponding
244/// element, but symmetry is always obeyed. This allows the user code to be much prettier,
245/// while maintaining a relatively small object size.
246///
247template<>
248class Matrix4<3,Sym::MajorMinor>
249{
250 Scalar data[21] = { NAN,NAN,NAN,NAN,NAN,NAN,NAN,
251 NAN,NAN,NAN,NAN,NAN,NAN,NAN,
252 NAN,NAN,NAN,NAN,NAN,NAN,NAN };
253
254public:
255 AMREX_GPU_HOST_DEVICE Matrix4() {};
256 AMREX_FORCE_INLINE
257 const Scalar & operator () (const int i, const int j, const int k, const int l) const
258 {
259 int uid = i + 3*j + 9*k + 27*l;
260 // [0000]
261 if (uid==0 ) return data[0];
262 // [0001] [0010] [0100] [1000]
263 else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
264 // [0002] [0020] [0200] [2000]
265 else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
266 // [0011] [1100]
267 else if (uid==36 || uid==4 ) return data[3];
268 // [0012] [0021] [1200] [2100]
269 else if (uid==63 || uid==45 || uid==7 || uid==5 ) return data[4];
270 // [0022] [2200]
271 else if (uid==72 || uid==8 ) return data[5];
272 // [0101] [0110] [1001] [1010]
273 else if (uid==30 || uid==12 || uid==28 || uid==10 ) return data[6];
274 // [0102] [0120] [1002] [1020] [0201] [2001] [0210] [2010]
275 else if (uid==57 || uid==21 || uid==33 || uid==15 || uid==55 || uid==19 || uid==29 || uid==11 ) return data[7];
276 // [0111] [1011] [1101] [1110]
277 else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[8];
278 // [0112] [1012] [0121] [1021] [1201] [1210] [2101] [2110]
279 else if (uid==66 || uid==48 || uid==64 || uid==46 || uid==34 || uid==16 || uid==32 || uid==14 ) return data[9];
280 // [0122] [1022] [2201] [2210]
281 else if (uid==75 || uid==73 || uid==35 || uid==17 ) return data[10];
282 // [0202] [2002] [0220] [2020]
283 else if (uid==60 || uid==24 || uid==56 || uid==20 ) return data[11];
284 // [0211] [2011] [1102] [1120]
285 else if (uid==42 || uid==58 || uid==22 || uid==38 ) return data[12];
286 // [0212] [2012] [0221] [2021] [1202] [1220] [2102] [2120]
287 else if (uid==69 || uid==51 || uid==61 || uid==25 || uid==65 || uid==47 || uid==59 || uid==23 ) return data[13];
288 // [0222] [2022] [2202] [2220]
289 else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[14];
290 // [1111]
291 else if (uid==40 ) return data[15];
292 // [1112] [1121] [1211] [2111]
293 else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[16];
294 // [1122] [2211]
295 else if (uid==76 || uid==44 ) return data[17];
296 // [1212] [2112] [1221] [2121]
297 else if (uid==70 || uid==52 || uid==68 || uid==50 ) return data[18];
298 // [1222] [2122] [2212] [2221]
299 else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[19];
300 // [2222]
301 else if (uid==80 ) return data[20];
302 else Util::Abort(INFO,"Index out of range");
303 return Set::Garbage;
304 }
305 AMREX_FORCE_INLINE
306 Scalar & operator () (const int i, const int j, const int k, const int l)
307 {
308 int uid = i + 3*j + 9*k + 27*l;
309 // [0000]
310 if (uid==0 ) return data[0];
311 // [0001] [0010] [0100] [1000]
312 else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
313 // [0002] [0020] [0200] [2000]
314 else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
315 // [0011] [1100]
316 else if (uid==36 || uid==4 ) return data[3];
317 // [0012] [0021] [1200] [2100]
318 else if (uid==63 || uid==45 || uid==7 || uid==5 ) return data[4];
319 // [0022] [2200]
320 else if (uid==72 || uid==8 ) return data[5];
321 // [0101] [0110] [1001] [1010]
322 else if (uid==30 || uid==12 || uid==28 || uid==10 ) return data[6];
323 // [0102] [0120] [1002] [1020] [0201] [2001] [0210] [2010]
324 else if (uid==57 || uid==21 || uid==33 || uid==15 || uid==55 || uid==19 || uid==29 || uid==11 ) return data[7];
325 // [0111] [1011] [1101] [1110]
326 else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[8];
327 // [0112] [1012] [0121] [1021] [1201] [1210] [2101] [2110]
328 else if (uid==66 || uid==48 || uid==64 || uid==46 || uid==34 || uid==16 || uid==32 || uid==14 ) return data[9];
329 // [0122] [1022] [2201] [2210]
330 else if (uid==75 || uid==73 || uid==35 || uid==17 ) return data[10];
331 // [0202] [2002] [0220] [2020]
332 else if (uid==60 || uid==24 || uid==56 || uid==20 ) return data[11];
333 // [0211] [2011] [1102] [1120]
334 else if (uid==42 || uid==58 || uid==22 || uid==38 ) return data[12];
335 // [0212] [2012] [0221] [2021] [1202] [1220] [2102] [2120]
336 else if (uid==69 || uid==51 || uid==61 || uid==25 || uid==65 || uid==47 || uid==59 || uid==23 ) return data[13];
337 // [0222] [2022] [2202] [2220]
338 else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[14];
339 // [1111]
340 else if (uid==40 ) return data[15];
341 // [1112] [1121] [1211] [2111]
342 else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[16];
343 // [1122] [2211]
344 else if (uid==76 || uid==44 ) return data[17];
345 // [1212] [2112] [1221] [2121]
346 else if (uid==70 || uid==52 || uid==68 || uid==50 ) return data[18];
347 // [1222] [2122] [2212] [2221]
348 else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[19];
349 // [2222]
350 else if (uid==80 ) return data[20];
351 else Util::Abort(INFO,"Index out of range");
352 return Set::Garbage;
353 }
354 void Print (std::ostream& os)
355 {
356 for (int i = 0; i < 21; i++)
357 os << "i = " << i << " " << data[i] << std::endl;
358
359 os.precision(4);
360 for (int k = 0; k < 3; k++)
361 {
362 for (int i = -1; i < 4; i++)
363 {
364 for (int l = 0; l < 3; l++)
365 {
366 if (i==-1) os << "┌ ┐ ";
367 else if (i == 3) os << "└ ┘ ";
368 else
369 {
370 os << "│ ";
371 for (int j = 0; j < 3; j++)
372 {
373 const Set::Scalar &val = (*this)(i,j,k,l);
374 os << std::scientific << std::setw(11) << val ; //(fabs(val)>1E-10 ? val : 0);
375 os << " ";
376 }
377 os << "│ ";
378 }
379 }
380 os<<std::endl;
381 }
382 }
383 }
384
385 Set::Scalar Norm()
386 {
387 Set::Scalar retsq = 0;
388 for (int i = 0; i < 21; i++) retsq += data[i];
389 return std::sqrt(retsq);
390 }
391 bool contains_nan() const
392 {
393 if (std::isnan(data[ 0])) return true;
394 if (std::isnan(data[ 1])) return true;
395 if (std::isnan(data[ 2])) return true;
396 if (std::isnan(data[ 3])) return true;
397 if (std::isnan(data[ 4])) return true;
398 if (std::isnan(data[ 5])) return true;
399 if (std::isnan(data[ 6])) return true;
400 if (std::isnan(data[ 7])) return true;
401 if (std::isnan(data[ 8])) return true;
402 if (std::isnan(data[ 9])) return true;
403 if (std::isnan(data[10])) return true;
404 if (std::isnan(data[11])) return true;
405 if (std::isnan(data[12])) return true;
406 if (std::isnan(data[13])) return true;
407 if (std::isnan(data[14])) return true;
408 if (std::isnan(data[15])) return true;
409 if (std::isnan(data[16])) return true;
410 if (std::isnan(data[17])) return true;
411 if (std::isnan(data[18])) return true;
412 if (std::isnan(data[19])) return true;
413 if (std::isnan(data[20])) return true;
414 return false;
415 }
416
417 //AMREX_GPU_HOST_DEVICE void operator = (Matrix4<3,Sym::MajorMinor> a) {for (int i = 0; i < 21; i++) data[i] = a.data[i];}
418 AMREX_GPU_HOST_DEVICE void operator += (Matrix4<3,Sym::MajorMinor> a) {for (int i = 0; i < 21; i++) data[i] += a.data[i];}
419 AMREX_GPU_HOST_DEVICE void operator -= (Matrix4<3,Sym::MajorMinor> a) {for (int i = 0; i < 21; i++) data[i] -= a.data[i];}
420 AMREX_GPU_HOST_DEVICE void operator *= (Matrix4<3,Sym::MajorMinor> a) {for (int i = 0; i < 21; i++) data[i] *= a.data[i];}
421 AMREX_GPU_HOST_DEVICE void operator /= (Matrix4<3,Sym::MajorMinor> a) {for (int i = 0; i < 21; i++) data[i] /= a.data[i];}
422 AMREX_GPU_HOST_DEVICE void operator *= (Set::Scalar alpha) {for (int i = 0; i < 21; i++) data[i] *= alpha;}
423 AMREX_GPU_HOST_DEVICE void operator /= (Set::Scalar alpha) {for (int i = 0; i < 21; i++) data[i] /= alpha;}
424
425 static Matrix4<3,Sym::MajorMinor> Increment()
426 {
427 Matrix4<3,Sym::MajorMinor> ret;
428 for (int i = 0 ; i < 21; i++) ret.data[i] = (Set::Scalar)i;
429 return ret;
430 }
431 static Matrix4<3,Sym::MajorMinor> Randomize()
432 {
433 Matrix4<3,Sym::MajorMinor> ret;
434 for (int i = 0 ; i < 21; i++) ret.data[i] = (Util::Random());
435 return ret;
436 }
437 static Matrix4<3,Sym::MajorMinor> Zero()
438 {
439 Matrix4<3,Sym::MajorMinor> ret;
440 for (int i = 0 ; i < 21; i++) ret.data[i] = 0.0;
441 return ret;
442 }
443 static Matrix4<3,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
444 Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
445 {
446 Set::Matrix4<3,Sym::MajorMinor> ret;
447 Set::Scalar Ctmp[3][3][3][3];
448
449 for(int i = 0; i < 3; i++)
450 for(int j = 0; j < 3; j++)
451 for(int k = 0; k < 3; k++)
452 for(int l = 0; l < 3; l++)
453 {
454 if(i == j && j == k && k == l) Ctmp[i][j][k][l] = C11;
455 else if (i==k && j==l) Ctmp[i][j][k][l] = C44;
456 else if (i==j && k==l) Ctmp[i][j][k][l] = C12;
457 else Ctmp[i][j][k][l] = 0.0;
458 }
459 for(int p = 0; p < 3; p++)
460 for(int q = 0; q < 3; q++)
461 for(int s = 0; s < 3; s++)
462 for(int t = 0; t < 3; t++)
463 {
464 ret(p,q,s,t) = 0.0;
465 for(int i = 0; i < 3; i++)
466 for(int j = 0; j < 3; j++)
467 for(int k = 0; k < 3; k++)
468 for(int l = 0; l < 3; l++)
469 ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
470 }
471 return ret;
472 }
473 static Matrix4<3,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
474 Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
475 {
476 Eigen::Matrix3d R;
477 R = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
478 Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
479 Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
480 return Matrix4<3,Sym::MajorMinor>::Cubic(C11,C12,C44,R);
481 }
482 static Matrix4<3,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
483 Set::Quaternion q)
484 {
485 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
486 return Cubic(C11,C12,C44,R);
487 }
488
489 friend bool operator == (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b);
490 friend Matrix4<3,Sym::MajorMinor> operator + (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b);
491 friend Matrix4<3,Sym::MajorMinor> operator - (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b);
492 friend Matrix4<3,Sym::MajorMinor> operator * (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b);
493 friend Matrix4<3,Sym::MajorMinor> operator * (Set::Scalar b, Matrix4<3,Sym::MajorMinor> a);
494 friend Matrix4<3,Sym::MajorMinor> operator / (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b);
495 friend Set::Vector operator * (Matrix4<3,Sym::MajorMinor> a, Set::Matrix3 b);
496 friend Eigen::Matrix<amrex::Real,3,3> operator * (Matrix4<3,Sym::MajorMinor> a, Eigen::Matrix<amrex::Real,3,3> b);
497};
498AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
499bool operator == (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
500{
501 for (int i = 0 ; i < 21 ; i++) if (a.data[i] != b.data[i]) return false;
502 return true;
503}
504AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
505Matrix4<3,Sym::MajorMinor> operator + (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
506{
507 Matrix4<3,Sym::MajorMinor> ret;
508 for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] + b.data[i];
509 return ret;
510}
511AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
512Matrix4<3,Sym::MajorMinor> operator - (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
513{
514 Matrix4<3,Sym::MajorMinor> ret;
515 for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] - b.data[i];
516 return ret;
517}
518AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
519Matrix4<3,Sym::MajorMinor> operator * (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b)
520{
521 Matrix4<3,Sym::MajorMinor> ret;
522 for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] * b;
523 return ret;
524}
525AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
530AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
531Matrix4<3,Sym::MajorMinor> operator / (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b)
532{
533 Matrix4<3,Sym::MajorMinor> ret;
534 for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] / b;
535 return ret;
536}
537AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
538Eigen::Matrix<amrex::Real,3,3> operator * (Matrix4<3,Sym::MajorMinor> a, Eigen::Matrix<amrex::Real,3,3> b)
539{
540 Eigen::Matrix<amrex::Real,3,3> ret = Eigen::Matrix<amrex::Real,3,3>::Zero();
541 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);
542 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);
543 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);
544 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);
545 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);
546 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);
547
548 ret(2,1) = ret(1,2);
549 ret(0,2) = ret(2,0);
550 ret(1,0) = ret(0,1);
551
552 return ret;
553}
554AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
555Eigen::Matrix<amrex::Real,2,2> operator * (Matrix4<3,Sym::MajorMinor> a, Eigen::Matrix<amrex::Real,2,2> b)
556{
557 Eigen::Matrix<amrex::Real,2,2> ret = Eigen::Matrix<amrex::Real,2,2>::Zero();
558 for (int i = 0; i < 2; i++)
559 for (int j = 0; j < 2; j++)
560 for (int k = 0; k < 2; k++)
561 for (int l = 0; l < 2; l++)
562 ret(i,j) += a(i,j,k,l)*b(k,l);
563 return ret;
564}
565AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
566Eigen::Matrix<amrex::Real,3,1> operator * (Matrix4<3,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,3,3>,3> b);
567
568AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
569Eigen::Matrix<amrex::Real,2,1> operator * (Matrix4<3,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,2,2>,2> b);
570
571AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
572Set::Vector operator * (Matrix4<3,Sym::MajorMinor> a, Set::Matrix3 b)
573{
574 // TODO: improve efficiency of this method
575 Set::Vector ret = Set::Vector::Zero();
576
577 ret(0) +=
578 a.data[ 0]*b(0,0,0) +
579 a.data[ 1]*(b(0,1,0) + b(1,0,0) + b(0,0,1)) +
580 a.data[ 2]*(b(0,2,0) + b(2,0,0) + b(0,0,2)) +
581 a.data[ 3]*b(1,1,0) +
582 a.data[ 4]*(b(1,2,0) + b(2,1,0)) +
583 a.data[ 5]*b(2,2,0) +
584 a.data[ 6]*(b(0,1,1) + b(1,0,1)) +
585 a.data[ 7]*(b(0,2,1) + b(2,0,1) + b(0,1,2) + b(1,0,2)) +
586 a.data[ 8]*b(1,1,1) +
587 a.data[ 9]*(b(1,2,1) + b(2,1,1)) +
588 a.data[ 10]*b(2,2,1) +
589 a.data[ 11]*(b(0,2,2) + b(2,0,2)) +
590 a.data[ 12]*b(1,1,2) +
591 a.data[ 13]*(b(1,2,2) + b(2,1,2)) +
592 a.data[ 14]*b(2,2,2);
593
594 ret(1) +=
595 a.data[ 1]*b(0,0,0) +
596 a.data[ 4]*b(0,0,2) +
597 a.data[ 3]*b(0,0,1) +
598 a.data[ 6]*(b(0,1,0) + b(1,0,0)) +
599 a.data[ 7]*(b(0,2,0) + b(2,0,0)) +
600 a.data[ 8]*(b(1,1,0) + b(0,1,1) + b(1,0,1)) +
601 a.data[ 9]*(b(1,2,0) + b(2,1,0) + b(0,1,2) + b(1,0,2)) +
602 a.data[ 10]*b(2,2,0) +
603 a.data[ 12]*(b(0,2,1) + b(2,0,1)) +
604 a.data[ 13]*(b(0,2,2) + b(2,0,2)) +
605 a.data[ 15]*b(1,1,1) +
606 a.data[ 16]*(b(1,2,1) + b(2,1,1) + b(1,1,2)) +
607 a.data[ 17]*b(2,2,1) +
608 a.data[ 18]*(b(1,2,2) + b(2,1,2)) +
609 a.data[ 19]*b(2,2,2);
610
611 ret(2) +=
612 a.data[ 2]*b(0,0,0) +
613 a.data[ 4]*b(0,0,1) +
614 a.data[ 5]*b(0,0,2) +
615 a.data[ 7]*(b(0,1,0) + b(1,0,0)) +
616 a.data[ 9]*(b(0,1,1) + b(1,0,1)) +
617 a.data[ 10]*(b(0,1,2) + b(1,0,2)) +
618 a.data[ 11]*(b(0,2,0) + b(2,0,0)) +
619 a.data[ 12]*b(1,1,0) +
620 a.data[ 13]*(b(1,2,0) + b(2,1,0) + b(0,2,1) + b(2,0,1)) +
621 a.data[ 14]*(b(2,2,0) + b(0,2,2) + b(2,0,2)) +
622 a.data[ 16]*b(1,1,1) +
623 a.data[ 17]*b(1,1,2) +
624 a.data[ 18]*(b(1,2,1) + b(2,1,1)) +
625 a.data[ 19]*(b(2,2,1) + b(1,2,2) + b(2,1,2)) +
626 a.data[ 20]*b(2,2,2);
627
628 return ret;
629}
630
631}
632#endif
operator*=
Matrix< _Scalar, _Rows, _Cols > & operator*=(const amrex::Vector< amrex::Real > &x)
Definition Eigen_Amrex.H:18
t
std::time_t t
Definition FileNameParse.cpp:9
operator+=
AMREX_FORCE_INLINE void operator+=(const OP_CLASS &rhs)
Definition InClassOperators.H:2
INFO
#define INFO
Definition Util.H:20
Set::Matrix4< 2, Sym::MajorMinor >::Cubic
static Matrix4< 2, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Quaternion q)
Definition Matrix4_MajorMinor.H:157
Set::Matrix4< 2, Sym::MajorMinor >::Cubic
static Matrix4< 2, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Eigen::Matrix3d R=Eigen::Matrix3d::Identity())
Definition Matrix4_MajorMinor.H:118
Set::Matrix4< 2, Sym::MajorMinor >::Zero
static Matrix4< 2, Sym::MajorMinor > Zero()
Definition Matrix4_MajorMinor.H:112
Set::Matrix4< 2, Sym::MajorMinor >::Randomize
static Matrix4< 2, Sym::MajorMinor > Randomize()
Definition Matrix4_MajorMinor.H:106
Set::Matrix4< 2, Sym::MajorMinor >::Increment
static Matrix4< 2, Sym::MajorMinor > Increment()
Definition Matrix4_MajorMinor.H:100
Set::Matrix4< 2, Sym::MajorMinor >::Cubic
static Matrix4< 2, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
Definition Matrix4_MajorMinor.H:148
Set::Matrix4< 3, Sym::MajorMinor >::Increment
static Matrix4< 3, Sym::MajorMinor > Increment()
Definition Matrix4_MajorMinor.H:425
Set::Matrix4< 3, Sym::MajorMinor >::Randomize
static Matrix4< 3, Sym::MajorMinor > Randomize()
Definition Matrix4_MajorMinor.H:431
Set::Matrix4< 3, Sym::MajorMinor >::Cubic
static Matrix4< 3, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Quaternion q)
Definition Matrix4_MajorMinor.H:482
Set::Matrix4< 3, Sym::MajorMinor >::Cubic
static Matrix4< 3, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
Definition Matrix4_MajorMinor.H:473
Set::Matrix4< 3, Sym::MajorMinor >::Zero
static Matrix4< 3, Sym::MajorMinor > Zero()
Definition Matrix4_MajorMinor.H:437
Set::Matrix4< 3, Sym::MajorMinor >::Cubic
static Matrix4< 3, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Eigen::Matrix3d R=Eigen::Matrix3d::Identity())
Definition Matrix4_MajorMinor.H:443
Set
A collection of data types and symmetry-reduced data structures.
Definition Base.H:18
Set::operator-
AMREX_FORCE_INLINE Quaternion operator-(const Quaternion a, const Quaternion b)
Definition Base.H:100
Set::Garbage
static Scalar Garbage
Definition Base.H:118
Set::Sym
Sym
Definition Base.H:197
Set::MajorMinor
@ MajorMinor
Definition Base.H:197
Set::Scalar
amrex::Real Scalar
Definition Base.H:19
Set::operator/
AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE Matrix4< AMREX_SPACEDIM, Sym::Diagonal > operator/(const Matrix4< AMREX_SPACEDIM, Sym::Diagonal > &a, const Set::Scalar &b)
Definition Matrix4_Diagonal.H:130
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:20
Set::operator==
AMREX_FORCE_INLINE bool operator==(const Quaternion a, const Quaternion b)
Definition Base.H:108
Set::operator+
AMREX_FORCE_INLINE Quaternion operator+(const Quaternion a, const Quaternion b)
Definition Base.H:92
Set::operator*
AMREX_FORCE_INLINE Quaternion operator*(const Set::Scalar alpha, const Quaternion b)
Definition Base.H:78
Util::Abort
void Abort(const char *msg)
Definition Util.cpp:170
Util::Random
Set::Scalar Random()
Definition Set.cpp:9
Generated by doxygen 1.9.8