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