Line data Source code
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 1362714 : 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 6815032 : int uid = i + 2*j + 4*k + 8*l;
31 : // [0, 0, 0, 0]
32 6814986 : if (uid==0 ) return data[0];
33 : // [0, 0, 0, 1]
34 5452228 : else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1];
35 : // [0, 0, 1, 1]
36 4089148 : else if (uid==12 || uid==10 || uid==6 || uid==9 || uid==5 || uid==3 ) return data[2];
37 : // [0, 1, 1, 1]
38 2725884 : else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[3];
39 : // [1, 1, 1, 1]
40 1362804 : else if (uid==15 ) return data[4];
41 0 : else Util::Abort(INFO,"Index out of range");
42 0 : return data[0]; // this code is unreachable
43 : }
44 6813560 : Scalar operator () (const int i, const int j, const int k, const int l) const
45 : {
46 6813560 : int uid = i + 2*j + 4*k + 8*l;
47 : // [0, 0, 0, 0]
48 6813560 : if (uid==0 ) return data[0];
49 : // [0, 0, 0, 1]
50 5450848 : else if (uid==8 || uid==4 || uid==2 || uid==1 ) return data[1];
51 : // [0, 0, 1, 1]
52 4088136 : else if (uid==12 || uid==10 || uid==6 || uid==9 || uid==5 || uid==3 ) return data[2];
53 : // [0, 1, 1, 1]
54 2725424 : else if (uid==14 || uid==13 || uid==11 || uid==7 ) return data[3];
55 : // [1, 1, 1, 1]
56 1362712 : else if (uid==15 ) return data[4];
57 0 : else Util::Abort(INFO,"Index out of range");
58 0 : return data[0]; // this code is unreachable
59 : }
60 2 : static Matrix4<2,Sym::Full> Randomize()
61 : {
62 2 : Matrix4<2,Sym::Full> ret;
63 12 : for (int i = 0 ; i < 5; i++) ret.data[i] = Util::Random();
64 2 : 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 2 : 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 7452 : int uid = i + 3*j + 9*k + 27*l;
105 : // [0, 0, 0, 0]
106 7452 : if (uid==0 ) return data[0];
107 : // [0, 0, 0, 1]
108 7360 : else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
109 : // [0, 0, 0, 2]
110 6992 : else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
111 : // [0, 0, 1, 1]
112 6624 : else if (uid==36 || uid==30 || uid==12 || uid==28 || uid==10 || uid==4 ) return data[3];
113 : // [0, 0, 1, 2]
114 6072 : 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 4968 : else if (uid==72 || uid==60 || uid==24 || uid==56 || uid==20 || uid==8 ) return data[5];
117 : // [0, 1, 1, 1]
118 4416 : else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[6];
119 : // [0, 1, 1, 2]
120 4048 : 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 2944 : 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 1840 : else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[9];
125 : // [1, 1, 1, 1]
126 1472 : else if (uid==40 ) return data[10];
127 : // [1, 1, 1, 2]
128 1380 : else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[11];
129 : // [1, 1, 2, 2]
130 1012 : else if (uid==76 || uid==70 || uid==52 || uid==68 || uid==50 || uid==44 ) return data[12];
131 : // [1, 2, 2, 2]
132 460 : else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[13];
133 : // [2, 2, 2, 2]
134 92 : else if (uid==80 ) return data[14];
135 0 : else Util::Abort(INFO,"uid not in range (uid=",uid,")");
136 0 : return data[0]; // this return statement is unreachable.
137 : }
138 0 : Scalar operator () (const int i, const int j, const int k, const int l) const
139 : {
140 0 : int uid = i + 3*j + 9*k + 27*l;
141 : // [0, 0, 0, 0]
142 0 : if (uid==0 ) return data[0];
143 : // [0, 0, 0, 1]
144 0 : else if (uid==27 || uid==9 || uid==3 || uid==1 ) return data[1];
145 : // [0, 0, 0, 2]
146 0 : else if (uid==54 || uid==18 || uid==6 || uid==2 ) return data[2];
147 : // [0, 0, 1, 1]
148 0 : else if (uid==36 || uid==30 || uid==12 || uid==28 || uid==10 || uid==4 ) return data[3];
149 : // [0, 0, 1, 2]
150 0 : 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 0 : else if (uid==72 || uid==60 || uid==24 || uid==56 || uid==20 || uid==8 ) return data[5];
153 : // [0, 1, 1, 1]
154 0 : else if (uid==39 || uid==37 || uid==31 || uid==13 ) return data[6];
155 : // [0, 1, 1, 2]
156 0 : 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 0 : 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 0 : else if (uid==78 || uid==74 || uid==62 || uid==26 ) return data[9];
161 : // [1, 1, 1, 1]
162 0 : else if (uid==40 ) return data[10];
163 : // [1, 1, 1, 2]
164 0 : else if (uid==67 || uid==49 || uid==43 || uid==41 ) return data[11];
165 : // [1, 1, 2, 2]
166 0 : else if (uid==76 || uid==70 || uid==52 || uid==68 || uid==50 || uid==44 ) return data[12];
167 : // [1, 2, 2, 2]
168 0 : else if (uid==79 || uid==77 || uid==71 || uid==53 ) return data[13];
169 : // [2, 2, 2, 2]
170 0 : else if (uid==80 ) return data[14];
171 0 : else Util::Abort(INFO,"uid not in range (uid=",uid,")");
172 0 : 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 2 : static Matrix4<3,Sym::Full> Randomize()
211 : {
212 2 : Matrix4<3,Sym::Full> ret;
213 32 : for (int i = 0 ; i < 15; i++) ret.data[i] = Util::Random();
214 2 : 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
|
