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