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 1360 : 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 1600 : int uid = i + 2*j + 4*k + 8*l;
28 1600 : if (uid==0 ) return data[0]; // [0000]
29 1500 : else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1]; // [0001] [0010] [0100] [1000]
30 1100 : else if (uid==12 || uid==3 ) return data[2]; // [0011] [1100]
31 900 : else if (uid==10 || uid==6 || uid==9 || uid==5 ) return data[3]; // [0101] [1001] [0110] [1010]
32 500 : else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[4]; // [0111] [1011] [1101] [1110]
33 100 : 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 1565280 : int uid = i + 2*j + 4*k + 8*l;
41 193545 : if (uid==0 ) return data[0]; // [0000]
42 1467450 : else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1]; // [0001] [0010] [0100] [1000]
43 1076130 : else if (uid==12 || uid==3 ) return data[2]; // [0011] [1100]
44 880470 : else if (uid==10 || uid==6 || uid==9 || uid==5 ) return data[3]; // [0101] [1001] [0110] [1010]
45 489150 : else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[4]; // [0111] [1011] [1101] [1110]
46 97830 : 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 96 : static Matrix4<2,Sym::MajorMinor> Zero()
113 : {
114 96 : Matrix4<2,Sym::MajorMinor> ret;
115 672 : for (int i = 0 ; i < 6; i++) ret.data[i] = 0.0;
116 96 : 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 23 : 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 : {
167 23 : Set::Matrix4<2,Sym::MajorMinor> ret;
168 23 : Set::Scalar Ctmp[3][3][3][3] = {{{{0.0}}}};
169 :
170 92 : for(int i = 0; i < 3; i++)
171 276 : for(int j = 0; j < 3; j++)
172 828 : for(int k = 0; k < 3; k++)
173 2484 : for(int l = 0; l < 3; l++)
174 : {
175 1863 : if(i==j && j==k && k==l)
176 : {
177 69 : if (i < 2) Ctmp[i][j][k][l] = C11;
178 23 : else Ctmp[i][j][k][l] = C33;
179 : }
180 1794 : else if(i==j && k==l)
181 : {
182 138 : if (i<2 && k<2 && i!=k) Ctmp[i][j][k][l] = C12;
183 92 : else if ((i<2 && k==2) || (i==2 && k<2)) Ctmp[i][j][k][l] = C13;
184 : }
185 1656 : else if((i==k && j==l) && (i!=j))
186 : {
187 138 : Ctmp[i][j][k][l] = C44;
188 : }
189 : else
190 : {
191 1518 : Ctmp[i][j][k][l] = 0.0;
192 : }
193 : }
194 :
195 : // C66 = 0.5*(C11 - C12)
196 23 : Ctmp[0][1][0][1] = 0.5 * (C11 - C12);
197 23 : Ctmp[1][0][1][0] = 0.5 * (C11 - C12);
198 :
199 69 : for(int p = 0; p < 2; p++)
200 138 : for(int q = 0; q < 2; q++)
201 276 : for(int s = 0; s < 2; s++)
202 552 : for(int t = 0; t < 2; t++)
203 : {
204 368 : ret(p,q,s,t) = 0.0;
205 1472 : for(int i = 0; i < 3; i++)
206 4416 : for(int j = 0; j < 3; j++)
207 13248 : for(int k = 0; k < 3; k++)
208 39744 : for(int l = 0; l < 3; l++)
209 59616 : ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
210 : }
211 46 : return ret;
212 : }
213 22 : 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)
214 : {
215 : Eigen::Quaterniond q =
216 22 : Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
217 22 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
218 44 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
219 22 : Eigen::Matrix3d R = q.toRotationMatrix();
220 22 : return Matrix4<2,Sym::MajorMinor>::Transverse(C11,C12,C13,C33,C44,R);
221 : }
222 : static Matrix4<2,Sym::MajorMinor> Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Quaternion q)
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);
230 : friend bool operator == (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b);
231 : friend Matrix4<2,Sym::MajorMinor> operator + (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b);
232 : friend Matrix4<2,Sym::MajorMinor> operator - (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b);
233 : friend Matrix4<2,Sym::MajorMinor> operator * (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b);
234 : friend Matrix4<2,Sym::MajorMinor> operator * (Set::Scalar b, Matrix4<2,Sym::MajorMinor> a);
235 : friend Matrix4<2,Sym::MajorMinor> operator / (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b);
236 : friend Set::Vector operator * (Matrix4<2,Sym::MajorMinor> a, Set::Matrix3 b);
237 : };
238 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
239 : bool operator == (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b)
240 : {
241 175 : for (int i = 0 ; i < 6 ; i++) if (a.data[i] != b.data[i]) return false;
242 25 : return true;
243 : }
244 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
245 : Matrix4<2,Sym::MajorMinor> operator + (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b)
246 : {
247 0 : Matrix4<2,Sym::MajorMinor> ret;
248 0 : for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] + b.data[i];
249 0 : return ret;
250 : }
251 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
252 : Matrix4<2,Sym::MajorMinor> operator - (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b)
253 : {
254 0 : Matrix4<2,Sym::MajorMinor> ret;
255 0 : for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] - b.data[i];
256 0 : return ret;
257 : }
258 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
259 : Matrix4<2,Sym::MajorMinor> operator * (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b)
260 : {
261 10 : Matrix4<2,Sym::MajorMinor> ret;
262 70 : for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] * b;
263 10 : return ret;
264 : }
265 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
266 : Matrix4<2,Sym::MajorMinor> operator * (Set::Scalar b,Matrix4<2,Sym::MajorMinor> a)
267 : {
268 0 : return a*b;
269 : }
270 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
271 : Matrix4<2,Sym::MajorMinor> operator / (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b)
272 : {
273 0 : Matrix4<2,Sym::MajorMinor> ret;
274 0 : for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] / b;
275 0 : return ret;
276 : }
277 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
278 : Eigen::Matrix<Set::Scalar,2,2> operator * (Matrix4<2,Sym::MajorMinor> a, Eigen::Matrix<Set::Scalar,2,2> b)
279 : {
280 2430 : Eigen::Matrix<Set::Scalar,2,2> ret = Eigen::Matrix<Set::Scalar,2,2>::Zero();
281 2430 : 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 2430 : 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 2430 : 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 2430 : ret(1,0) = ret(0,1);
285 2430 : return ret;
286 : }
287 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
288 : Eigen::Matrix<amrex::Real,2,1> operator * (Matrix4<2,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,2,2>,2> b);
289 :
290 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
291 : Set::Vector operator * (Matrix4<2,Sym::MajorMinor> a, Set::Matrix3 b)
292 : {
293 0 : Set::Vector ret = Set::Vector::Zero();
294 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);
295 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);
296 0 : 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 : ///
311 : template<>
312 : class Matrix4<3,Sym::MajorMinor>
313 : {
314 : Scalar data[21] = { NAN,NAN,NAN,NAN,NAN,NAN,NAN,
315 : NAN,NAN,NAN,NAN,NAN,NAN,NAN,
316 : NAN,NAN,NAN,NAN,NAN,NAN,NAN };
317 :
318 : public:
319 1425086 : AMREX_GPU_HOST_DEVICE Matrix4() {};
320 : AMREX_FORCE_INLINE
321 : const Scalar & operator () (const int i, const int j, const int k, const int l) const
322 : {
323 9475 : int uid = i + 3*j + 9*k + 27*l;
324 : // [0000]
325 9475 : if (uid==0 ) return data[0];
326 : // [0001] [0010] [0100] [1000]
327 8000 : else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
328 : // [0002] [0020] [0200] [2000]
329 7600 : else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
330 : // [0011] [1100]
331 7200 : else if (uid==36 || uid==4 ) return data[3];
332 : // [0012] [0021] [1200] [2100]
333 7000 : else if (uid==63 || uid==45 || uid==7 || uid==5 ) return data[4];
334 : // [0022] [2200]
335 6600 : else if (uid==72 || uid==8 ) return data[5];
336 : // [0101] [0110] [1001] [1010]
337 6400 : else if (uid==30 || uid==12 || uid==28 || uid==10 ) return data[6];
338 : // [0102] [0120] [1002] [1020] [0201] [2001] [0210] [2010]
339 6000 : 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 5200 : else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[8];
342 : // [0112] [1012] [0121] [1021] [1201] [1210] [2101] [2110]
343 4800 : 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 4000 : else if (uid==75 || uid==73 || uid==35 || uid==17 ) return data[10];
346 : // [0202] [2002] [0220] [2020]
347 3600 : else if (uid==60 || uid==24 || uid==56 || uid==20 ) return data[11];
348 : // [0211] [2011] [1102] [1120]
349 3200 : else if (uid==42 || uid==58 || uid==22 || uid==38 ) return data[12];
350 : // [0212] [2012] [0221] [2021] [1202] [1220] [2102] [2120]
351 2800 : 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 2000 : else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[14];
354 : // [1111]
355 1600 : else if (uid==40 ) return data[15];
356 : // [1112] [1121] [1211] [2111]
357 1500 : else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[16];
358 : // [1122] [2211]
359 1100 : else if (uid==76 || uid==44 ) return data[17];
360 : // [1212] [2112] [1221] [2121]
361 900 : else if (uid==70 || uid==52 || uid==68 || uid==50 ) return data[18];
362 : // [1222] [2122] [2212] [2221]
363 500 : else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[19];
364 : // [2222]
365 100 : else if (uid==80 ) return data[20];
366 0 : else Util::Abort(INFO,"Index out of range");
367 0 : return Set::Garbage;
368 : }
369 : AMREX_FORCE_INLINE
370 : Scalar & operator () (const int i, const int j, const int k, const int l)
371 : {
372 28067796 : int uid = i + 3*j + 9*k + 27*l;
373 : // [0000]
374 2607870 : if (uid==0 ) return data[0];
375 : // [0001] [0010] [0100] [1000]
376 27721280 : else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
377 : // [0002] [0020] [0200] [2000]
378 26335216 : else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
379 : // [0011] [1100]
380 24949152 : else if (uid==36 || uid==4 ) return data[3];
381 : // [0012] [0021] [1200] [2100]
382 24256120 : else if (uid==63 || uid==45 || uid==7 || uid==5 ) return data[4];
383 : // [0022] [2200]
384 22870056 : else if (uid==72 || uid==8 ) return data[5];
385 : // [0101] [0110] [1001] [1010]
386 22177024 : else if (uid==30 || uid==12 || uid==28 || uid==10 ) return data[6];
387 : // [0102] [0120] [1002] [1020] [0201] [2001] [0210] [2010]
388 20790960 : 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 18018832 : else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[8];
391 : // [0112] [1012] [0121] [1021] [1201] [1210] [2101] [2110]
392 16632768 : 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 13860640 : else if (uid==75 || uid==73 || uid==35 || uid==17 ) return data[10];
395 : // [0202] [2002] [0220] [2020]
396 12474576 : else if (uid==60 || uid==24 || uid==56 || uid==20 ) return data[11];
397 : // [0211] [2011] [1102] [1120]
398 11088512 : else if (uid==42 || uid==58 || uid==22 || uid==38 ) return data[12];
399 : // [0212] [2012] [0221] [2021] [1202] [1220] [2102] [2120]
400 9702448 : 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 6930320 : else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[14];
403 : // [1111]
404 5544256 : else if (uid==40 ) return data[15];
405 : // [1112] [1121] [1211] [2111]
406 5197740 : else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[16];
407 : // [1122] [2211]
408 3811676 : else if (uid==76 || uid==44 ) return data[17];
409 : // [1212] [2112] [1221] [2121]
410 3118644 : else if (uid==70 || uid==52 || uid==68 || uid==50 ) return data[18];
411 : // [1222] [2122] [2212] [2221]
412 1732580 : else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[19];
413 : // [2222]
414 346516 : else if (uid==80 ) return data[20];
415 0 : else Util::Abort(INFO,"Index out of range");
416 0 : return Set::Garbage;
417 : }
418 0 : void Print (std::ostream& os)
419 : {
420 0 : for (int i = 0; i < 21; i++)
421 0 : os << "i = " << i << " " << data[i] << std::endl;
422 :
423 0 : os.precision(4);
424 0 : for (int k = 0; k < 3; k++)
425 : {
426 0 : for (int i = -1; i < 4; i++)
427 : {
428 0 : for (int l = 0; l < 3; l++)
429 : {
430 0 : if (i==-1) os << "┌ ┐ ";
431 0 : else if (i == 3) os << "└ ┘ ";
432 : else
433 : {
434 0 : os << "│ ";
435 0 : for (int j = 0; j < 3; j++)
436 : {
437 0 : const Set::Scalar &val = (*this)(i,j,k,l);
438 0 : os << std::scientific << std::setw(11) << val ; //(fabs(val)>1E-10 ? val : 0);
439 0 : os << " ";
440 : }
441 0 : os << "│ ";
442 : }
443 : }
444 0 : os<<std::endl;
445 : }
446 : }
447 0 : }
448 :
449 : Set::Scalar Norm()
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 2750 : 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 :
489 : static Matrix4<3,Sym::MajorMinor> Increment()
490 : {
491 : Matrix4<3,Sym::MajorMinor> ret;
492 : for (int i = 0 ; i < 21; i++) ret.data[i] = (Set::Scalar)i;
493 : return ret;
494 : }
495 2 : static Matrix4<3,Sym::MajorMinor> Randomize()
496 : {
497 2 : Matrix4<3,Sym::MajorMinor> ret;
498 44 : for (int i = 0 ; i < 21; i++) ret.data[i] = (Util::Random());
499 2 : return ret;
500 : }
501 323 : static Matrix4<3,Sym::MajorMinor> Zero()
502 : {
503 323 : Matrix4<3,Sym::MajorMinor> ret;
504 7106 : for (int i = 0 ; i < 21; i++) ret.data[i] = 0.0;
505 323 : return ret;
506 : }
507 3906 : static Matrix4<3,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
508 : Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
509 : {
510 3906 : Set::Matrix4<3,Sym::MajorMinor> ret;
511 : Set::Scalar Ctmp[3][3][3][3];
512 :
513 15624 : for(int i = 0; i < 3; i++)
514 46872 : for(int j = 0; j < 3; j++)
515 140616 : for(int k = 0; k < 3; k++)
516 421848 : for(int l = 0; l < 3; l++)
517 : {
518 316386 : if(i == j && j == k && k == l) Ctmp[i][j][k][l] = C11;
519 304668 : else if (i==k && j==l) Ctmp[i][j][k][l] = C44;
520 281232 : else if (i==j && k==l) Ctmp[i][j][k][l] = C12;
521 257796 : else Ctmp[i][j][k][l] = 0.0;
522 : }
523 15624 : for(int p = 0; p < 3; p++)
524 46872 : for(int q = 0; q < 3; q++)
525 140616 : for(int s = 0; s < 3; s++)
526 421848 : for(int t = 0; t < 3; t++)
527 : {
528 316386 : ret(p,q,s,t) = 0.0;
529 1265544 : for(int i = 0; i < 3; i++)
530 3796632 : for(int j = 0; j < 3; j++)
531 11389896 : for(int k = 0; k < 3; k++)
532 34169688 : for(int l = 0; l < 3; l++)
533 51254532 : ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
534 : }
535 7812 : return ret;
536 : }
537 46 : static Matrix4<3,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
538 : Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
539 : {
540 46 : Eigen::Matrix3d R;
541 46 : R = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
542 92 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
543 92 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
544 46 : return Matrix4<3,Sym::MajorMinor>::Cubic(C11,C12,C44,R);
545 : }
546 3860 : static Matrix4<3,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
547 : Set::Quaternion q)
548 : {
549 3860 : Eigen::Matrix3d R = q.normalized().toRotationMatrix();
550 3860 : return Cubic(C11,C12,C44,R);
551 : }
552 :
553 23 : 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 : {
555 23 : Set::Matrix4<3,Sym::MajorMinor> ret;
556 23 : Set::Scalar Ctmp[3][3][3][3] = {{{{0.0}}}};
557 :
558 92 : for(int i = 0; i < 3; i++)
559 276 : for(int j = 0; j < 3; j++)
560 828 : for(int k = 0; k < 3; k++)
561 2484 : for(int l = 0; l < 3; l++)
562 : {
563 1863 : if(i==j && j==k && k==l)
564 : {
565 69 : if (i < 2) Ctmp[i][j][k][l] = C11;
566 23 : else Ctmp[i][j][k][l] = C33;
567 : }
568 1794 : else if(i==j && k==l)
569 : {
570 138 : if (i<2 && k<2 && i!=k) Ctmp[i][j][k][l] = C12;
571 92 : else if ((i<2 && k==2) || (i==2 && k<2)) Ctmp[i][j][k][l] = C13;
572 : }
573 1656 : else if((i==k && j==l) && (i!=j))
574 : {
575 138 : Ctmp[i][j][k][l] = C44;
576 : }
577 : else
578 : {
579 1518 : Ctmp[i][j][k][l] = 0.0;
580 : }
581 : }
582 :
583 : // C66 = 0.5*(C11 - C12)
584 23 : Ctmp[0][1][0][1] = 0.5 * (C11 - C12);
585 23 : Ctmp[1][0][1][0] = 0.5 * (C11 - C12);
586 :
587 92 : for(int p = 0; p < 3; p++)
588 276 : for(int q = 0; q < 3; q++)
589 828 : for(int s = 0; s < 3; s++)
590 2484 : for(int t = 0; t < 3; t++)
591 : {
592 1863 : ret(p,q,s,t) = 0.0;
593 7452 : for(int i = 0; i < 3; i++)
594 22356 : for(int j = 0; j < 3; j++)
595 67068 : for(int k = 0; k < 3; k++)
596 201204 : for(int l = 0; l < 3; l++)
597 301806 : ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
598 : }
599 46 : return ret;
600 : }
601 22 : 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)
602 : {
603 : Eigen::Quaterniond q =
604 22 : Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
605 22 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
606 44 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
607 22 : Eigen::Matrix3d R = q.toRotationMatrix();
608 22 : return Matrix4<3,Sym::MajorMinor>::Transverse(C11,C12,C13,C33,C44,R);
609 : }
610 : static Matrix4<3,Sym::MajorMinor> Trans(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Quaternion q)
611 : {
612 : Eigen::Matrix3d R = q.normalized().toRotationMatrix();
613 : return Transverse(C11,C12,C13,C33,C44,R);
614 : }
615 :
616 : friend bool operator == (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b);
617 : friend Matrix4<3,Sym::MajorMinor> operator + (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b);
618 : friend Matrix4<3,Sym::MajorMinor> operator - (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b);
619 : friend Matrix4<3,Sym::MajorMinor> operator * (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b);
620 : friend Matrix4<3,Sym::MajorMinor> operator * (Set::Scalar b, Matrix4<3,Sym::MajorMinor> a);
621 : friend Matrix4<3,Sym::MajorMinor> operator / (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b);
622 : friend Set::Vector operator * (Matrix4<3,Sym::MajorMinor> a, Set::Matrix3 b);
623 : friend Eigen::Matrix<amrex::Real,3,3> operator * (Matrix4<3,Sym::MajorMinor> a, Eigen::Matrix<amrex::Real,3,3> b);
624 : };
625 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
626 : bool operator == (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
627 : {
628 550 : for (int i = 0 ; i < 21 ; i++) if (a.data[i] != b.data[i]) return false;
629 25 : return true;
630 : }
631 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
632 : Matrix4<3,Sym::MajorMinor> operator + (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
633 : {
634 9800 : Matrix4<3,Sym::MajorMinor> ret;
635 215600 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] + b.data[i];
636 9800 : return ret;
637 : }
638 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
639 : Matrix4<3,Sym::MajorMinor> operator - (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
640 : {
641 294597 : Matrix4<3,Sym::MajorMinor> ret;
642 6481134 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] - b.data[i];
643 294597 : return ret;
644 : }
645 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
646 : Matrix4<3,Sym::MajorMinor> operator * (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b)
647 : {
648 301332 : Matrix4<3,Sym::MajorMinor> ret;
649 6629304 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] * b;
650 301332 : return ret;
651 : }
652 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
653 : Matrix4<3,Sym::MajorMinor> operator * (Set::Scalar b,Matrix4<3,Sym::MajorMinor> a)
654 : {
655 0 : return a*b;
656 : }
657 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
658 : Matrix4<3,Sym::MajorMinor> operator / (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b)
659 : {
660 295597 : Matrix4<3,Sym::MajorMinor> ret;
661 6503134 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] / b;
662 295597 : return ret;
663 : }
664 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
665 : Eigen::Matrix<amrex::Real,3,3> operator * (Matrix4<3,Sym::MajorMinor> a, Eigen::Matrix<amrex::Real,3,3> b)
666 : {
667 1413576 : Eigen::Matrix<amrex::Real,3,3> ret = Eigen::Matrix<amrex::Real,3,3>::Zero();
668 1413576 : 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 1413576 : 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 1413576 : 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 1413576 : 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 1413576 : 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 1413576 : 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 1413576 : ret(2,1) = ret(1,2);
676 1413576 : ret(0,2) = ret(2,0);
677 1413576 : ret(1,0) = ret(0,1);
678 :
679 1413576 : return ret;
680 : }
681 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
682 : Eigen::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 : }
692 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
693 : Eigen::Matrix<amrex::Real,3,1> operator * (Matrix4<3,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,3,3>,3> b);
694 :
695 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
696 : Eigen::Matrix<amrex::Real,2,1> operator * (Matrix4<3,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,2,2>,2> b);
697 :
698 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
699 : Set::Vector operator * (Matrix4<3,Sym::MajorMinor> a, Set::Matrix3 b)
700 : {
701 : // TODO: improve efficiency of this method
702 106599 : Set::Vector ret = Set::Vector::Zero();
703 :
704 213198 : ret(0) +=
705 106599 : a.data[ 0]*b(0,0,0) +
706 319797 : a.data[ 1]*(b(0,1,0) + b(1,0,0) + b(0,0,1)) +
707 319797 : a.data[ 2]*(b(0,2,0) + b(2,0,0) + b(0,0,2)) +
708 106599 : a.data[ 3]*b(1,1,0) +
709 213198 : a.data[ 4]*(b(1,2,0) + b(2,1,0)) +
710 106599 : a.data[ 5]*b(2,2,0) +
711 213198 : a.data[ 6]*(b(0,1,1) + b(1,0,1)) +
712 426396 : a.data[ 7]*(b(0,2,1) + b(2,0,1) + b(0,1,2) + b(1,0,2)) +
713 106599 : a.data[ 8]*b(1,1,1) +
714 213198 : a.data[ 9]*(b(1,2,1) + b(2,1,1)) +
715 106599 : a.data[ 10]*b(2,2,1) +
716 213198 : a.data[ 11]*(b(0,2,2) + b(2,0,2)) +
717 106599 : a.data[ 12]*b(1,1,2) +
718 213198 : a.data[ 13]*(b(1,2,2) + b(2,1,2)) +
719 213198 : a.data[ 14]*b(2,2,2);
720 :
721 213198 : ret(1) +=
722 106599 : a.data[ 1]*b(0,0,0) +
723 106599 : a.data[ 4]*b(0,0,2) +
724 106599 : a.data[ 3]*b(0,0,1) +
725 213198 : a.data[ 6]*(b(0,1,0) + b(1,0,0)) +
726 213198 : a.data[ 7]*(b(0,2,0) + b(2,0,0)) +
727 319797 : a.data[ 8]*(b(1,1,0) + b(0,1,1) + b(1,0,1)) +
728 426396 : a.data[ 9]*(b(1,2,0) + b(2,1,0) + b(0,1,2) + b(1,0,2)) +
729 106599 : a.data[ 10]*b(2,2,0) +
730 213198 : a.data[ 12]*(b(0,2,1) + b(2,0,1)) +
731 213198 : a.data[ 13]*(b(0,2,2) + b(2,0,2)) +
732 106599 : a.data[ 15]*b(1,1,1) +
733 319797 : a.data[ 16]*(b(1,2,1) + b(2,1,1) + b(1,1,2)) +
734 106599 : a.data[ 17]*b(2,2,1) +
735 213198 : a.data[ 18]*(b(1,2,2) + b(2,1,2)) +
736 213198 : a.data[ 19]*b(2,2,2);
737 :
738 213198 : ret(2) +=
739 106599 : a.data[ 2]*b(0,0,0) +
740 106599 : a.data[ 4]*b(0,0,1) +
741 106599 : a.data[ 5]*b(0,0,2) +
742 213198 : a.data[ 7]*(b(0,1,0) + b(1,0,0)) +
743 213198 : a.data[ 9]*(b(0,1,1) + b(1,0,1)) +
744 213198 : a.data[ 10]*(b(0,1,2) + b(1,0,2)) +
745 213198 : a.data[ 11]*(b(0,2,0) + b(2,0,0)) +
746 106599 : a.data[ 12]*b(1,1,0) +
747 426396 : a.data[ 13]*(b(1,2,0) + b(2,1,0) + b(0,2,1) + b(2,0,1)) +
748 319797 : a.data[ 14]*(b(2,2,0) + b(0,2,2) + b(2,0,2)) +
749 106599 : a.data[ 16]*b(1,1,1) +
750 106599 : a.data[ 17]*b(1,1,2) +
751 213198 : a.data[ 18]*(b(1,2,1) + b(2,1,1)) +
752 319797 : a.data[ 19]*(b(2,2,1) + b(1,2,2) + b(2,1,2)) +
753 213198 : a.data[ 20]*b(2,2,2);
754 :
755 106599 : return ret;
756 : }
757 :
758 : }
759 : #endif
|
