Alamo
Matrix4_Full.H
Go to the documentation of this file.
1 #ifndef SET_MATRIX4_FULL_H
2 #define SET_MATRIX4_FULL_H
3 #include "Util/Util.H"
4 #include "Base.H"
5 namespace Set
6 {
7 
8 /// \brief Data structure for a symmetrix 4th order 3D tensor
9 ///
10 /// Let the tensor
11 /// \f$\mathbb{C}\in\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{R}^3\f$
12 /// be fully symmetrix such that \f$\mathbb{C}_{ijkl}=\mathbb{C}_{\sigma(i,j,k,l)\f$ where
13 /// \f$\sigma\f$ is any permutation.
14 /// Then there are only 5 (2D) or 15(3D) unique elements (rather than 81).
15 ///
16 /// This object acts like a 4D array such that `C(i,j,k,l)` returns the corresponding
17 /// element, but symmetry is always obeyed. This allows the user code to be much prettier.
18 ///
19 
20 template<>
21 class Matrix4<2,Sym::Full>
22 {
23  Scalar data[5] = {NAN,NAN,NAN,NAN,NAN};
24 public:
25  AMREX_GPU_HOST_DEVICE Matrix4() {};
26  AMREX_FORCE_INLINE
27  AMREX_GPU_HOST_DEVICE
28  Scalar & operator () (const int i, const int j, const int k, const int l)
29  {
30  int uid = i + 2*j + 4*k + 8*l;
31  // [0, 0, 0, 0]
32  if (uid==0 ) return data[0];
33  // [0, 0, 0, 1]
34  else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1];
35  // [0, 0, 1, 1]
36  else if (uid==12 || uid==10 || uid==6 || uid==9 || uid==5 || uid==3 ) return data[2];
37  // [0, 1, 1, 1]
38  else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[3];
39  // [1, 1, 1, 1]
40  else if (uid==15 ) return data[4];
41  else Util::Abort(INFO,"Index out of range");
42  return data[0]; // this code is unreachable
43  }
44  Scalar operator () (const int i, const int j, const int k, const int l) const
45  {
46  int uid = i + 2*j + 4*k + 8*l;
47  // [0, 0, 0, 0]
48  if (uid==0 ) return data[0];
49  // [0, 0, 0, 1]
50  else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1];
51  // [0, 0, 1, 1]
52  else if (uid==12 || uid==10 || uid==6 || uid==9 || uid==5 || uid==3 ) return data[2];
53  // [0, 1, 1, 1]
54  else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[3];
55  // [1, 1, 1, 1]
56  else if (uid==15 ) return data[4];
57  else Util::Abort(INFO,"Index out of range");
58  return data[0]; // this code is unreachable
59  }
60  static Matrix4<2,Sym::Full> Randomize()
61  {
62  Matrix4<2,Sym::Full> ret;
63  for (int i = 0 ; i < 5; i++) ret.data[i] = Util::Random();
64  return ret;
65  }
66  void Print (std::ostream& /* os */)
67  {
68  Util::Abort(INFO,"Not yet implemented");
69  }
70  static Matrix4<2,Sym::Full> Zero()
71  {
72  Matrix4<2,Sym::Full> zero;
73  zero.data[0] = 0.0;
74  zero.data[1] = 0.0;
75  zero.data[2] = 0.0;
76  zero.data[3] = 0.0;
77  zero.data[4] = 0.0;
78  return zero;
79  }
80  bool contains_nan() const
81  {
82  if (std::isnan(data[0])) return true;
83  if (std::isnan(data[1])) return true;
84  if (std::isnan(data[2])) return true;
85  if (std::isnan(data[3])) return true;
86  if (std::isnan(data[4])) return true;
87  return false;
88  }
89 };
90 
91 template<>
92 class Matrix4<3,Sym::Full>
93 {
94  Scalar data[15] = { NAN,NAN,NAN,NAN,NAN,
95  NAN,NAN,NAN,NAN,NAN,
96  NAN,NAN,NAN,NAN,NAN};
97 
98 public:
99  AMREX_GPU_HOST_DEVICE Matrix4() {};
100  AMREX_FORCE_INLINE
101  AMREX_GPU_HOST_DEVICE
102  Scalar & operator () (const int i, const int j, const int k, const int l)
103  {
104  int uid = i + 3*j + 9*k + 27*l;
105  // [0, 0, 0, 0]
106  if (uid==0 ) return data[0];
107  // [0, 0, 0, 1]
108  else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
109  // [0, 0, 0, 2]
110  else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
111  // [0, 0, 1, 1]
112  else if (uid==36 || uid==30 || uid==12 || uid==28 || uid==10 || uid==4 ) return data[3];
113  // [0, 0, 1, 2]
114  else if (uid==63 || uid==45 || uid==57 || uid==21 || uid==33 || uid==15 || uid==55 || uid==19 || uid==7 || uid==29 || uid==11 || uid==5 ) return data[4];
115  // [0, 0, 2, 2]
116  else if (uid==72 || uid==60 || uid==24 || uid==56 || uid==20 || uid==8 ) return data[5];
117  // [0, 1, 1, 1]
118  else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[6];
119  // [0, 1, 1, 2]
120  else if (uid==66 || uid==48 || uid==42 || uid==64 || uid==46 || uid==58 || uid==22 || uid==34 || uid==16 || uid==38 || uid==32 || uid==14 ) return data[7];
121  // [0, 1, 2, 2]
122  else if (uid==75 || uid==69 || uid==51 || uid==73 || uid==61 || uid==25 || uid==65 || uid==47 || uid==59 || uid==23 || uid==35 || uid==17 ) return data[8];
123  // [0, 2, 2, 2]
124  else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[9];
125  // [1, 1, 1, 1]
126  else if (uid==40 ) return data[10];
127  // [1, 1, 1, 2]
128  else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[11];
129  // [1, 1, 2, 2]
130  else if (uid==76 || uid==70 || uid==52 || uid==68 || uid==50 || uid==44 ) return data[12];
131  // [1, 2, 2, 2]
132  else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[13];
133  // [2, 2, 2, 2]
134  else if (uid==80 ) return data[14];
135  else Util::Abort(INFO,"uid not in range (uid=",uid,")");
136  return data[0]; // this return statement is unreachable.
137  }
138  Scalar operator () (const int i, const int j, const int k, const int l) const
139  {
140  int uid = i + 3*j + 9*k + 27*l;
141  // [0, 0, 0, 0]
142  if (uid==0 ) return data[0];
143  // [0, 0, 0, 1]
144  else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
145  // [0, 0, 0, 2]
146  else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
147  // [0, 0, 1, 1]
148  else if (uid==36 || uid==30 || uid==12 || uid==28 || uid==10 || uid==4 ) return data[3];
149  // [0, 0, 1, 2]
150  else if (uid==63 || uid==45 || uid==57 || uid==21 || uid==33 || uid==15 || uid==55 || uid==19 || uid==7 || uid==29 || uid==11 || uid==5 ) return data[4];
151  // [0, 0, 2, 2]
152  else if (uid==72 || uid==60 || uid==24 || uid==56 || uid==20 || uid==8 ) return data[5];
153  // [0, 1, 1, 1]
154  else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[6];
155  // [0, 1, 1, 2]
156  else if (uid==66 || uid==48 || uid==42 || uid==64 || uid==46 || uid==58 || uid==22 || uid==34 || uid==16 || uid==38 || uid==32 || uid==14 ) return data[7];
157  // [0, 1, 2, 2]
158  else if (uid==75 || uid==69 || uid==51 || uid==73 || uid==61 || uid==25 || uid==65 || uid==47 || uid==59 || uid==23 || uid==35 || uid==17 ) return data[8];
159  // [0, 2, 2, 2]
160  else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[9];
161  // [1, 1, 1, 1]
162  else if (uid==40 ) return data[10];
163  // [1, 1, 1, 2]
164  else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[11];
165  // [1, 1, 2, 2]
166  else if (uid==76 || uid==70 || uid==52 || uid==68 || uid==50 || uid==44 ) return data[12];
167  // [1, 2, 2, 2]
168  else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[13];
169  // [2, 2, 2, 2]
170  else if (uid==80 ) return data[14];
171  else Util::Abort(INFO,"uid not in range (uid=",uid,")");
172  return data[0]; // this return statement is unreachable.
173  }
174  void Print (std::ostream& os)
175  {
176  for (int i = 0; i < 14; i++)
177  os << "i = " << i << " " << data[i] << std::endl;
178 
179  os.precision(4);
180  for (int k = 0; k < 3; k++)
181  {
182  for (int i = -1; i < 4; i++)
183  {
184  for (int l = 0; l < 3; l++)
185  {
186  if (i==-1) os << "┌ ┐ ";
187  else if (i == 3) os << "└ ┘ ";
188  else
189  {
190  os << "│ ";
191  for (int j = 0; j < 3; j++)
192  {
193  const Set::Scalar &val = (*this)(i,j,k,l);
194  os << std::scientific << std::setw(11) << val ; //(fabs(val)>1E-10 ? val : 0);
195  os << " ";
196  }
197  os << "│ ";
198  }
199  }
200  os<<std::endl;
201  }
202  }
203  }
204  static Matrix4<3,Sym::Full> Increment()
205  {
206  Matrix4<3,Sym::Full> ret;
207  for (int i = 0 ; i < 15; i++) ret.data[i] = (Set::Scalar)i;
208  return ret;
209  }
210  static Matrix4<3,Sym::Full> Randomize()
211  {
212  Matrix4<3,Sym::Full> ret;
213  for (int i = 0 ; i < 15; i++) ret.data[i] = Util::Random();
214  return ret;
215  }
216  static Matrix4<3,Sym::Full> Zero()
217  {
218  Matrix4<3,Sym::Full> ret;
219  for (int i = 0 ; i < 15; i++) ret.data[i] = 0.0;
220  return ret;
221  }
222 
223  bool contains_nan() const
224  {
225  if (std::isnan(data[ 0])) return true;
226  if (std::isnan(data[ 1])) return true;
227  if (std::isnan(data[ 2])) return true;
228  if (std::isnan(data[ 3])) return true;
229  if (std::isnan(data[ 4])) return true;
230  if (std::isnan(data[ 5])) return true;
231  if (std::isnan(data[ 6])) return true;
232  if (std::isnan(data[ 7])) return true;
233  if (std::isnan(data[ 8])) return true;
234  if (std::isnan(data[ 9])) return true;
235  if (std::isnan(data[10])) return true;
236  if (std::isnan(data[11])) return true;
237  if (std::isnan(data[12])) return true;
238  if (std::isnan(data[13])) return true;
239  if (std::isnan(data[14])) return true;
240  return false;
241  }
242 
243 };
244 std::ostream&
245 operator<< (std::ostream& os, const Matrix4<3,Sym::Full>& b);
246 }
247 #endif
Set::Full
@ Full
Definition: Base.H:197
Set::Matrix4< 2, Sym::Full >::contains_nan
bool contains_nan() const
Definition: Matrix4_Full.H:80
Set::Sym
Sym
Definition: Base.H:197
Util.H
Set::Matrix4< 3, Sym::Full >::contains_nan
bool contains_nan() const
Definition: Matrix4_Full.H:223
Set::Matrix4< 3, Sym::Full >::Print
void Print(std::ostream &os)
Definition: Matrix4_Full.H:174
Set::Matrix4< 2, Sym::Full >
Data structure for a symmetrix 4th order 3D tensor.
Definition: Matrix4_Full.H:21
Set::Matrix4< 3, Sym::Full >::Randomize
static Matrix4< 3, Sym::Full > Randomize()
Definition: Matrix4_Full.H:210
Set::Matrix4< 2, Sym::Full >::Print
void Print(std::ostream &)
Definition: Matrix4_Full.H:66
Base.H
Set::Matrix4< 2, Sym::Full >::Zero
static Matrix4< 2, Sym::Full > Zero()
Definition: Matrix4_Full.H:70
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< 3, Sym::Full >::Increment
static Matrix4< 3, Sym::Full > Increment()
Definition: Matrix4_Full.H:204
Set::operator<<
std::ostream & operator<<(std::ostream &os, const Matrix4< dim, sym > &b)
Definition: Base.H:201
Set::Matrix4< 3, Sym::Full >
Definition: Matrix4_Full.H:92
Set::Matrix4< 3, Sym::Full >::Zero
static Matrix4< 3, Sym::Full > Zero()
Definition: Matrix4_Full.H:216
Set::Matrix4< 3, Sym::Full >::Matrix4
AMREX_GPU_HOST_DEVICE Matrix4()
Definition: Matrix4_Full.H:99
Set::Matrix4< 2, Sym::Full >::Randomize
static Matrix4< 2, Sym::Full > Randomize()
Definition: Matrix4_Full.H:60
Util::Abort
void Abort(const char *msg)
Definition: Util.cpp:165
Set::Matrix4
Definition: Base.H:198
Set::Matrix4< 2, Sym::Full >::data
Scalar data[5]
Definition: Matrix4_Full.H:23
Set::Matrix4< 3, Sym::Full >::data
Scalar data[15]
Definition: Matrix4_Full.H:94
Set::Matrix4< 2, Sym::Full >::Matrix4
AMREX_GPU_HOST_DEVICE Matrix4()
Definition: Matrix4_Full.H:25
INFO
#define INFO
Definition: Util.H:20
Generated on Fri Aug 16 2024 16:30:29 for Alamo by   doxygen 1.8.17