Line data Source code
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 :
10 : namespace 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 : ///
16 : template<>
17 : class Matrix4<2,Sym::MajorMinor>
18 : {
19 : Scalar data[6] = {NAN,NAN,NAN,NAN,NAN,NAN};
20 :
21 : public:
22 1373 : AMREX_GPU_HOST_DEVICE Matrix4() {};
23 : //Matrix4(const Matrix4<2,Sym::MajorMinor> &in)
24 : //{
25 : // for (int i = 0; i < 6; i++) data[i] = in.data[i];
26 : //}
27 : AMREX_FORCE_INLINE
28 : const Scalar & operator () (const int i, const int j, const int k, const int l) const
29 : {
30 1600 : int uid = i + 2*j + 4*k + 8*l;
31 1600 : if (uid==0 ) return data[0]; // [0000]
32 1500 : else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1]; // [0001] [0010] [0100] [1000]
33 1100 : else if (uid==12 || uid==3 ) return data[2]; // [0011] [1100]
34 900 : else if (uid==10 || uid==6 || uid==9 || uid==5 ) return data[3]; // [0101] [1001] [0110] [1010]
35 500 : else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[4]; // [0111] [1011] [1101] [1110]
36 100 : else if (uid==15 ) return data[5]; // [1111]
37 0 : else Util::Abort(INFO,"Index out of range");
38 0 : return Set::Garbage;
39 : }
40 : AMREX_FORCE_INLINE
41 : Scalar & operator () (const int i, const int j, const int k, const int l)
42 : {
43 1537696 : int uid = i + 2*j + 4*k + 8*l;
44 190723 : if (uid==0 ) return data[0]; // [0000]
45 1461390 : else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1]; // [0001] [0010] [0100] [1000]
46 1071686 : else if (uid==12 || uid==3 ) return data[2]; // [0011] [1100]
47 876834 : else if (uid==10 || uid==6 || uid==9 || uid==5 ) return data[3]; // [0101] [1001] [0110] [1010]
48 487130 : else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[4]; // [0111] [1011] [1101] [1110]
49 97426 : else if (uid==15 ) return data[5]; // [1111]
50 0 : else Util::Abort(INFO,"Index out of range");
51 0 : 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 :
81 : Set::Scalar Norm()
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 0 : 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 :
103 : static Matrix4<2,Sym::MajorMinor> Increment()
104 : {
105 : Matrix4<2,Sym::MajorMinor> ret;
106 : for (int i = 0 ; i < 6; i++) ret.data[i] = (Set::Scalar)i;
107 : return ret;
108 : }
109 1 : void Randomize()
110 : {
111 7 : for (int i = 0 ; i < 6; i++) data[i] = (Util::Random());
112 1 : }
113 1 : static Matrix4<2,Sym::MajorMinor> Random()
114 : {
115 1 : Matrix4<2,Sym::MajorMinor> ret;
116 1 : ret.Randomize();
117 1 : return ret;
118 : }
119 94 : static Matrix4<2,Sym::MajorMinor> Zero()
120 : {
121 94 : Matrix4<2,Sym::MajorMinor> ret;
122 658 : for (int i = 0 ; i < 6; i++) ret.data[i] = 0.0;
123 94 : return ret;
124 : }
125 1104 : static Matrix4<2,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
126 : Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
127 : {
128 1104 : Set::Matrix4<2,Sym::MajorMinor> ret;
129 : Set::Scalar Ctmp[3][3][3][3];
130 :
131 4416 : for(int i = 0; i < 3; i++)
132 13248 : for(int j = 0; j < 3; j++)
133 39744 : for(int k = 0; k < 3; k++)
134 119232 : for(int l = 0; l < 3; l++)
135 : {
136 89424 : if(i == j && j == k && k == l) Ctmp[i][j][k][l] = C11;
137 86112 : else if (i==k && j==l) Ctmp[i][j][k][l] = C44;
138 79488 : else if (i==j && k==l) Ctmp[i][j][k][l] = C12;
139 72864 : else Ctmp[i][j][k][l] = 0.0;
140 : }
141 3312 : for(int p = 0; p < 2; p++)
142 6624 : for(int q = 0; q < 2; q++)
143 13248 : for(int s = 0; s < 2; s++)
144 26496 : for(int t = 0; t < 2; t++)
145 : {
146 17664 : ret(p,q,s,t) = 0.0;
147 70656 : for(int i = 0; i < 3; i++)
148 211968 : for(int j = 0; j < 3; j++)
149 635904 : for(int k = 0; k < 3; k++)
150 1907712 : for(int l = 0; l < 3; l++)
151 2861568 : ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
152 : }
153 2208 : return ret;
154 : }
155 44 : static Matrix4<2,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
156 : Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
157 : {
158 44 : Eigen::Matrix3d R;
159 44 : R = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
160 88 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
161 88 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
162 44 : return Matrix4<2,Sym::MajorMinor>::Cubic(C11,C12,C44,R);
163 : }
164 1060 : static Matrix4<2,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
165 : Set::Quaternion q)
166 : {
167 1060 : Eigen::Matrix3d R = q.normalized().toRotationMatrix();
168 1060 : return Cubic(C11,C12,C44,R);
169 : }
170 :
171 :
172 22 : 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 : {
174 22 : Set::Matrix4<2,Sym::MajorMinor> ret;
175 22 : Set::Scalar Ctmp[3][3][3][3] = {{{{0.0}}}};
176 :
177 88 : for(int i = 0; i < 3; i++)
178 264 : for(int j = 0; j < 3; j++)
179 792 : for(int k = 0; k < 3; k++)
180 2376 : for(int l = 0; l < 3; l++)
181 : {
182 1782 : if(i==j && j==k && k==l)
183 : {
184 66 : if (i < 2) Ctmp[i][j][k][l] = C11;
185 22 : else Ctmp[i][j][k][l] = C33;
186 : }
187 1716 : else if(i==j && k==l)
188 : {
189 132 : if (i<2 && k<2 && i!=k) Ctmp[i][j][k][l] = C12;
190 88 : else if ((i<2 && k==2) || (i==2 && k<2)) Ctmp[i][j][k][l] = C13;
191 : }
192 1584 : else if((i==k && j==l) && (i!=j))
193 : {
194 132 : Ctmp[i][j][k][l] = C44;
195 : }
196 : else
197 : {
198 1452 : Ctmp[i][j][k][l] = 0.0;
199 : }
200 : }
201 :
202 : // C66 = 0.5*(C11 - C12)
203 22 : Ctmp[0][1][0][1] = 0.5 * (C11 - C12);
204 22 : Ctmp[1][0][1][0] = 0.5 * (C11 - C12);
205 :
206 66 : for(int p = 0; p < 2; p++)
207 132 : for(int q = 0; q < 2; q++)
208 264 : for(int s = 0; s < 2; s++)
209 528 : for(int t = 0; t < 2; t++)
210 : {
211 352 : ret(p,q,s,t) = 0.0;
212 1408 : for(int i = 0; i < 3; i++)
213 4224 : for(int j = 0; j < 3; j++)
214 12672 : for(int k = 0; k < 3; k++)
215 38016 : for(int l = 0; l < 3; l++)
216 57024 : ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
217 : }
218 44 : return ret;
219 : }
220 21 : 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)
221 : {
222 : Eigen::Quaterniond q =
223 21 : Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
224 21 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
225 42 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
226 21 : Eigen::Matrix3d R = q.toRotationMatrix();
227 21 : return Matrix4<2,Sym::MajorMinor>::Transverse(C11,C12,C13,C33,C44,R);
228 : }
229 : static Matrix4<2,Sym::MajorMinor> Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Quaternion q)
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);
237 : friend bool operator == (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b);
238 : friend Matrix4<2,Sym::MajorMinor> operator + (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b);
239 : friend Matrix4<2,Sym::MajorMinor> operator - (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b);
240 : friend Matrix4<2,Sym::MajorMinor> operator * (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b);
241 : friend Matrix4<2,Sym::MajorMinor> operator * (Set::Scalar b, Matrix4<2,Sym::MajorMinor> a);
242 : friend Matrix4<2,Sym::MajorMinor> operator / (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b);
243 : friend Set::Vector operator * (Matrix4<2,Sym::MajorMinor> a, Set::Matrix3 b);
244 : };
245 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
246 : bool operator == (Matrix4<2,Sym::MajorMinor> a, Matrix4<2,Sym::MajorMinor> b)
247 : {
248 175 : for (int i = 0 ; i < 6 ; i++) if (a.data[i] != b.data[i]) return false;
249 25 : return true;
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, Matrix4<2,Sym::MajorMinor> b)
260 : {
261 0 : Matrix4<2,Sym::MajorMinor> ret;
262 0 : for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] - b.data[i];
263 0 : return ret;
264 : }
265 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
266 : Matrix4<2,Sym::MajorMinor> operator * (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b)
267 : {
268 10 : Matrix4<2,Sym::MajorMinor> ret;
269 70 : for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] * b;
270 10 : return ret;
271 : }
272 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
273 : Matrix4<2,Sym::MajorMinor> operator * (Set::Scalar b,Matrix4<2,Sym::MajorMinor> a)
274 : {
275 0 : return a*b;
276 : }
277 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
278 : Matrix4<2,Sym::MajorMinor> operator / (Matrix4<2,Sym::MajorMinor> a, Set::Scalar b)
279 : {
280 0 : Matrix4<2,Sym::MajorMinor> ret;
281 0 : for (int i = 0 ; i < 6 ; i++) ret.data[i] = a.data[i] / b;
282 0 : return ret;
283 : }
284 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
285 : Eigen::Matrix<Set::Scalar,2,2> operator * (Matrix4<2,Sym::MajorMinor> a, Eigen::Matrix<Set::Scalar,2,2> b)
286 : {
287 2430 : Eigen::Matrix<Set::Scalar,2,2> ret = Eigen::Matrix<Set::Scalar,2,2>::Zero();
288 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);
289 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);
290 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);
291 2430 : ret(1,0) = ret(0,1);
292 2430 : return ret;
293 : }
294 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
295 : Eigen::Matrix<amrex::Real,2,1> operator * (Matrix4<2,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,2,2>,2> b);
296 :
297 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
298 : Set::Vector operator * (Matrix4<2,Sym::MajorMinor> a, Set::Matrix3 b)
299 : {
300 0 : Set::Vector ret = Set::Vector::Zero();
301 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);
302 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);
303 0 : 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 : ///
318 : template<>
319 : class Matrix4<3,Sym::MajorMinor>
320 : {
321 : Scalar data[21] = { NAN,NAN,NAN,NAN,NAN,NAN,NAN,
322 : NAN,NAN,NAN,NAN,NAN,NAN,NAN,
323 : NAN,NAN,NAN,NAN,NAN,NAN,NAN };
324 :
325 : public:
326 1425097 : AMREX_GPU_HOST_DEVICE Matrix4() {};
327 : AMREX_FORCE_INLINE
328 : const Scalar & operator () (const int i, const int j, const int k, const int l) const
329 : {
330 9475 : int uid = i + 3*j + 9*k + 27*l;
331 : // [0000]
332 9475 : if (uid==0 ) return data[0];
333 : // [0001] [0010] [0100] [1000]
334 8000 : else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
335 : // [0002] [0020] [0200] [2000]
336 7600 : else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
337 : // [0011] [1100]
338 7200 : else if (uid==36 || uid==4 ) return data[3];
339 : // [0012] [0021] [1200] [2100]
340 7000 : else if (uid==63 || uid==45 || uid==7 || uid==5 ) return data[4];
341 : // [0022] [2200]
342 6600 : else if (uid==72 || uid==8 ) return data[5];
343 : // [0101] [0110] [1001] [1010]
344 6400 : else if (uid==30 || uid==12 || uid==28 || uid==10 ) return data[6];
345 : // [0102] [0120] [1002] [1020] [0201] [2001] [0210] [2010]
346 6000 : 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 5200 : else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[8];
349 : // [0112] [1012] [0121] [1021] [1201] [1210] [2101] [2110]
350 4800 : 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 4000 : else if (uid==75 || uid==73 || uid==35 || uid==17 ) return data[10];
353 : // [0202] [2002] [0220] [2020]
354 3600 : else if (uid==60 || uid==24 || uid==56 || uid==20 ) return data[11];
355 : // [0211] [2011] [1102] [1120]
356 3200 : else if (uid==42 || uid==58 || uid==22 || uid==38 ) return data[12];
357 : // [0212] [2012] [0221] [2021] [1202] [1220] [2102] [2120]
358 2800 : 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 2000 : else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[14];
361 : // [1111]
362 1600 : else if (uid==40 ) return data[15];
363 : // [1112] [1121] [1211] [2111]
364 1500 : else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[16];
365 : // [1122] [2211]
366 1100 : else if (uid==76 || uid==44 ) return data[17];
367 : // [1212] [2112] [1221] [2121]
368 900 : else if (uid==70 || uid==52 || uid==68 || uid==50 ) return data[18];
369 : // [1222] [2122] [2212] [2221]
370 500 : else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[19];
371 : // [2222]
372 100 : else if (uid==80 ) return data[20];
373 0 : else Util::Abort(INFO,"Index out of range");
374 0 : return Set::Garbage;
375 : }
376 : AMREX_FORCE_INLINE
377 : Scalar & operator () (const int i, const int j, const int k, const int l)
378 : {
379 28034100 : int uid = i + 3*j + 9*k + 27*l;
380 : // [0000]
381 2593617 : if (uid==0 ) return data[0];
382 : // [0001] [0010] [0100] [1000]
383 27688000 : else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
384 : // [0002] [0020] [0200] [2000]
385 26303600 : else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
386 : // [0011] [1100]
387 24919200 : else if (uid==36 || uid==4 ) return data[3];
388 : // [0012] [0021] [1200] [2100]
389 24227000 : else if (uid==63 || uid==45 || uid==7 || uid==5 ) return data[4];
390 : // [0022] [2200]
391 22842600 : else if (uid==72 || uid==8 ) return data[5];
392 : // [0101] [0110] [1001] [1010]
393 22150400 : else if (uid==30 || uid==12 || uid==28 || uid==10 ) return data[6];
394 : // [0102] [0120] [1002] [1020] [0201] [2001] [0210] [2010]
395 20766000 : 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 17997200 : else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[8];
398 : // [0112] [1012] [0121] [1021] [1201] [1210] [2101] [2110]
399 16612800 : 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 13844000 : else if (uid==75 || uid==73 || uid==35 || uid==17 ) return data[10];
402 : // [0202] [2002] [0220] [2020]
403 12459600 : else if (uid==60 || uid==24 || uid==56 || uid==20 ) return data[11];
404 : // [0211] [2011] [1102] [1120]
405 11075200 : else if (uid==42 || uid==58 || uid==22 || uid==38 ) return data[12];
406 : // [0212] [2012] [0221] [2021] [1202] [1220] [2102] [2120]
407 9690800 : 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 6922000 : else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[14];
410 : // [1111]
411 5537600 : else if (uid==40 ) return data[15];
412 : // [1112] [1121] [1211] [2111]
413 5191500 : else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[16];
414 : // [1122] [2211]
415 3807100 : else if (uid==76 || uid==44 ) return data[17];
416 : // [1212] [2112] [1221] [2121]
417 3114900 : else if (uid==70 || uid==52 || uid==68 || uid==50 ) return data[18];
418 : // [1222] [2122] [2212] [2221]
419 1730500 : else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[19];
420 : // [2222]
421 346100 : else if (uid==80 ) return data[20];
422 0 : else Util::Abort(INFO,"Index out of range");
423 0 : 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 :
456 : Set::Scalar Norm()
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 2750 : 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 :
496 : static Matrix4<3,Sym::MajorMinor> Increment()
497 : {
498 : Matrix4<3,Sym::MajorMinor> ret;
499 : for (int i = 0 ; i < 21; i++) ret.data[i] = (Set::Scalar)i;
500 : return ret;
501 : }
502 1 : void Randomize()
503 : {
504 22 : for (int i = 0 ; i < 21; i++) data[i] = (Util::Random());
505 1 : }
506 1 : static Matrix4<3,Sym::MajorMinor> Random()
507 : {
508 1 : Matrix4<3,Sym::MajorMinor> ret;
509 1 : ret.Randomize();
510 1 : return ret;
511 : }
512 321 : static Matrix4<3,Sym::MajorMinor> Zero()
513 : {
514 321 : Matrix4<3,Sym::MajorMinor> ret;
515 7062 : for (int i = 0 ; i < 21; i++) ret.data[i] = 0.0;
516 321 : return ret;
517 : }
518 3904 : static Matrix4<3,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
519 : Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
520 : {
521 3904 : Set::Matrix4<3,Sym::MajorMinor> ret;
522 : Set::Scalar Ctmp[3][3][3][3];
523 :
524 15616 : for(int i = 0; i < 3; i++)
525 46848 : for(int j = 0; j < 3; j++)
526 140544 : for(int k = 0; k < 3; k++)
527 421632 : for(int l = 0; l < 3; l++)
528 : {
529 316224 : if(i == j && j == k && k == l) Ctmp[i][j][k][l] = C11;
530 304512 : else if (i==k && j==l) Ctmp[i][j][k][l] = C44;
531 281088 : else if (i==j && k==l) Ctmp[i][j][k][l] = C12;
532 257664 : else Ctmp[i][j][k][l] = 0.0;
533 : }
534 15616 : for(int p = 0; p < 3; p++)
535 46848 : for(int q = 0; q < 3; q++)
536 140544 : for(int s = 0; s < 3; s++)
537 421632 : for(int t = 0; t < 3; t++)
538 : {
539 316224 : ret(p,q,s,t) = 0.0;
540 1264896 : for(int i = 0; i < 3; i++)
541 3794688 : for(int j = 0; j < 3; j++)
542 11384064 : for(int k = 0; k < 3; k++)
543 34152192 : for(int l = 0; l < 3; l++)
544 51228288 : ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
545 : }
546 7808 : return ret;
547 : }
548 44 : static Matrix4<3,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
549 : Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
550 : {
551 44 : Eigen::Matrix3d R;
552 44 : R = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
553 88 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
554 88 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
555 44 : return Matrix4<3,Sym::MajorMinor>::Cubic(C11,C12,C44,R);
556 : }
557 3860 : static Matrix4<3,Sym::MajorMinor> Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44,
558 : Set::Quaternion q)
559 : {
560 3860 : Eigen::Matrix3d R = q.normalized().toRotationMatrix();
561 3860 : return Cubic(C11,C12,C44,R);
562 : }
563 :
564 22 : 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 : {
566 22 : Set::Matrix4<3,Sym::MajorMinor> ret;
567 22 : Set::Scalar Ctmp[3][3][3][3] = {{{{0.0}}}};
568 :
569 88 : for(int i = 0; i < 3; i++)
570 264 : for(int j = 0; j < 3; j++)
571 792 : for(int k = 0; k < 3; k++)
572 2376 : for(int l = 0; l < 3; l++)
573 : {
574 1782 : if(i==j && j==k && k==l)
575 : {
576 66 : if (i < 2) Ctmp[i][j][k][l] = C11;
577 22 : else Ctmp[i][j][k][l] = C33;
578 : }
579 1716 : else if(i==j && k==l)
580 : {
581 132 : if (i<2 && k<2 && i!=k) Ctmp[i][j][k][l] = C12;
582 88 : else if ((i<2 && k==2) || (i==2 && k<2)) Ctmp[i][j][k][l] = C13;
583 : }
584 1584 : else if((i==k && j==l) && (i!=j))
585 : {
586 132 : Ctmp[i][j][k][l] = C44;
587 : }
588 : else
589 : {
590 1452 : Ctmp[i][j][k][l] = 0.0;
591 : }
592 : }
593 :
594 : // C66 = 0.5*(C11 - C12)
595 22 : Ctmp[0][1][0][1] = 0.5 * (C11 - C12);
596 22 : Ctmp[1][0][1][0] = 0.5 * (C11 - C12);
597 :
598 88 : for(int p = 0; p < 3; p++)
599 264 : for(int q = 0; q < 3; q++)
600 792 : for(int s = 0; s < 3; s++)
601 2376 : for(int t = 0; t < 3; t++)
602 : {
603 1782 : ret(p,q,s,t) = 0.0;
604 7128 : for(int i = 0; i < 3; i++)
605 21384 : for(int j = 0; j < 3; j++)
606 64152 : for(int k = 0; k < 3; k++)
607 192456 : for(int l = 0; l < 3; l++)
608 288684 : ret(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
609 : }
610 44 : return ret;
611 : }
612 21 : 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)
613 : {
614 : Eigen::Quaterniond q =
615 21 : Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
616 21 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
617 42 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
618 21 : Eigen::Matrix3d R = q.toRotationMatrix();
619 21 : return Matrix4<3,Sym::MajorMinor>::Transverse(C11,C12,C13,C33,C44,R);
620 : }
621 : static Matrix4<3,Sym::MajorMinor> Trans(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Quaternion q)
622 : {
623 : Eigen::Matrix3d R = q.normalized().toRotationMatrix();
624 : return Transverse(C11,C12,C13,C33,C44,R);
625 : }
626 :
627 : friend bool operator == (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b);
628 : friend Matrix4<3,Sym::MajorMinor> operator + (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b);
629 : friend Matrix4<3,Sym::MajorMinor> operator - (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b);
630 : friend Matrix4<3,Sym::MajorMinor> operator * (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b);
631 : friend Matrix4<3,Sym::MajorMinor> operator * (Set::Scalar b, Matrix4<3,Sym::MajorMinor> a);
632 : friend Matrix4<3,Sym::MajorMinor> operator / (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b);
633 : friend Set::Vector operator * (Matrix4<3,Sym::MajorMinor> a, Set::Matrix3 b);
634 : friend Eigen::Matrix<amrex::Real,3,3> operator * (Matrix4<3,Sym::MajorMinor> a, Eigen::Matrix<amrex::Real,3,3> b);
635 : };
636 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
637 : bool operator == (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
638 : {
639 550 : for (int i = 0 ; i < 21 ; i++) if (a.data[i] != b.data[i]) return false;
640 25 : return true;
641 : }
642 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
643 : Matrix4<3,Sym::MajorMinor> operator + (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
644 : {
645 9800 : Matrix4<3,Sym::MajorMinor> ret;
646 215600 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] + b.data[i];
647 9800 : return ret;
648 : }
649 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
650 : Matrix4<3,Sym::MajorMinor> operator - (Matrix4<3,Sym::MajorMinor> a, Matrix4<3,Sym::MajorMinor> b)
651 : {
652 294597 : Matrix4<3,Sym::MajorMinor> ret;
653 6481134 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] - b.data[i];
654 294597 : return ret;
655 : }
656 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
657 : Matrix4<3,Sym::MajorMinor> operator * (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b)
658 : {
659 301332 : Matrix4<3,Sym::MajorMinor> ret;
660 6629304 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] * b;
661 301332 : return ret;
662 : }
663 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
664 : Matrix4<3,Sym::MajorMinor> operator * (Set::Scalar b,Matrix4<3,Sym::MajorMinor> a)
665 : {
666 0 : return a*b;
667 : }
668 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
669 : Matrix4<3,Sym::MajorMinor> operator / (Matrix4<3,Sym::MajorMinor> a, Set::Scalar b)
670 : {
671 295597 : Matrix4<3,Sym::MajorMinor> ret;
672 6503134 : for (int i = 0 ; i < 21 ; i++) ret.data[i] = a.data[i] / b;
673 295597 : return ret;
674 : }
675 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
676 : Eigen::Matrix<amrex::Real,3,3> operator * (Matrix4<3,Sym::MajorMinor> a, Eigen::Matrix<amrex::Real,3,3> b)
677 : {
678 1413576 : Eigen::Matrix<amrex::Real,3,3> ret = Eigen::Matrix<amrex::Real,3,3>::Zero();
679 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);
680 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);
681 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);
682 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);
683 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);
684 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);
685 :
686 1413576 : ret(2,1) = ret(1,2);
687 1413576 : ret(0,2) = ret(2,0);
688 1413576 : ret(1,0) = ret(0,1);
689 :
690 1413576 : return ret;
691 : }
692 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
693 : Eigen::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 : }
703 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
704 : Eigen::Matrix<amrex::Real,3,1> operator * (Matrix4<3,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,3,3>,3> b);
705 :
706 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
707 : Eigen::Matrix<amrex::Real,2,1> operator * (Matrix4<3,Sym::MajorMinor> a, std::array<Eigen::Matrix<amrex::Real,2,2>,2> b);
708 :
709 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
710 : Set::Vector operator * (Matrix4<3,Sym::MajorMinor> a, Set::Matrix3 b)
711 : {
712 : // TODO: improve efficiency of this method
713 106599 : Set::Vector ret = Set::Vector::Zero();
714 :
715 213198 : ret(0) +=
716 106599 : a.data[ 0]*b(0,0,0) +
717 319797 : a.data[ 1]*(b(0,1,0) + b(1,0,0) + b(0,0,1)) +
718 319797 : a.data[ 2]*(b(0,2,0) + b(2,0,0) + b(0,0,2)) +
719 106599 : a.data[ 3]*b(1,1,0) +
720 213198 : a.data[ 4]*(b(1,2,0) + b(2,1,0)) +
721 106599 : a.data[ 5]*b(2,2,0) +
722 213198 : a.data[ 6]*(b(0,1,1) + b(1,0,1)) +
723 426396 : a.data[ 7]*(b(0,2,1) + b(2,0,1) + b(0,1,2) + b(1,0,2)) +
724 106599 : a.data[ 8]*b(1,1,1) +
725 213198 : a.data[ 9]*(b(1,2,1) + b(2,1,1)) +
726 106599 : a.data[ 10]*b(2,2,1) +
727 213198 : a.data[ 11]*(b(0,2,2) + b(2,0,2)) +
728 106599 : a.data[ 12]*b(1,1,2) +
729 213198 : a.data[ 13]*(b(1,2,2) + b(2,1,2)) +
730 213198 : a.data[ 14]*b(2,2,2);
731 :
732 213198 : ret(1) +=
733 106599 : a.data[ 1]*b(0,0,0) +
734 106599 : a.data[ 4]*b(0,0,2) +
735 106599 : a.data[ 3]*b(0,0,1) +
736 213198 : a.data[ 6]*(b(0,1,0) + b(1,0,0)) +
737 213198 : a.data[ 7]*(b(0,2,0) + b(2,0,0)) +
738 319797 : a.data[ 8]*(b(1,1,0) + b(0,1,1) + b(1,0,1)) +
739 426396 : a.data[ 9]*(b(1,2,0) + b(2,1,0) + b(0,1,2) + b(1,0,2)) +
740 106599 : a.data[ 10]*b(2,2,0) +
741 213198 : a.data[ 12]*(b(0,2,1) + b(2,0,1)) +
742 213198 : a.data[ 13]*(b(0,2,2) + b(2,0,2)) +
743 106599 : a.data[ 15]*b(1,1,1) +
744 319797 : a.data[ 16]*(b(1,2,1) + b(2,1,1) + b(1,1,2)) +
745 106599 : a.data[ 17]*b(2,2,1) +
746 213198 : a.data[ 18]*(b(1,2,2) + b(2,1,2)) +
747 213198 : a.data[ 19]*b(2,2,2);
748 :
749 213198 : ret(2) +=
750 106599 : a.data[ 2]*b(0,0,0) +
751 106599 : a.data[ 4]*b(0,0,1) +
752 106599 : a.data[ 5]*b(0,0,2) +
753 213198 : a.data[ 7]*(b(0,1,0) + b(1,0,0)) +
754 213198 : a.data[ 9]*(b(0,1,1) + b(1,0,1)) +
755 213198 : a.data[ 10]*(b(0,1,2) + b(1,0,2)) +
756 213198 : a.data[ 11]*(b(0,2,0) + b(2,0,0)) +
757 106599 : a.data[ 12]*b(1,1,0) +
758 426396 : a.data[ 13]*(b(1,2,0) + b(2,1,0) + b(0,2,1) + b(2,0,1)) +
759 319797 : a.data[ 14]*(b(2,2,0) + b(0,2,2) + b(2,0,2)) +
760 106599 : a.data[ 16]*b(1,1,1) +
761 106599 : a.data[ 17]*b(1,1,2) +
762 213198 : a.data[ 18]*(b(1,2,1) + b(2,1,1)) +
763 319797 : a.data[ 19]*(b(2,2,1) + b(1,2,2) + b(2,1,2)) +
764 213198 : a.data[ 20]*b(2,2,2);
765 :
766 106599 : return ret;
767 : }
768 :
769 : }
770 : #endif
|
