Line data Source code
1 : #ifndef SET_MATRIX4_MAJORMINOR_H
2 : #define SET_MATRIX4_MAJORMINOR_H
3 :
4 : #include "Base.H"
5 : #include <math.h>
6 :
7 : namespace 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 : ///
13 : template<>
14 : class Matrix4<2,Sym::MajorMinor>
15 : {
16 : Scalar data[6] = {NAN,NAN,NAN,NAN,NAN,NAN};
17 :
18 : public:
19 1300 : 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 1280 : int uid = i + 2*j + 4*k + 8*l;
28 1280 : if (uid==0 ) return data[0]; // [0000]
29 1200 : else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1]; // [0001] [0010] [0100] [1000]
30 880 : else if (uid==12 || uid==3 ) return data[2]; // [0011] [1100]
31 720 : else if (uid==10 || uid==6 || uid==9 || uid==5 ) return data[3]; // [0101] [1001] [0110] [1010]
32 400 : else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[4]; // [0111] [1011] [1101] [1110]
33 80 : else if (uid==15 ) return data[5]; // [1111]
34 0 : else Util::Abort(INFO,"Index out of range");
35 0 : 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 1535104 : int uid = i + 2*j + 4*k + 8*l;
41 191314 : if (uid==0 ) return data[0]; // [0000]
42 1439160 : else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1]; // [0001] [0010] [0100] [1000]
43 1055384 : else if (uid==12 || uid==3 ) return data[2]; // [0011] [1100]
44 863496 : else if (uid==10 || uid==6 || uid==9 || uid==5 ) return data[3]; // [0101] [1001] [0110] [1010]
45 479720 : else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[4]; // [0111] [1011] [1101] [1110]
46 95944 : else if (uid==15 ) return data[5]; // [1111]
47 0 : else Util::Abort(INFO,"Index out of range");
48 0 : return Set::Garbage;
49 : }
50 0 : void Print (std::ostream& os)
51 : {
52 0 : os.precision(4);
53 0 : for (int k = 0; k < 2; k++)
54 : {
55 0 : for (int i = -1; i < 3; i++)
56 : {
57 0 : for (int l = 0; l < 2; l++)
58 : {
59 0 : if (i==-1) os << "┌ ┐ ";
60 0 : else if (i == 2) os << "└ ┘ ";
61 : else
62 : {
63 0 : os << "│ ";
64 0 : for (int j = 0; j < 2; j++)
65 : {
66 0 : const Set::Scalar &val = (*this)(i,j,k,l);
67 0 : os << std::scientific << std::setw(11) << val ; //(fabs(val)>1E-10 ? val : 0);
68 0 : os << " ";
69 : }
70 0 : os << "│ ";
71 : }
72 : }
73 0 : os<<std::endl;
74 : }
75 : }
76 0 : }
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 0 : 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 86 : static Matrix4<2,Sym::MajorMinor> Zero()
113 : {
114 86 : Matrix4<2,Sym::MajorMinor> ret;
115 602 : for (int i = 0 ; i < 6; i++) ret.data[i] = 0.0;
116 86 : return ret;
117 : }
118 1106 : static Matrix4<2,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
119 : Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
120 : {
121 1106 : Set::Matrix4<2,Sym::MajorMinor> ret;
122 : Set::Scalar Ctmp[3][3][3][3];
123 :
124 4424 : for(int i = 0; i < 3; i++)
125 13272 : for(int j = 0; j < 3; j++)
126 39816 : for(int k = 0; k < 3; k++)
127 119448 : for(int l = 0; l < 3; l++)
128 : {
129 89586 : if(i == j && j == k && k == l) Ctmp[i][j][k][l] = C11;
130 86268 : else if (i==k && j==l) Ctmp[i][j][k][l] = C44;
131 79632 : else if (i==j && k==l) Ctmp[i][j][k][l] = C12;
132 72996 : else Ctmp[i][j][k][l] = 0.0;
133 : }
134 3318 : for(int p = 0; p < 2; p++)
135 6636 : for(int q = 0; q < 2; q++)
136 13272 : for(int s = 0; s < 2; s++)
137 26544 : for(int t = 0; t < 2; t++)
138 : {
139 17696 : ret(p,q,s,t) = 0.0;
140 70784 : for(int i = 0; i < 3; i++)
141 212352 : for(int j = 0; j < 3; j++)
142 637056 : for(int k = 0; k < 3; k++)
143 1911168 : for(int l = 0; l < 3; l++)
144 2866752 : ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
145 : }
146 2212 : return ret;
147 : }
148 46 : 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 46 : Eigen::Matrix3d R;
152 46 : R = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
153 92 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
154 92 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
155 46 : return Matrix4<2,Sym::MajorMinor>::Cubic(C11,C12,C44,R);
156 : }
157 1060 : static Matrix4<2,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
158 : Set::Quaternion q)
159 : {
160 1060 : Eigen::Matrix3d R = q.normalized().toRotationMatrix();
161 1060 : 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 : };
174 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
175 : bool operator == (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b)
176 : {
177 140 : for (int i = 0 ; i < 6 ; i++) if (a.data[i] != b.data[i]) return false;
178 20 : return true;
179 : }
180 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
181 : Matrix4<2,Sym::MajorMinor> operator + (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b)
182 : {
183 0 : Matrix4<2,Sym::MajorMinor> ret;
184 0 : for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] + b.data[i];
185 0 : return ret;
186 : }
187 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
188 : Matrix4<2,Sym::MajorMinor> operator - (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b)
189 : {
190 0 : Matrix4<2,Sym::MajorMinor> ret;
191 0 : for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] - b.data[i];
192 0 : return ret;
193 : }
194 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
195 : Matrix4<2,Sym::MajorMinor> operator * (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b)
196 : {
197 8 : Matrix4<2,Sym::MajorMinor> ret;
198 56 : for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] * b;
199 8 : return ret;
200 : }
201 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
202 : Matrix4<2,Sym::MajorMinor> operator * (Set::Scalar b,Matrix4<2,Sym::MajorMinor> a)
203 : {
204 0 : return a*b;
205 : }
206 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
207 : Matrix4<2,Sym::MajorMinor> operator / (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b)
208 : {
209 0 : Matrix4<2,Sym::MajorMinor> ret;
210 0 : for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] / b;
211 0 : return ret;
212 : }
213 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
214 : Eigen::Matrix<Set::Scalar,2,2> operator * (Matrix4<2,Sym::MajorMinor> a, Eigen::Matrix<Set::Scalar,2,2> b)
215 : {
216 2020 : Eigen::Matrix<Set::Scalar,2,2> ret = Eigen::Matrix<Set::Scalar,2,2>::Zero();
217 2020 : 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 2020 : 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 2020 : 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 2020 : ret(1,0) = ret(0,1);
221 2020 : return ret;
222 : }
223 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
224 : Eigen::Matrix<amrex::Real,2,1> operator * (Matrix4<2,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,2,2>,2> b);
225 :
226 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
227 : Set::Vector operator * (Matrix4<2,Sym::MajorMinor> a, Set::Matrix3 b)
228 : {
229 0 : Set::Vector ret = Set::Vector::Zero();
230 0 : 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 0 : 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 0 : 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 : ///
247 : template<>
248 : class 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 :
254 : public:
255 1425026 : 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 7855 : int uid = i + 3*j + 9*k + 27*l;
260 : // [0000]
261 7855 : if (uid==0 ) return data[0];
262 : // [0001] [0010] [0100] [1000]
263 6400 : else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
264 : // [0002] [0020] [0200] [2000]
265 6080 : else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
266 : // [0011] [1100]
267 5760 : else if (uid==36 || uid==4 ) return data[3];
268 : // [0012] [0021] [1200] [2100]
269 5600 : else if (uid==63 || uid==45 || uid==7 || uid==5 ) return data[4];
270 : // [0022] [2200]
271 5280 : else if (uid==72 || uid==8 ) return data[5];
272 : // [0101] [0110] [1001] [1010]
273 5120 : else if (uid==30 || uid==12 || uid==28 || uid==10 ) return data[6];
274 : // [0102] [0120] [1002] [1020] [0201] [2001] [0210] [2010]
275 4800 : 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 4160 : else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[8];
278 : // [0112] [1012] [0121] [1021] [1201] [1210] [2101] [2110]
279 3840 : 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 3200 : else if (uid==75 || uid==73 || uid==35 || uid==17 ) return data[10];
282 : // [0202] [2002] [0220] [2020]
283 2880 : else if (uid==60 || uid==24 || uid==56 || uid==20 ) return data[11];
284 : // [0211] [2011] [1102] [1120]
285 2560 : else if (uid==42 || uid==58 || uid==22 || uid==38 ) return data[12];
286 : // [0212] [2012] [0221] [2021] [1202] [1220] [2102] [2120]
287 2240 : 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 1600 : else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[14];
290 : // [1111]
291 1280 : else if (uid==40 ) return data[15];
292 : // [1112] [1121] [1211] [2111]
293 1200 : else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[16];
294 : // [1122] [2211]
295 880 : else if (uid==76 || uid==44 ) return data[17];
296 : // [1212] [2112] [1221] [2121]
297 720 : else if (uid==70 || uid==52 || uid==68 || uid==50 ) return data[18];
298 : // [1222] [2122] [2212] [2221]
299 400 : else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[19];
300 : // [2222]
301 80 : else if (uid==80 ) return data[20];
302 0 : else Util::Abort(INFO,"Index out of range");
303 0 : 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 27915030 : int uid = i + 3*j + 9*k + 27*l;
309 : // [0000]
310 2604144 : if (uid==0 ) return data[0];
311 : // [0001] [0010] [0100] [1000]
312 27570400 : else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
313 : // [0002] [0020] [0200] [2000]
314 26191880 : else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
315 : // [0011] [1100]
316 24813360 : else if (uid==36 || uid==4 ) return data[3];
317 : // [0012] [0021] [1200] [2100]
318 24124100 : else if (uid==63 || uid==45 || uid==7 || uid==5 ) return data[4];
319 : // [0022] [2200]
320 22745580 : else if (uid==72 || uid==8 ) return data[5];
321 : // [0101] [0110] [1001] [1010]
322 22056320 : else if (uid==30 || uid==12 || uid==28 || uid==10 ) return data[6];
323 : // [0102] [0120] [1002] [1020] [0201] [2001] [0210] [2010]
324 20677800 : 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 17920760 : else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[8];
327 : // [0112] [1012] [0121] [1021] [1201] [1210] [2101] [2110]
328 16542240 : 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 13785200 : else if (uid==75 || uid==73 || uid==35 || uid==17 ) return data[10];
331 : // [0202] [2002] [0220] [2020]
332 12406680 : else if (uid==60 || uid==24 || uid==56 || uid==20 ) return data[11];
333 : // [0211] [2011] [1102] [1120]
334 11028160 : else if (uid==42 || uid==58 || uid==22 || uid==38 ) return data[12];
335 : // [0212] [2012] [0221] [2021] [1202] [1220] [2102] [2120]
336 9649640 : 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 6892600 : else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[14];
339 : // [1111]
340 5514080 : else if (uid==40 ) return data[15];
341 : // [1112] [1121] [1211] [2111]
342 5169450 : else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[16];
343 : // [1122] [2211]
344 3790930 : else if (uid==76 || uid==44 ) return data[17];
345 : // [1212] [2112] [1221] [2121]
346 3101670 : else if (uid==70 || uid==52 || uid==68 || uid==50 ) return data[18];
347 : // [1222] [2122] [2212] [2221]
348 1723150 : else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[19];
349 : // [2222]
350 344630 : else if (uid==80 ) return data[20];
351 0 : else Util::Abort(INFO,"Index out of range");
352 0 : return Set::Garbage;
353 : }
354 0 : void Print (std::ostream& os)
355 : {
356 0 : for (int i = 0; i < 21; i++)
357 0 : os << "i = " << i << " " << data[i] << std::endl;
358 :
359 0 : os.precision(4);
360 0 : for (int k = 0; k < 3; k++)
361 : {
362 0 : for (int i = -1; i < 4; i++)
363 : {
364 0 : for (int l = 0; l < 3; l++)
365 : {
366 0 : if (i==-1) os << "┌ ┐ ";
367 0 : else if (i == 3) os << "└ ┘ ";
368 : else
369 : {
370 0 : os << "│ ";
371 0 : for (int j = 0; j < 3; j++)
372 : {
373 0 : const Set::Scalar &val = (*this)(i,j,k,l);
374 0 : os << std::scientific << std::setw(11) << val ; //(fabs(val)>1E-10 ? val : 0);
375 0 : os << " ";
376 : }
377 0 : os << "│ ";
378 : }
379 : }
380 0 : os<<std::endl;
381 : }
382 : }
383 0 : }
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 2750 : 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 2 : static Matrix4<3,Sym::MajorMinor> Randomize()
432 : {
433 2 : Matrix4<3,Sym::MajorMinor> ret;
434 44 : for (int i = 0 ; i < 21; i++) ret.data[i] = (Util::Random());
435 2 : return ret;
436 : }
437 313 : static Matrix4<3,Sym::MajorMinor> Zero()
438 : {
439 313 : Matrix4<3,Sym::MajorMinor> ret;
440 6886 : for (int i = 0 ; i < 21; i++) ret.data[i] = 0.0;
441 313 : return ret;
442 : }
443 3906 : static Matrix4<3,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
444 : Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
445 : {
446 3906 : Set::Matrix4<3,Sym::MajorMinor> ret;
447 : Set::Scalar Ctmp[3][3][3][3];
448 :
449 15624 : for(int i = 0; i < 3; i++)
450 46872 : for(int j = 0; j < 3; j++)
451 140616 : for(int k = 0; k < 3; k++)
452 421848 : for(int l = 0; l < 3; l++)
453 : {
454 316386 : if(i == j && j == k && k == l) Ctmp[i][j][k][l] = C11;
455 304668 : else if (i==k && j==l) Ctmp[i][j][k][l] = C44;
456 281232 : else if (i==j && k==l) Ctmp[i][j][k][l] = C12;
457 257796 : else Ctmp[i][j][k][l] = 0.0;
458 : }
459 15624 : for(int p = 0; p < 3; p++)
460 46872 : for(int q = 0; q < 3; q++)
461 140616 : for(int s = 0; s < 3; s++)
462 421848 : for(int t = 0; t < 3; t++)
463 : {
464 316386 : ret(p,q,s,t) = 0.0;
465 1265544 : for(int i = 0; i < 3; i++)
466 3796632 : for(int j = 0; j < 3; j++)
467 11389896 : for(int k = 0; k < 3; k++)
468 34169688 : for(int l = 0; l < 3; l++)
469 51254532 : ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
470 : }
471 7812 : return ret;
472 : }
473 46 : 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 46 : Eigen::Matrix3d R;
477 46 : R = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
478 92 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
479 92 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
480 46 : return Matrix4<3,Sym::MajorMinor>::Cubic(C11,C12,C44,R);
481 : }
482 3860 : static Matrix4<3,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
483 : Set::Quaternion q)
484 : {
485 3860 : Eigen::Matrix3d R = q.normalized().toRotationMatrix();
486 3860 : 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 : };
498 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
499 : bool operator == (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
500 : {
501 440 : for (int i = 0 ; i < 21 ; i++) if (a.data[i] != b.data[i]) return false;
502 20 : return true;
503 : }
504 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
505 : Matrix4<3,Sym::MajorMinor> operator + (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
506 : {
507 9800 : Matrix4<3,Sym::MajorMinor> ret;
508 215600 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] + b.data[i];
509 9800 : return ret;
510 : }
511 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
512 : Matrix4<3,Sym::MajorMinor> operator - (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
513 : {
514 294597 : Matrix4<3,Sym::MajorMinor> ret;
515 6481134 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] - b.data[i];
516 294597 : return ret;
517 : }
518 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
519 : Matrix4<3,Sym::MajorMinor> operator * (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b)
520 : {
521 301330 : Matrix4<3,Sym::MajorMinor> ret;
522 6629260 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] * b;
523 301330 : return ret;
524 : }
525 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
526 : Matrix4<3,Sym::MajorMinor> operator * (Set::Scalar b,Matrix4<3,Sym::MajorMinor> a)
527 : {
528 0 : return a*b;
529 : }
530 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
531 : Matrix4<3,Sym::MajorMinor> operator / (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b)
532 : {
533 295597 : Matrix4<3,Sym::MajorMinor> ret;
534 6503134 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] / b;
535 295597 : return ret;
536 : }
537 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
538 : Eigen::Matrix<amrex::Real,3,3> operator * (Matrix4<3,Sym::MajorMinor> a, Eigen::Matrix<amrex::Real,3,3> b)
539 : {
540 1411766 : Eigen::Matrix<amrex::Real,3,3> ret = Eigen::Matrix<amrex::Real,3,3>::Zero();
541 1411766 : 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 1411766 : 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 1411766 : 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 1411766 : 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 1411766 : 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 1411766 : 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 1411766 : ret(2,1) = ret(1,2);
549 1411766 : ret(0,2) = ret(2,0);
550 1411766 : ret(1,0) = ret(0,1);
551 :
552 1411766 : return ret;
553 : }
554 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
555 : Eigen::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 : }
565 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
566 : Eigen::Matrix<amrex::Real,3,1> operator * (Matrix4<3,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,3,3>,3> b);
567 :
568 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
569 : Eigen::Matrix<amrex::Real,2,1> operator * (Matrix4<3,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,2,2>,2> b);
570 :
571 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
572 : Set::Vector operator * (Matrix4<3,Sym::MajorMinor> a, Set::Matrix3 b)
573 : {
574 : // TODO: improve efficiency of this method
575 106599 : Set::Vector ret = Set::Vector::Zero();
576 :
577 213198 : ret(0) +=
578 106599 : a.data[ 0]*b(0,0,0) +
579 319797 : a.data[ 1]*(b(0,1,0) + b(1,0,0) + b(0,0,1)) +
580 319797 : a.data[ 2]*(b(0,2,0) + b(2,0,0) + b(0,0,2)) +
581 106599 : a.data[ 3]*b(1,1,0) +
582 213198 : a.data[ 4]*(b(1,2,0) + b(2,1,0)) +
583 106599 : a.data[ 5]*b(2,2,0) +
584 213198 : a.data[ 6]*(b(0,1,1) + b(1,0,1)) +
585 426396 : a.data[ 7]*(b(0,2,1) + b(2,0,1) + b(0,1,2) + b(1,0,2)) +
586 106599 : a.data[ 8]*b(1,1,1) +
587 213198 : a.data[ 9]*(b(1,2,1) + b(2,1,1)) +
588 106599 : a.data[ 10]*b(2,2,1) +
589 213198 : a.data[ 11]*(b(0,2,2) + b(2,0,2)) +
590 106599 : a.data[ 12]*b(1,1,2) +
591 213198 : a.data[ 13]*(b(1,2,2) + b(2,1,2)) +
592 213198 : a.data[ 14]*b(2,2,2);
593 :
594 213198 : ret(1) +=
595 106599 : a.data[ 1]*b(0,0,0) +
596 106599 : a.data[ 4]*b(0,0,2) +
597 106599 : a.data[ 3]*b(0,0,1) +
598 213198 : a.data[ 6]*(b(0,1,0) + b(1,0,0)) +
599 213198 : a.data[ 7]*(b(0,2,0) + b(2,0,0)) +
600 319797 : a.data[ 8]*(b(1,1,0) + b(0,1,1) + b(1,0,1)) +
601 426396 : a.data[ 9]*(b(1,2,0) + b(2,1,0) + b(0,1,2) + b(1,0,2)) +
602 106599 : a.data[ 10]*b(2,2,0) +
603 213198 : a.data[ 12]*(b(0,2,1) + b(2,0,1)) +
604 213198 : a.data[ 13]*(b(0,2,2) + b(2,0,2)) +
605 106599 : a.data[ 15]*b(1,1,1) +
606 319797 : a.data[ 16]*(b(1,2,1) + b(2,1,1) + b(1,1,2)) +
607 106599 : a.data[ 17]*b(2,2,1) +
608 213198 : a.data[ 18]*(b(1,2,2) + b(2,1,2)) +
609 213198 : a.data[ 19]*b(2,2,2);
610 :
611 213198 : ret(2) +=
612 106599 : a.data[ 2]*b(0,0,0) +
613 106599 : a.data[ 4]*b(0,0,1) +
614 106599 : a.data[ 5]*b(0,0,2) +
615 213198 : a.data[ 7]*(b(0,1,0) + b(1,0,0)) +
616 213198 : a.data[ 9]*(b(0,1,1) + b(1,0,1)) +
617 213198 : a.data[ 10]*(b(0,1,2) + b(1,0,2)) +
618 213198 : a.data[ 11]*(b(0,2,0) + b(2,0,0)) +
619 106599 : a.data[ 12]*b(1,1,0) +
620 426396 : a.data[ 13]*(b(1,2,0) + b(2,1,0) + b(0,2,1) + b(2,0,1)) +
621 319797 : a.data[ 14]*(b(2,2,0) + b(0,2,2) + b(2,0,2)) +
622 106599 : a.data[ 16]*b(1,1,1) +
623 106599 : a.data[ 17]*b(1,1,2) +
624 213198 : a.data[ 18]*(b(1,2,1) + b(2,1,1)) +
625 319797 : a.data[ 19]*(b(2,2,1) + b(1,2,2) + b(2,1,2)) +
626 213198 : a.data[ 20]*b(2,2,2);
627 :
628 106599 : return ret;
629 : }
630 :
631 : }
632 : #endif
|
