FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
struct_index_mapping.hpp
1// FEAT3: Finite Element Analysis Toolbox, Version 3
2// Copyright (C) 2010 by Stefan Turek & the FEAT group
3// FEAT3 is released under the GNU General Public License version 3,
4// see the file 'copyright.txt' in the top level directory for details.
5
6#pragma once
7
8// includes, FEAT
9#include <kernel/base_header.hpp>
10
11namespace FEAT
12{
13 namespace Geometry
14 {
16 namespace Intern
17 {
18 template<
19 int shape_dim_,
20 int cell_dim_,
21 int face_dim_>
22 struct StructIndexMapping;
23
24
25 // one-dimensional
26
28 template<>
29 struct StructIndexMapping<1, 1, 0>
30 {
31 static Index compute(Index i, Index j, const Index /*num_slices*/[])
32 {
33 return i + j;
34 }
35 }; //<1, 1, 0>
36
37
38 // quad mesh
39
41 template<>
42 struct StructIndexMapping<2, 1, 0>
43 {
44 static Index compute(Index i, Index j, const Index num_slices[])
45 {
46 // compute index
47 if(i < num_slices[0] * (num_slices[1] + 1))
48 {
49 return i + j + i / num_slices[0];
50 }
51 else
52 {
53 return ((i - num_slices[0] * (num_slices[1] + 1)) / num_slices[1])
54 + ((i - num_slices[0] * (num_slices[1] + 1)) % num_slices[1])
55 * (num_slices[0] + 1) + j * (num_slices[0] + 1);
56 }
57 }
58 }; //<2, 1, 0>
59
61 template<>
62 struct StructIndexMapping<2, 2, 0>
63 {
64 static Index compute(Index i, Index j, const Index num_slices[])
65 {
66 // compute index
67 return ((i / num_slices[0]) + ((j >> 1) & 0x1)) * (num_slices[0] + 1) + (i % num_slices[0]) + (j & 0x1);
68 }
69 }; //<2, 2, 0>
70
72 template<>
73 struct StructIndexMapping<2, 2, 1>
74 {
75 static Index compute(Index i, Index j, const Index num_slices[])
76 {
77 // compute index
78 return ((3 - j) / 2) * (i + j * num_slices[0]) + (j / 2) * (num_slices[0] * (num_slices[1] + 1)
79 + ((i % num_slices[0]) + (j / 3)) * num_slices[1] + (i / num_slices[0]));
80 }
81 }; //<2, 2, 1>
82
83
84 // hexahedron mesh
85
87 template<>
88 struct StructIndexMapping<3, 1, 0>
89 {
90 static Index compute(Index i, Index j, const Index num_slices[])
91 {
92 // auxiliary index
93 Index pos;
94
95 // compute auxiliary variable
96
97 //z-direction
98 if( i >= num_slices[1] * (num_slices[0] + 1) * (num_slices[2] + 1) +
99 (num_slices[2] + 1) * (num_slices[1] + 1) * num_slices[0])
100 {
101 pos = i - num_slices[1] * (num_slices[0] + 1) * (num_slices[2] + 1) - (num_slices[2] + 1)
102 * (num_slices[1] + 1) * num_slices[0];
103 return (pos / num_slices[2]) + (pos % num_slices[2] + j) * (num_slices[0] + 1)*(num_slices[1] + 1);
104 }
105 //y-direction
106 else if( i >= (num_slices[2] + 1) * (num_slices[1] + 1) * num_slices[0])
107 {
108 pos = i - (num_slices[2] + 1) * (num_slices[1] + 1) * num_slices[0];
109 Index x = (pos / num_slices[1]) % (num_slices[0] + 1);
110 Index y = pos % num_slices[1];
111 Index z = (pos / num_slices[1]) / (num_slices[0] + 1);
112 return z * (num_slices[0] + 1) * (num_slices[1] + 1) + (y + j) * (num_slices[0] + 1) + x;
113 }
114 // x-direction
115 else
116 {
117 return i + j + i / num_slices[0];
118 }
119 }
120 }; //<3, 1, 0>
121
123 template<>
124 struct StructIndexMapping<3, 2, 0>
125 {
126 static Index compute(Index i, Index j, const Index num_slices[])
127 {
128 // compute index
129 Index pos;
130 if( i < num_slices[0] * num_slices[1] * (num_slices[2] + 1))
131 {
132 return i + (i / num_slices[0]) + (i / (num_slices[0] * num_slices[1])) * (num_slices[0] + 1)
133 + (j % 2) + (j / 2) * (num_slices[0] + 1);
134 }
135 else if( i >= num_slices[0] * num_slices[1] * (num_slices[2] + 1)
136 && i < num_slices[0] * (num_slices[1] * (num_slices[2] + 1) + num_slices[2] * (num_slices[1] + 1)))
137 {
138 pos = i - num_slices[0] * num_slices[1] * (num_slices[2] + 1);
139 return pos % num_slices[0] + (pos / (num_slices[0] * num_slices[2])) * (num_slices[0] + 1)
140 + (pos % (num_slices[0] * num_slices[2])) / num_slices[0]
141 *(num_slices[0] + 1)*(num_slices[1] + 1) + (j % 2) + (j / 2) * (num_slices[0] + 1) * (num_slices[1] + 1);
142 }
143 else
144 {
145 pos = i - num_slices[0] * (num_slices[1] * (num_slices[2] + 1) + num_slices[2] * (num_slices[1] + 1));
146 return pos / (num_slices[1] * num_slices[2]) + (pos % num_slices[1]) * (num_slices[0] + 1) + ((pos / num_slices[1]) % num_slices[2])
147 * (num_slices[0] + 1) * (num_slices[1] + 1) + (j % 2) * (num_slices[0] + 1) + (j / 2)
148 * (num_slices[0] + 1) * (num_slices[1] + 1);
149 }
150 }
151 }; //<3, 2, 0>
152
154 template<>
155 struct StructIndexMapping<3, 3, 0>
156 {
157 static Index compute(Index i, Index j, const Index num_slices[])
158 {
159 // compute index
160 return (i / (num_slices[0] * num_slices[1]) + j / 4) * ((num_slices[0] + 1)*(num_slices[1] + 1))
161 + ((i % (num_slices[0] * num_slices[1])) / num_slices[0] + (j % 4) / 2) * (num_slices[0] + 1)
162 + (i % (num_slices[0] * num_slices[1])) % num_slices[0] + j % 2;
163 }
164 }; //<3, 3, 0>
165
167 template<>
168 struct StructIndexMapping<3, 2, 1>
169 {
170 static Index compute(Index i, Index j, const Index num_slices[])
171 {
172 // auxiliary indices
173 Index pos, node_0, node_0_x, node_0_y, node_0_z;
174
175 // xy plain
176 if( i < num_slices[0] * num_slices[1] * (num_slices[2] + 1))
177 {
178 // node 0
179 node_0 = i + (i / num_slices[0]) + (i / (num_slices[0] * num_slices[1])) * (num_slices[0] + 1);
180
181 if (j/2 == 0) //x direction edge
182 {
183 return node_0 - (node_0 / (num_slices[0] + 1)) + (j % 2) * num_slices[0];
184 }
185 else // y direction edge
186 {
187 // mesh coordinates
188 node_0_x = node_0 % (num_slices[0] + 1);
189 node_0_y = (node_0 / (num_slices[0] + 1)) % (num_slices[1] + 1);
190 node_0_z = node_0 / ((num_slices[0] + 1) * (num_slices[1] + 1));
191 return num_slices[0] * (num_slices[1] + 1) * (num_slices[2] + 1) + node_0_y + (node_0_x + j - 2)
192 * num_slices[1] + node_0_z * (num_slices[0] + 1) * num_slices[1];
193 }
194 }
195 // xz plain
196 else if( i >= num_slices[0] * num_slices[1] * (num_slices[2]+1)
197 && i < num_slices[0] * (num_slices[1] * (num_slices[2] + 1) + num_slices[2] * (num_slices[1] + 1)))
198 {
199 // position and node 0
200 pos = i - num_slices[0] * num_slices[1] * (num_slices[2] + 1);
201 node_0 = pos % num_slices[0] + (pos / (num_slices[0] * num_slices[2])) * (num_slices[0] + 1)
202 + (pos % (num_slices[0] * num_slices[2])) / num_slices[0]
203 * (num_slices[0] + 1) * (num_slices[1] + 1);
204
205 if (j/2 == 0) //x direction edge
206 {
207 return node_0 - (node_0 / (num_slices[0] + 1)) + (j % 2) * num_slices[0] * (num_slices[1] + 1);
208 }
209 else // z direction edge
210 {
211 // coordinates
212 node_0_x = node_0 % (num_slices[0] + 1);
213 node_0_y = (node_0 / (num_slices[0] + 1)) % (num_slices[1] + 1);
214 node_0_z = node_0 / ((num_slices[0] + 1) * (num_slices[1] + 1));
215 return num_slices[0] * (num_slices[1] + 1) * (num_slices[2] + 1)
216 + num_slices[1] * (num_slices[0] + 1) * (num_slices[2] + 1)
217 + node_0_z + (node_0_x + j - 2) * num_slices[2] + node_0_y * (num_slices[0] + 1) * num_slices[2];
218 }
219 }
220 // yz plain
221 else
222 {
223 // position and node 0
224 pos = i - num_slices[0] * (num_slices[1] * (num_slices[2] + 1) + num_slices[2] * (num_slices[1] + 1));
225 node_0 = pos / (num_slices[1] * num_slices[2]) + (pos % num_slices[1]) * (num_slices[0] + 1)
226 + ((pos / num_slices[1]) % num_slices[2]) * (num_slices[0] + 1) * (num_slices[1] + 1);
227
228 // coordinates
229 node_0_x = node_0 % (num_slices[0] + 1);
230 node_0_y = (node_0 / (num_slices[0] + 1)) % (num_slices[1] + 1);
231 node_0_z = node_0 / ((num_slices[0] + 1) * (num_slices[1] + 1));
232
233 if (j/2 == 0) //y direction edge
234 {
235 return num_slices[0] * (num_slices[1] + 1) * (num_slices[2] + 1) + node_0_y + node_0_x * num_slices[1]
236 + (node_0_z + j) * (num_slices[0] + 1) * num_slices[1];
237 }
238 else // z direction edge
239 {
240 return num_slices[0] * (num_slices[1] + 1) * (num_slices[2] + 1) + num_slices[1] * (num_slices[0] + 1) * (num_slices[2] + 1)
241 + node_0_z + node_0_x * num_slices[2] + (node_0_y + j - 2) * (num_slices[0] + 1) * num_slices[2];
242 }
243 }//if
244
245 }
246 }; //<3, 2, 1>
247
248
250 template<>
251 struct StructIndexMapping<3, 3, 1>
252 {
253 static Index compute(Index i, Index j, const Index num_slices[])
254 {
255 // compute index
256 Index x = i % (num_slices[0] * num_slices[1]) % num_slices[0];
257 Index y = i % (num_slices[0] * num_slices[1]) / num_slices[0];
258 Index z = i / (num_slices[0] * num_slices[1]);
259 if( j / 4 == 0)
260 {
261 return (z + j / 2) * num_slices[0] * (num_slices[1] + 1) + (y + j % 2) * num_slices[0] + x;
262 }
263 else if( j / 4 == 1)
264 {
265 return (num_slices[2] + 1) * (num_slices[1] + 1) * num_slices[0] + (z + j / 6) * num_slices[1] * (num_slices[0]+1)
266 + (x + j % 2) * num_slices[1] + y;
267 }
268 else
269 {
270 return (num_slices[2] + 1) * (num_slices[1] + 1) * num_slices[0] + num_slices[1] * (num_slices[0] + 1)
271 * (num_slices[2] + 1) + (y + j / 10) * num_slices[2] * (num_slices[0] + 1) + (x + j % 2) * num_slices[2] + z;
272 }
273 }
274 }; //<3, 3, 1>
275
277 template<>
278 struct StructIndexMapping<3, 3, 2>
279 {
280 static Index compute(Index i, Index j, const Index num_slices[])
281 {
282
283 if(j/2 == 0) // x surface
284 {
285 return i + j*num_slices[0] * num_slices[1];
286 }
287 else if( j/4 == 0) // y surface
288 {
289 // compute index
290 Index x = i % (num_slices[0] * num_slices[1]) % num_slices[0];
291 Index y = i % (num_slices[0] * num_slices[1]) / num_slices[0];
292 Index z = i / (num_slices[0] * num_slices[1]);
293 return num_slices[0] * num_slices[1] * (num_slices[2] + 1)
294 + x + z * num_slices[0] + (y + j - 2) * num_slices[0] * num_slices[2];
295 }
296 else // z surface
297 {
298 // compute index
299 Index x = i % (num_slices[0] * num_slices[1]) % num_slices[0];
300 Index y = i % (num_slices[0] * num_slices[1]) / num_slices[0];
301 Index z = i / (num_slices[0] * num_slices[1]);
302 return num_slices[0] * (num_slices[1] * (num_slices[2]+1) + num_slices[2] * (num_slices[1]+1))
303 + y + z * num_slices[1] + (x + j - 4) * num_slices[1] * num_slices[2];
304 }
305 }
306 }; //<3, 3, 2>
307
308 } // namespace Intern
310 } // namespace Geometry
311} // namespace FEAT
base_header.hpp
FEAT Kernel base header.
FEAT
FEAT namespace.
Definition: adjactor.hpp:12
FEAT::Index
std::uint64_t Index
Index data type.
Definition: base_header.hpp:122