FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
shape_convert_index.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/geometry/index_set.hpp>
10#include <kernel/geometry/intern/shape_convert_traits.hpp>
11#include <kernel/geometry/intern/entity_counter.hpp>
12
13namespace FEAT
14{
15 namespace Geometry
16 {
18 namespace Intern
19 {
20 template<
21 typename Shape_,
22 int cell_dim_>
23 struct ShapeConvertIndex;
24
25 /* *************************************************************************************** */
26 /* *************************************************************************************** */
27 /* ***************************** Hypercube --> Simplex *********************************** */
28 /* *************************************************************************************** */
29 /* *************************************************************************************** */
30
36 template<>
37 struct ShapeConvertIndex<Shape::Simplex<1>, 1>
38 {
39 static constexpr int shape_dim = 1;
40 static constexpr int cell_dim = 1;
41 typedef Shape::Simplex<shape_dim> ShapeType;
42 typedef Shape::Hypercube<shape_dim> OtherShapeType;
43 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
44 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
45 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
46
47 static Index refine(
48 IndexSetType& index_set_out,
49 const Index offset,
50 const Index index_offsets[],
51 const IndexSetHolderType& index_set_holder_in)
52 {
53 // fetch vertex index offsets
54 const Index iov = index_offsets[0];
55
56 // typedef index set vector reference for output index set
57 typedef IndexSetType::IndexTupleType IndexTupleType;
58
59 // typedef index set types
60 typedef IndexSet<2> IndexSetTypeIn;
61
62 // typedef index vector type
63 typedef IndexSetTypeIn::IndexTupleType IndexTupleTypeIn;
64
65 // fetch the vertices-at-edge index set
66 const IndexSetTypeIn& index_set_in = index_set_holder_in.get_index_set<1,0>();
67
68 // fetch number of coarse mesh edges
69 const Index num_edges = index_set_in.get_num_entities();
70
71 // loop over all coarse mesh quads
72 FEAT_PRAGMA_OMP(parallel for)
73 for(Index i = 0; i < num_edges; ++i)
74 {
75 // fetch coarse mesh vertices-at-quad index vector
76 const IndexTupleTypeIn& v_e = index_set_in[i];
77
78 // fetch fine mesh vertices-at-edge index vectors
79 IndexTupleType& e_0 = index_set_out[offset + i];
80
81 e_0[0] = iov + v_e[0];
82 e_0[1] = iov + v_e[1];
83 }
84
85 // return inner edge count
86 return num_edges;
87 }
88 };
89
95 template<>
96 struct ShapeConvertIndex<Shape::Simplex<2>, 1>
97 {
98 static constexpr int shape_dim = 2;
99 static constexpr int cell_dim = 1;
100 typedef Shape::Simplex<shape_dim> ShapeType;
101 typedef Shape::Hypercube<shape_dim> OtherShapeType;
102 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
103 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
104 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
105
106 static Index refine(
107 IndexSetType& index_set_out,
108 const Index offset,
109 const Index index_offsets[],
110 const IndexSetHolderType& index_set_holder_in)
111 {
112 // fetch vertex index offsets
113 const Index iov = index_offsets[0];
114 const Index ioq = index_offsets[2];
115
116 // typedef index set vector reference for output index set
117 typedef IndexSetType::IndexTupleType IndexTupleType;
118
119 // typedef index set types
120 typedef IndexSet<4> IndexSetTypeIn;
121
122 // typedef index vector type
123 typedef IndexSetTypeIn::IndexTupleType IndexTupleTypeIn;
124
125 // fetch the vertices-at-quad index set
126 const IndexSetTypeIn& index_set_in = index_set_holder_in.get_index_set<2,0>();
127
128 // fetch number of coarse mesh quads
129 const Index num_quads = index_set_in.get_num_entities();
130
131 // Each quad generates four inner edges (e_i) upon conversion,
132 // with each one connecting one of the quad's corner vertices (v_i)
133 // with the vertex generated from the quad midpoint (X).
134 //
135 // v_2-------------------------v_3
136 // |\__ __/|
137 // | \ / |
138 // | e_2__ __e_3 |
139 // | \__ __/ |
140 // | \ / |
141 // | X |
142 // | __/ \__ |
143 // | __/ \__ |
144 // | e_0 e_1 |
145 // | __/ \__ |
146 // |/ \|
147 // v_0-------------------------v_1
148
149
150 // loop over all coarse mesh quads
151 FEAT_PRAGMA_OMP(parallel for)
152 for(Index i = 0; i < num_quads; ++i)
153 {
154 // fetch coarse mesh vertices-at-quad index vector
155 const IndexTupleTypeIn& v_q = index_set_in[i];
156
157 // fetch fine mesh vertices-at-edge index vectors
158 IndexTupleType& e_0 = index_set_out[offset + 4*i + 0];
159 IndexTupleType& e_1 = index_set_out[offset + 4*i + 1];
160 IndexTupleType& e_2 = index_set_out[offset + 4*i + 2];
161 IndexTupleType& e_3 = index_set_out[offset + 4*i + 3];
162
163 e_0[0] = iov + v_q[0]; // v_0
164 e_0[1] = ioq + i; // X
165
166 e_1[0] = iov + v_q[1]; // v_1
167 e_1[1] = ioq + i; // X
168
169 e_2[0] = iov + v_q[2]; // v_2
170 e_2[1] = ioq + i; // X
171
172 e_3[0] = iov + v_q[3]; // v_3
173 e_3[1] = ioq + i; // X
174 }
175
176 // return inner edge count
177 return 4*num_quads;
178 }
179 };
180
186 template<>
187 struct ShapeConvertIndex<Shape::Simplex<3>, 1>
188 {
189 static constexpr int shape_dim = 3;
190 static constexpr int cell_dim = 1;
191 typedef Shape::Simplex<shape_dim> ShapeType;
192 typedef Shape::Hypercube<shape_dim> OtherShapeType;
193 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
194 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
195 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
196
197 static Index refine(
198 IndexSetType& index_set_out,
199 const Index offset,
200 const Index index_offsets[],
201 const IndexSetHolderType& index_set_holder_in)
202 {
203 // fetch vertex index offsets
204 const Index iov = index_offsets[0];
205 const Index ioq = index_offsets[2];
206 const Index ioc = index_offsets[3];
207
208 // typedef index set vector reference for output index set
209 typedef IndexSetType::IndexTupleType IndexTupleType;
210
211 // typedef index set types
212 typedef IndexSet<8> IndexSetTypeIn_v_c;
213 typedef IndexSet<6> IndexSetTypeIn_q_c;
214
215 // typedef index vector type
216 typedef IndexSetTypeIn_v_c::IndexTupleType IndexTupleType_v_c;
217 typedef IndexSetTypeIn_q_c::IndexTupleType IndexTupleType_q_c;
218
219 // fetch the vertices-at-hexa and the quad-at-hexa index set
220 const IndexSetTypeIn_v_c& index_set_in_v_c = index_set_holder_in.get_index_set<3,0>();
221 const IndexSetTypeIn_q_c& index_set_in_q_c = index_set_holder_in.get_index_set<3,2>();
222
223 // fetch number of coarse mesh hexas
224 const Index num_hexas = index_set_in_v_c.get_num_entities();
225
226 // loop over all coarse mesh hexas
227 FEAT_PRAGMA_OMP(parallel for)
228 for(Index i = 0; i < num_hexas; ++i)
229 {
230 // fetch coarse mesh vertices-at-hexa index vector
231 const IndexTupleType_v_c& v_c = index_set_in_v_c[i];
232 const IndexTupleType_q_c& q_c = index_set_in_q_c[i];
233
234 // fetch new mesh vertices-at-edge index vectors
235 IndexTupleType& e_0 = index_set_out[offset + 14*i + 0];
236 IndexTupleType& e_1 = index_set_out[offset + 14*i + 1];
237 IndexTupleType& e_2 = index_set_out[offset + 14*i + 2];
238 IndexTupleType& e_3 = index_set_out[offset + 14*i + 3];
239 IndexTupleType& e_4 = index_set_out[offset + 14*i + 4];
240 IndexTupleType& e_5 = index_set_out[offset + 14*i + 5];
241 IndexTupleType& e_6 = index_set_out[offset + 14*i + 6];
242 IndexTupleType& e_7 = index_set_out[offset + 14*i + 7];
243 IndexTupleType& e_8 = index_set_out[offset + 14*i + 8];
244 IndexTupleType& e_9 = index_set_out[offset + 14*i + 9];
245 IndexTupleType& e_10 = index_set_out[offset + 14*i + 10];
246 IndexTupleType& e_11 = index_set_out[offset + 14*i + 11];
247 IndexTupleType& e_12 = index_set_out[offset + 14*i + 12];
248 IndexTupleType& e_13 = index_set_out[offset + 14*i + 13];
249
250 // corner-midpoint edges
251 e_0[0] = iov + v_c[0]; // v_0
252 e_0[1] = ioc + i; // X
253
254 e_1[0] = iov + v_c[1]; // v_1
255 e_1[1] = ioc + i; // X
256
257 e_2[0] = iov + v_c[2]; // v_2
258 e_2[1] = ioc + i; // X
259
260 e_3[0] = iov + v_c[3]; // v_3
261 e_3[1] = ioc + i; // X
262
263 e_4[0] = iov + v_c[4]; // v_4
264 e_4[1] = ioc + i; // X
265
266 e_5[0] = iov + v_c[5]; // v_5
267 e_5[1] = ioc + i; // X
268
269 e_6[0] = iov + v_c[6]; // v_6
270 e_6[1] = ioc + i; // X
271
272 e_7[0] = iov + v_c[7]; // v_7
273 e_7[1] = ioc + i; // X
274
275
276 // face-midpoint edges
277 e_8[0] = ioq + q_c[0]; // q_0
278 e_8[1] = ioc + i; // X
279
280 e_9[0] = ioq + q_c[1]; // q_1
281 e_9[1] = ioc + i; // X
282
283 e_10[0] = ioq + q_c[2]; // q_2
284 e_10[1] = ioc + i; // X
285
286 e_11[0] = ioq + q_c[3]; // q_3
287 e_11[1] = ioc + i; // X
288
289 e_12[0] = ioq + q_c[4]; // q_4
290 e_12[1] = ioc + i; // X
291
292 e_13[0] = ioq + q_c[5]; // q_5
293 e_13[1] = ioc + i; // X
294
295 }
296
297 // return inner edge count
298 return 14*num_hexas;
299 }
300 };
301
307 template<>
308 struct ShapeConvertIndex<Shape::Simplex<2>, 2>
309 {
310 static constexpr int shape_dim = 2;
311 static constexpr int cell_dim = 2;
312 typedef Shape::Simplex<shape_dim> ShapeType;
313 typedef Shape::Hypercube<shape_dim> OtherShapeType;
314 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
315 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
316 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
317
318 static Index refine(
319 IndexSetType& index_set_out,
320 const Index offset,
321 const Index index_offsets[],
322 const IndexSetHolderType& index_set_holder_in)
323 {
324 // fetch vertex index offsets
325 const Index iov = index_offsets[0];
326 const Index ioq = index_offsets[2];
327
328 // typedef index set vector reference for output index set
329 typedef IndexSetType::IndexTupleType IndexTupleType;
330
331 // typedef index set types
332 typedef IndexSet<4> IndexSetTypeIn;
333
334 // typedef index vector type
335 typedef IndexSetTypeIn::IndexTupleType IndexTupleTypeIn;
336
337 // fetch the vertices-at-quad index set
338 const IndexSetTypeIn& index_set_in = index_set_holder_in.get_index_set<2,0>();
339
340 // fetch number of coarse mesh quads
341 const Index num_quads = index_set_in.get_num_entities();
342
343 // Each quad generates four inner triangles (t_i) upon conversion,
344 // with each one connecting two of the quad's corner vertices (v_i)
345 // with the vertex generated from the quad midpoint (X).
346 //
347 // v_2-------------------------v_3
348 // |\__ __/|
349 // | \__ t_1 __/ |
350 // | \__ __/ |
351 // | \__ __/ |
352 // | \ / |
353 // | t_2 X t_3 |
354 // | __/ \__ |
355 // | __/ \__ |
356 // | __/ \__ |
357 // | __/ t_0 \__ |
358 // |/ \|
359 // v_0-------------------------v_1
360
361
362 // loop over all coarse mesh quads
363 FEAT_PRAGMA_OMP(parallel for)
364 for(Index i = 0; i < num_quads; ++i)
365 {
366 // fetch coarse mesh vertices-at-quad index vector
367 const IndexTupleTypeIn& v_q = index_set_in[i];
368
369 // fetch fine mesh vertices-at-edge index vectors
370 IndexTupleType& t_0 = index_set_out[offset + 4*i + 0];
371 IndexTupleType& t_1 = index_set_out[offset + 4*i + 1];
372 IndexTupleType& t_2 = index_set_out[offset + 4*i + 2];
373 IndexTupleType& t_3 = index_set_out[offset + 4*i + 3];
374
375 t_0[0] = ioq + i; // X
376 t_0[1] = iov + v_q[0]; // v_0
377 t_0[2] = iov + v_q[1]; // v_1
378
379 t_1[0] = ioq + i; // X
380 t_1[1] = iov + v_q[3]; // v_3
381 t_1[2] = iov + v_q[2]; // v_2
382
383 t_2[0] = ioq + i; // X
384 t_2[1] = iov + v_q[2]; // v_2
385 t_2[2] = iov + v_q[0]; // v_0
386
387 t_3[0] = ioq + i; // X
388 t_3[1] = iov + v_q[1]; // v_1
389 t_3[2] = iov + v_q[3]; // v_3
390 }
391
392 // return triangle count
393 return 4*num_quads;
394 }
395 };
396
397
403 template<>
404 struct ShapeConvertIndex<Shape::Simplex<3>, 2>
405 {
406 static constexpr int shape_dim = 3;
407 static constexpr int cell_dim = 2;
408 typedef Shape::Simplex<shape_dim> ShapeType;
409 typedef Shape::Hypercube<shape_dim> OtherShapeType;
410 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
411 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
412 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
413
414 static Index refine(
415 IndexSetType& index_set_out,
416 const Index offset,
417 const Index index_offsets[],
418 const IndexSetHolderType& index_set_holder_in)
419 {
420 // fetch vertex index offsets
421 const Index iov = index_offsets[0];
422 const Index ioq = index_offsets[2];
423 const Index ioc = index_offsets[3];
424
425 // typedef index set vector reference for output index set
426 typedef IndexSetType::IndexTupleType IndexTupleType;
427
428 // typedef index set types
429 typedef IndexSet<8> IndexSetTypeIn_v_c;
430 typedef IndexSet<6> IndexSetTypeIn_q_c;
431
432 // typedef index vector type
433 typedef IndexSetTypeIn_v_c::IndexTupleType IndexTupleType_v_c;
434 typedef IndexSetTypeIn_q_c::IndexTupleType IndexTupleType_q_c;
435
436 // fetch the vertices-at-hexa and quad-at-hexa index set
437 const IndexSetTypeIn_v_c& index_set_in_v_c = index_set_holder_in.get_index_set<3,0>();
438 const IndexSetTypeIn_q_c& index_set_in_q_c = index_set_holder_in.get_index_set<3,2>();
439
440 // fetch number of coarse mesh hexas
441 const Index num_cubes = index_set_in_v_c.get_num_entities();
442
443 // loop over all coarse mesh hexas
444 FEAT_PRAGMA_OMP(parallel for)
445 for(Index i = 0; i < num_cubes; ++i)
446 {
447 // fetch coarse mesh vertices-at-hexas and quad-at-hexa index vector
448 const IndexTupleType_v_c& v_c = index_set_in_v_c[i];
449 const IndexTupleType_q_c& q_c = index_set_in_q_c[i];
450
451 // fetch fine mesh vertices-at-triangle index vectors
452 IndexTupleType& t_0 = index_set_out[offset + 36*i + 0];
453 IndexTupleType& t_1 = index_set_out[offset + 36*i + 1];
454 IndexTupleType& t_2 = index_set_out[offset + 36*i + 2];
455 IndexTupleType& t_3 = index_set_out[offset + 36*i + 3];
456 IndexTupleType& t_4 = index_set_out[offset + 36*i + 4];
457 IndexTupleType& t_5 = index_set_out[offset + 36*i + 5];
458 IndexTupleType& t_6 = index_set_out[offset + 36*i + 6];
459 IndexTupleType& t_7 = index_set_out[offset + 36*i + 7];
460 IndexTupleType& t_8 = index_set_out[offset + 36*i + 8];
461 IndexTupleType& t_9 = index_set_out[offset + 36*i + 9];
462 IndexTupleType& t_10 = index_set_out[offset + 36*i + 10];
463 IndexTupleType& t_11 = index_set_out[offset + 36*i + 11];
464 IndexTupleType& t_12 = index_set_out[offset + 36*i + 12];
465 IndexTupleType& t_13 = index_set_out[offset + 36*i + 13];
466 IndexTupleType& t_14 = index_set_out[offset + 36*i + 14];
467 IndexTupleType& t_15 = index_set_out[offset + 36*i + 15];
468 IndexTupleType& t_16 = index_set_out[offset + 36*i + 16];
469 IndexTupleType& t_17 = index_set_out[offset + 36*i + 17];
470 IndexTupleType& t_18 = index_set_out[offset + 36*i + 18];
471 IndexTupleType& t_19 = index_set_out[offset + 36*i + 19];
472 IndexTupleType& t_20 = index_set_out[offset + 36*i + 20];
473 IndexTupleType& t_21 = index_set_out[offset + 36*i + 21];
474 IndexTupleType& t_22 = index_set_out[offset + 36*i + 22];
475 IndexTupleType& t_23 = index_set_out[offset + 36*i + 23];
476 IndexTupleType& t_24 = index_set_out[offset + 36*i + 24];
477 IndexTupleType& t_25 = index_set_out[offset + 36*i + 25];
478 IndexTupleType& t_26 = index_set_out[offset + 36*i + 26];
479 IndexTupleType& t_27 = index_set_out[offset + 36*i + 27];
480 IndexTupleType& t_28 = index_set_out[offset + 36*i + 28];
481 IndexTupleType& t_29 = index_set_out[offset + 36*i + 29];
482 IndexTupleType& t_30 = index_set_out[offset + 36*i + 30];
483 IndexTupleType& t_31 = index_set_out[offset + 36*i + 31];
484 IndexTupleType& t_32 = index_set_out[offset + 36*i + 32];
485 IndexTupleType& t_33 = index_set_out[offset + 36*i + 33];
486 IndexTupleType& t_34 = index_set_out[offset + 36*i + 34];
487 IndexTupleType& t_35 = index_set_out[offset + 36*i + 35];
488
489 // coarse edge triangles
490 t_0[0] = ioc + i; // X
491 t_0[1] = iov + v_c[0]; // v_0
492 t_0[2] = iov + v_c[1]; // v_1
493
494 t_1[0] = ioc + i; // X
495 t_1[1] = iov + v_c[3]; // v_3
496 t_1[2] = iov + v_c[2]; // v_2
497
498 t_2[0] = ioc + i; // X
499 t_2[1] = iov + v_c[4]; // v_4
500 t_2[2] = iov + v_c[5]; // v_5
501
502 t_3[0] = ioc + i; // X
503 t_3[1] = iov + v_c[7]; // v_7
504 t_3[2] = iov + v_c[6]; // v_6
505
506 t_4[0] = ioc + i; // X
507 t_4[1] = iov + v_c[2]; // v_2
508 t_4[2] = iov + v_c[0]; // v_0
509
510 t_5[0] = ioc + i; // X
511 t_5[1] = iov + v_c[1]; // v_1
512 t_5[2] = iov + v_c[3]; // v_3
513
514 t_6[0] = ioc + i; // X
515 t_6[1] = iov + v_c[6]; // v_6
516 t_6[2] = iov + v_c[4]; // v_4
517
518 t_7[0] = ioc + i; // X
519 t_7[1] = iov + v_c[5]; // v_5
520 t_7[2] = iov + v_c[7]; // v_7
521
522 t_8[0] = ioc + i; // X
523 t_8[1] = iov + v_c[4]; // v_4
524 t_8[2] = iov + v_c[0]; // v_0
525
526 t_9[0] = ioc + i; // X
527 t_9[1] = iov + v_c[1]; // v_1
528 t_9[2] = iov + v_c[5]; // v_5
529
530 t_10[0] = ioc + i; // X
531 t_10[1] = iov + v_c[6]; // v_6
532 t_10[2] = iov + v_c[2]; // v_2
533
534 t_11[0] = ioc + i; // X
535 t_11[1] = iov + v_c[3]; // v_3
536 t_11[2] = iov + v_c[7]; // v_7
537
538 // face triangles
539 // face 0
540 t_12[0] = ioc + i; // X
541 t_12[1] = ioq + q_c[0]; // q_0
542 t_12[2] = iov + v_c[0]; // v_0
543
544 t_13[0] = ioc + i; // X
545 t_13[1] = ioq + q_c[0]; // q_0
546 t_13[2] = iov + v_c[1]; // v_1
547
548 t_14[0] = ioc + i; // X
549 t_14[1] = ioq + q_c[0]; // q_0
550 t_14[2] = iov + v_c[2]; // v_2
551
552 t_15[0] = ioc + i; // X
553 t_15[1] = ioq + q_c[0]; // q_0
554 t_15[2] = iov + v_c[3]; // v_3
555
556 // face 1
557 t_16[0] = ioc + i; // X
558 t_16[1] = ioq + q_c[1]; // q_1
559 t_16[2] = iov + v_c[4]; // v_4
560
561 t_17[0] = ioc + i; // X
562 t_17[1] = ioq + q_c[1]; // q_1
563 t_17[2] = iov + v_c[5]; // v_5
564
565 t_18[0] = ioc + i; // X
566 t_18[1] = ioq + q_c[1]; // q_1
567 t_18[2] = iov + v_c[6]; // v_6
568
569 t_19[0] = ioc + i; // X
570 t_19[1] = ioq + q_c[1]; // q_1
571 t_19[2] = iov + v_c[7]; // v_7
572
573 // face 2
574 t_20[0] = ioc + i; // X
575 t_20[1] = ioq + q_c[2]; // q_2
576 t_20[2] = iov + v_c[0]; // v_0
577
578 t_21[0] = ioc + i; // X
579 t_21[1] = ioq + q_c[2]; // q_2
580 t_21[2] = iov + v_c[1]; // v_1
581
582 t_22[0] = ioc + i; // X
583 t_22[1] = ioq + q_c[2]; // q_2
584 t_22[2] = iov + v_c[4]; // v_4
585
586 t_23[0] = ioc + i; // X
587 t_23[1] = ioq + q_c[2]; // q_2
588 t_23[2] = iov + v_c[5]; // v_5
589
590 // face 3
591 t_24[0] = ioc + i; // X
592 t_24[1] = ioq + q_c[3]; // q_3
593 t_24[2] = iov + v_c[2]; // v_2
594
595 t_25[0] = ioc + i; // X
596 t_25[1] = ioq + q_c[3]; // q_3
597 t_25[2] = iov + v_c[3]; // v_3
598
599 t_26[0] = ioc + i; // X
600 t_26[1] = ioq + q_c[3]; // q_3
601 t_26[2] = iov + v_c[6]; // v_6
602
603 t_27[0] = ioc + i; // X
604 t_27[1] = ioq + q_c[3]; // q_3
605 t_27[2] = iov + v_c[7]; // v_7
606
607 // face 4
608 t_28[0] = ioc + i; // X
609 t_28[1] = ioq + q_c[4]; // q_4
610 t_28[2] = iov + v_c[0]; // v_0
611
612 t_29[0] = ioc + i; // X
613 t_29[1] = ioq + q_c[4]; // q_4
614 t_29[2] = iov + v_c[2]; // v_2
615
616 t_30[0] = ioc + i; // X
617 t_30[1] = ioq + q_c[4]; // q_4
618 t_30[2] = iov + v_c[4]; // v_4
619
620 t_31[0] = ioc + i; // X
621 t_31[1] = ioq + q_c[4]; // q_4
622 t_31[2] = iov + v_c[6]; // v_6
623
624 // face 5
625 t_32[0] = ioc + i; // X
626 t_32[1] = ioq + q_c[5]; // q_5
627 t_32[2] = iov + v_c[1]; // v_1
628
629 t_33[0] = ioc + i; // X
630 t_33[1] = ioq + q_c[5]; // q_5
631 t_33[2] = iov + v_c[3]; // v_3
632
633 t_34[0] = ioc + i; // X
634 t_34[1] = ioq + q_c[5]; // q_5
635 t_34[2] = iov + v_c[5]; // v_5
636
637 t_35[0] = ioc + i; // X
638 t_35[1] = ioq + q_c[5]; // q_5
639 t_35[2] = iov + v_c[7]; // v_7
640
641 }
642
643 // return triangle count
644 return 36*num_cubes;
645 }
646 };
647
648
655 template<>
656 struct ShapeConvertIndex<Shape::Simplex<3>, 3>
657 {
658 static constexpr int shape_dim = 3;
659 static constexpr int cell_dim = 3;
660 typedef Shape::Simplex<shape_dim> ShapeType;
661 typedef Shape::Hypercube<shape_dim> OtherShapeType;
662 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
663 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
664 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
665
666 static Index refine(
667 IndexSetType& index_set_out,
668 const Index offset,
669 const Index index_offsets[],
670 const IndexSetHolderType& index_set_holder_in)
671 {
672 // fetch vertex index offsets
673 const Index iov = index_offsets[0];
674 const Index ioq = index_offsets[2];
675 const Index ioc = index_offsets[3];
676
677 // typedef index set vector reference for output index set
678 typedef IndexSetType::IndexTupleType IndexTupleType;
679
680 // typedef index set types
681 typedef IndexSet<8> IndexSetTypeIn_v_c;
682 typedef IndexSet<6> IndexSetTypeIn_q_c;
683
684 // typedef index vector type
685 typedef IndexSetTypeIn_v_c::IndexTupleType IndexTupleType_v_c;
686 typedef IndexSetTypeIn_q_c::IndexTupleType IndexTupleType_q_c;
687
688 // fetch the vertices-at-hexas and quad-at-hexas index set
689 const IndexSetTypeIn_v_c& index_set_in_v_c = index_set_holder_in.get_index_set<3,0>();
690 const IndexSetTypeIn_q_c& index_set_in_q_c = index_set_holder_in.get_index_set<3,2>();
691
692 // fetch number of coarse mesh hexas
693 const Index num_cubes = index_set_in_v_c.get_num_entities();
694
695 // loop over all coarse mesh hexas
696 FEAT_PRAGMA_OMP(parallel for)
697 for(Index i = 0; i < num_cubes; ++i)
698 {
699 // fetch coarse mesh vertices-at-hexas and quad-at-hexas index vector
700 const IndexTupleType_v_c& v_c = index_set_in_v_c[i];
701 const IndexTupleType_q_c& q_c = index_set_in_q_c[i];
702
703 // fetch fine mesh vertices-at-tetra index vectors
704 IndexTupleType& t_0 = index_set_out[offset + 24*i + 0];
705 IndexTupleType& t_1 = index_set_out[offset + 24*i + 1];
706 IndexTupleType& t_2 = index_set_out[offset + 24*i + 2];
707 IndexTupleType& t_3 = index_set_out[offset + 24*i + 3];
708 IndexTupleType& t_4 = index_set_out[offset + 24*i + 4];
709 IndexTupleType& t_5 = index_set_out[offset + 24*i + 5];
710 IndexTupleType& t_6 = index_set_out[offset + 24*i + 6];
711 IndexTupleType& t_7 = index_set_out[offset + 24*i + 7];
712 IndexTupleType& t_8 = index_set_out[offset + 24*i + 8];
713 IndexTupleType& t_9 = index_set_out[offset + 24*i + 9];
714 IndexTupleType& t_10 = index_set_out[offset + 24*i + 10];
715 IndexTupleType& t_11 = index_set_out[offset + 24*i + 11];
716 IndexTupleType& t_12 = index_set_out[offset + 24*i + 12];
717 IndexTupleType& t_13 = index_set_out[offset + 24*i + 13];
718 IndexTupleType& t_14 = index_set_out[offset + 24*i + 14];
719 IndexTupleType& t_15 = index_set_out[offset + 24*i + 15];
720 IndexTupleType& t_16 = index_set_out[offset + 24*i + 16];
721 IndexTupleType& t_17 = index_set_out[offset + 24*i + 17];
722 IndexTupleType& t_18 = index_set_out[offset + 24*i + 18];
723 IndexTupleType& t_19 = index_set_out[offset + 24*i + 19];
724 IndexTupleType& t_20 = index_set_out[offset + 24*i + 20];
725 IndexTupleType& t_21 = index_set_out[offset + 24*i + 21];
726 IndexTupleType& t_22 = index_set_out[offset + 24*i + 22];
727 IndexTupleType& t_23 = index_set_out[offset + 24*i + 23];
728
729 // end up with X and use convention for tetrahedron index numbering (see wiki)
730 // face 0
731 t_0[0] = iov + v_c[0]; // v_0
732 t_0[1] = iov + v_c[1]; // v_1
733 t_0[2] = ioq + q_c[0]; // q_0
734 t_0[3] = ioc + i; // X
735
736 t_1[0] = iov + v_c[3]; // v_3
737 t_1[1] = iov + v_c[2]; // v_2
738 t_1[2] = ioq + q_c[0]; // q_0
739 t_1[3] = ioc + i; // X
740
741 t_2[0] = iov + v_c[2]; // v_2
742 t_2[1] = iov + v_c[0]; // v_0
743 t_2[2] = ioq + q_c[0]; // q_0
744 t_2[3] = ioc + i; // X
745
746 t_3[0] = iov + v_c[1]; // v_1
747 t_3[1] = iov + v_c[3]; // v_3
748 t_3[2] = ioq + q_c[0]; // q_0
749 t_3[3] = ioc + i; // X
750
751 // face 1
752 t_4[0] = iov + v_c[5]; // v_5
753 t_4[1] = iov + v_c[4]; // v_4
754 t_4[2] = ioq + q_c[1]; // q_1
755 t_4[3] = ioc + i; // X
756
757 t_5[0] = iov + v_c[6]; // v_6
758 t_5[1] = iov + v_c[7]; // v_7
759 t_5[2] = ioq + q_c[1]; // q_1
760 t_5[3] = ioc + i; // X
761
762 t_6[0] = iov + v_c[4]; // v_4
763 t_6[1] = iov + v_c[6]; // v_6
764 t_6[2] = ioq + q_c[1]; // q_1
765 t_6[3] = ioc + i; // X
766
767 t_7[0] = iov + v_c[7]; // v_7
768 t_7[1] = iov + v_c[5]; // v_5
769 t_7[2] = ioq + q_c[1]; // q_1
770 t_7[3] = ioc + i; // X
771
772 // face 2
773 t_8[0] = iov + v_c[1]; // v_1
774 t_8[1] = iov + v_c[0]; // v_0
775 t_8[2] = ioq + q_c[2]; // q_2
776 t_8[3] = ioc + i; // X
777
778 t_9[0] = iov + v_c[4]; // v_4
779 t_9[1] = iov + v_c[5]; // v_5
780 t_9[2] = ioq + q_c[2]; // q_2
781 t_9[3] = ioc + i; // X
782
783 t_10[0] = iov + v_c[0]; // v_0
784 t_10[1] = iov + v_c[4]; // v_4
785 t_10[2] = ioq + q_c[2]; // q_2
786 t_10[3] = ioc + i; // X
787
788 t_11[0] = iov + v_c[5]; // v_5
789 t_11[1] = iov + v_c[1]; // v_1
790 t_11[2] = ioq + q_c[2]; // q_2
791 t_11[3] = ioc + i; // X
792
793 // face 3
794 t_12[0] = iov + v_c[2]; // v_2
795 t_12[1] = iov + v_c[3]; // v_3
796 t_12[2] = ioq + q_c[3]; // q_3
797 t_12[3] = ioc + i; // X
798
799 t_13[0] = iov + v_c[7]; // v_7
800 t_13[1] = iov + v_c[6]; // v_6
801 t_13[2] = ioq + q_c[3]; // q_3
802 t_13[3] = ioc + i; // X
803
804 t_14[0] = iov + v_c[6]; // v_6
805 t_14[1] = iov + v_c[2]; // v_2
806 t_14[2] = ioq + q_c[3]; // q_3
807 t_14[3] = ioc + i; // X
808
809 t_15[0] = iov + v_c[3]; // v_3
810 t_15[1] = iov + v_c[7]; // v_7
811 t_15[2] = ioq + q_c[3]; // q_3
812 t_15[3] = ioc + i; // X
813
814 // face 4
815 t_16[0] = iov + v_c[0]; // v_0
816 t_16[1] = iov + v_c[2]; // v_2
817 t_16[2] = ioq + q_c[4]; // q_4
818 t_16[3] = ioc + i; // X
819
820 t_17[0] = iov + v_c[6]; // v_6
821 t_17[1] = iov + v_c[4]; // v_4
822 t_17[2] = ioq + q_c[4]; // q_4
823 t_17[3] = ioc + i; // X
824
825 t_18[0] = iov + v_c[4]; // v_4
826 t_18[1] = iov + v_c[0]; // v_0
827 t_18[2] = ioq + q_c[4]; // q_4
828 t_18[3] = ioc + i; // X
829
830 t_19[0] = iov + v_c[2]; // v_2
831 t_19[1] = iov + v_c[6]; // v_6
832 t_19[2] = ioq + q_c[4]; // q_4
833 t_19[3] = ioc + i; // X
834
835 // face 5
836 t_20[0] = iov + v_c[3]; // v_3
837 t_20[1] = iov + v_c[1]; // v_1
838 t_20[2] = ioq + q_c[5]; // q_5
839 t_20[3] = ioc + i; // X
840
841 t_21[0] = iov + v_c[5]; // v_5
842 t_21[1] = iov + v_c[7]; // v_7
843 t_21[2] = ioq + q_c[5]; // q_5
844 t_21[3] = ioc + i; // X
845
846 t_22[0] = iov + v_c[1]; // v_1
847 t_22[1] = iov + v_c[5]; // v_5
848 t_22[2] = ioq + q_c[5]; // q_5
849 t_22[3] = ioc + i; // X
850
851 t_23[0] = iov + v_c[7]; // v_7
852 t_23[1] = iov + v_c[3]; // v_3
853 t_23[2] = ioq + q_c[5]; // q_5
854 t_23[3] = ioc + i; // X
855 }
856
857 // return triangle count
858 return 24*num_cubes;
859 }
860 };
861
862 /* *************************************************************************************** */
863 /* *************************************************************************************** */
864 /* ***************************** Simplex --> Hypercube *********************************** */
865 /* *************************************************************************************** */
866 /* *************************************************************************************** */
867
873 template<>
874 struct ShapeConvertIndex<Shape::Hypercube<1>, 1>
875 {
876 static constexpr int shape_dim = 1;
877 static constexpr int cell_dim = 1;
878 typedef Shape::Hypercube<shape_dim> ShapeType;
879 typedef Shape::Simplex<shape_dim> OtherShapeType;
880 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
881 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
882 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
883
884 static Index refine(
885 IndexSetType& index_set_out,
886 const Index offset,
887 const Index index_offsets[],
888 const IndexSetHolderType& index_set_holder_in)
889 {
890 // fetch vertex index offsets
891 const Index iov = index_offsets[0];
892 const Index ioe = index_offsets[1];
893
894 // typedef index set vector reference for output index set
895 typedef IndexSetType::IndexTupleType IndexTupleType;
896
897 // typedef index set types
898 typedef IndexSet<2> IndexSetTypeIn;
899
900 // typedef index vector type
901 typedef IndexSetTypeIn::IndexTupleType IndexTupleTypeIn;
902
903 // fetch the vertices-at-edge index set
904 const IndexSetTypeIn& index_set_in = index_set_holder_in.get_index_set<1,0>();
905
906 // fetch number of coarse mesh edges
907 const Index num_edges = index_set_in.get_num_entities();
908
909 // loop over all coarse mesh edges
910 FEAT_PRAGMA_OMP(parallel for)
911 for(Index i = 0; i < num_edges; ++i)
912 {
913 // fetch coarse mesh vertices-at-edge index vector
914 const IndexTupleTypeIn& v_e = index_set_in[i];
915
916 // fetch fine mesh vertices-at-edge index vectors
917 IndexTupleType& e_0 = index_set_out[offset + 2*i+0];
918 IndexTupleType& e_1 = index_set_out[offset + 2*i+1];
919
920 e_0[0] = iov + v_e[0];
921 e_0[1] = ioe + i;
922
923 e_1[0] = ioe + i;
924 e_1[1] = iov + v_e[1];
925 }
926
927 // return inner edge count
928 return 2*num_edges;
929 }
930 };
931
937 template<>
938 struct ShapeConvertIndex<Shape::Hypercube<2>, 1>
939 {
940 static constexpr int shape_dim = 2;
941 static constexpr int cell_dim = 1;
942 typedef Shape::Hypercube<shape_dim> ShapeType;
943 typedef Shape::Simplex<shape_dim> OtherShapeType;
944 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
945 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
946 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
947
948 static Index refine(
949 IndexSetType& index_set_out,
950 const Index offset,
951 const Index index_offsets[],
952 const IndexSetHolderType& index_set_holder_in)
953 {
954 // fetch vertex index offsets
955 const Index ioe = index_offsets[1];
956 const Index iot = index_offsets[2];
957
958 // typedef index set vector reference for output index set
959 typedef IndexSetType::IndexTupleType IndexTupleType;
960
961 // typedef index set types
962 typedef IndexSet<3> IndexSetTypeIn;
963
964 // typedef index vector type
965 typedef IndexSetTypeIn::IndexTupleType IndexTupleTypeIn;
966
967 // fetch the edge-at-triangle index set
968 const IndexSetTypeIn& index_set_in = index_set_holder_in.get_index_set<2,1>();
969
970 // fetch number of coarse mesh triangle
971 const Index num_trias = index_set_in.get_num_entities();
972
973 // loop over all coarse mesh triangle
974 FEAT_PRAGMA_OMP(parallel for)
975 for(Index i = 0; i < num_trias; ++i)
976 {
977 // fetch coarse mesh edge-at-triangle index vector
978 const IndexTupleTypeIn& e_t = index_set_in[i];
979
980 // fetch fine mesh vertices-at-edge index vectors
981 IndexTupleType& e_0 = index_set_out[offset + 3*i + 0];
982 IndexTupleType& e_1 = index_set_out[offset + 3*i + 1];
983 IndexTupleType& e_2 = index_set_out[offset + 3*i + 2];
984
985 // end with X
986 e_0[0] = ioe + e_t[0]; // v_0
987 e_0[1] = iot + i; // X
988
989 e_1[0] = ioe + e_t[1]; // v_1
990 e_1[1] = iot + i; // X
991
992 e_2[0] = ioe + e_t[2]; // v_2
993 e_2[1] = iot + i; // X
994 }
995
996 // return inner edge count
997 return 3*num_trias;
998 }
999 };
1000
1006 template<>
1007 struct ShapeConvertIndex<Shape::Hypercube<3>, 1>
1008 {
1009 static constexpr int shape_dim = 3;
1010 static constexpr int cell_dim = 1;
1011 typedef Shape::Hypercube<shape_dim> ShapeType;
1012 typedef Shape::Simplex<shape_dim> OtherShapeType;
1013 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1014 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
1015 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
1016
1017 static Index refine(
1018 IndexSetType& index_set_out,
1019 const Index offset,
1020 const Index index_offsets[],
1021 const IndexSetHolderType& index_set_holder_in)
1022 {
1023 // fetch vertex index offsets
1024 const Index iot = index_offsets[2];
1025 const Index iotet = index_offsets[3];
1026
1027 // typedef index set vector reference for output index set
1028 typedef IndexSetType::IndexTupleType IndexTupleType;
1029
1030 // typedef index set types
1031 typedef IndexSet<4> IndexSetTypeIn_t_tet;
1032
1033 // typedef index vector type
1034 typedef IndexSetTypeIn_t_tet::IndexTupleType IndexTupleType_t_tet;
1035
1036 // fetch the triangle-at-tetra
1037 const IndexSetTypeIn_t_tet& index_set_in_t_tet = index_set_holder_in.get_index_set<3,2>();
1038
1039 // fetch number of coarse mesh tetras
1040 const Index num_tetras = index_set_in_t_tet.get_num_entities();
1041
1042 // loop over all coarse mesh tetras
1043 FEAT_PRAGMA_OMP(parallel for)
1044 for(Index i = 0; i < num_tetras; ++i)
1045 {
1046 // fetch coarse mesh triangle-at-tetra index vector
1047 const IndexTupleType_t_tet& t_tet = index_set_in_t_tet[i];
1048
1049 // fetch new mesh vertices-at-edge index vectors
1050 IndexTupleType& e_0 = index_set_out[offset + 4*i + 0];
1051 IndexTupleType& e_1 = index_set_out[offset + 4*i + 1];
1052 IndexTupleType& e_2 = index_set_out[offset + 4*i + 2];
1053 IndexTupleType& e_3 = index_set_out[offset + 4*i + 3];
1054
1055
1056 // end with X
1057 e_0[0] = iot + t_tet[0]; // v_0
1058 e_0[1] = iotet + i; // X
1059
1060 e_1[0] = iot + t_tet[1]; // v_1
1061 e_1[1] = iotet + i; // X
1062
1063 e_2[0] = iot + t_tet[2]; // v_2
1064 e_2[1] = iotet + i; // X
1065
1066 e_3[0] = iot + t_tet[3]; // v_3
1067 e_3[1] = iotet + i; // X
1068
1069 }
1070
1071 // return inner edge count
1072 return 4*num_tetras;
1073 }
1074 };
1075
1081 template<>
1082 struct ShapeConvertIndex<Shape::Hypercube<2>, 2>
1083 {
1084 static constexpr int shape_dim = 2;
1085 static constexpr int cell_dim = 2;
1086 typedef Shape::Hypercube<shape_dim> ShapeType;
1087 typedef Shape::Simplex<shape_dim> OtherShapeType;
1088 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1089 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
1090 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
1091
1092 static Index refine(
1093 IndexSetType& index_set_out,
1094 const Index offset,
1095 const Index index_offsets[],
1096 const IndexSetHolderType& index_set_holder_in)
1097 {
1098 // fetch vertex index offsets
1099 const Index iov = index_offsets[0];
1100 const Index ioe = index_offsets[1];
1101 const Index iot = index_offsets[2];
1102
1103 // typedef index set vector reference for output index set
1104 typedef IndexSetType::IndexTupleType IndexTupleType;
1105
1106 // typedef index set types
1107 typedef IndexSet<3> IndexSetTypeIn_v_t;
1108 typedef IndexSet<3> IndexSetTypeIn_e_t;
1109
1110 // typedef index vector type
1111 typedef IndexSetTypeIn_v_t::IndexTupleType IndexTupleType_v_t;
1112 typedef IndexSetTypeIn_e_t::IndexTupleType IndexTupleType_e_t;
1113
1114 // fetch the vertices-at-triangle and edge-at-triangle index set
1115 const IndexSetTypeIn_v_t& index_set_in_v_t = index_set_holder_in.get_index_set<2,0>();
1116 const IndexSetTypeIn_e_t& index_set_in_e_t = index_set_holder_in.get_index_set<2,1>();
1117
1118 // fetch number of coarse mesh triangle
1119 const Index num_trias = index_set_in_e_t.get_num_entities();
1120
1121 // loop over all coarse mesh triangles
1122 FEAT_PRAGMA_OMP(parallel for)
1123 for(Index i = 0; i < num_trias; ++i)
1124 {
1125 // fetch coarse mesh vertices-at-triangle and edge-at-triangle index vector
1126 const IndexTupleType_v_t& v_t = index_set_in_v_t[i];
1127 const IndexTupleType_e_t& e_t = index_set_in_e_t[i];
1128
1129 // fetch fine mesh vertices-at-quads index vectors
1130 IndexTupleType& q_0 = index_set_out[offset + 3*i + 0];
1131 IndexTupleType& q_1 = index_set_out[offset + 3*i + 1];
1132 IndexTupleType& q_2 = index_set_out[offset + 3*i + 2];
1133
1134 // start with X and end with the corresponding corner
1135 // follow the convention in between (wiki)
1136 q_0[0] = iot + i; // X
1137 q_0[1] = ioe + e_t[1]; // e_1
1138 q_0[2] = ioe + e_t[2]; // e_2
1139 q_0[3] = iov + v_t[0]; // v_0
1140
1141 q_1[0] = iot + i; // X
1142 q_1[1] = ioe + e_t[2]; // e_2
1143 q_1[2] = ioe + e_t[0]; // e_0
1144 q_1[3] = iov + v_t[1]; // v_1
1145
1146 q_2[0] = iot + i; // X
1147 q_2[1] = ioe + e_t[0]; // e_0
1148 q_2[2] = ioe + e_t[1]; // e_1
1149 q_2[3] = iov + v_t[2]; // v_2
1150
1151 }
1152
1153 // return quads count
1154 return 3*num_trias;
1155 }
1156 };
1157
1163 template<>
1164 struct ShapeConvertIndex<Shape::Hypercube<3>, 2>
1165 {
1166 static constexpr int shape_dim = 3;
1167 static constexpr int cell_dim = 2;
1168 typedef Shape::Hypercube<shape_dim> ShapeType;
1169 typedef Shape::Simplex<shape_dim> OtherShapeType;
1170 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1171 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
1172 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
1173
1174 static Index refine(
1175 IndexSetType& index_set_out,
1176 const Index offset,
1177 const Index index_offsets[],
1178 const IndexSetHolderType& index_set_holder_in)
1179 {
1180 // fetch vertex index offsets
1181 const Index ioe = index_offsets[1];
1182 const Index iot = index_offsets[2];
1183 const Index iotet = index_offsets[3];
1184
1185 // typedef index set vector reference for output index set
1186 typedef IndexSetType::IndexTupleType IndexTupleType;
1187
1188 // typedef index set types
1189 typedef IndexSet<6> IndexSetTypeIn_e_tet;
1190 typedef IndexSet<4> IndexSetTypeIn_t_tet;
1191
1192 // typedef index vector type
1193 typedef IndexSetTypeIn_e_tet::IndexTupleType IndexTupleType_e_tet;
1194 typedef IndexSetTypeIn_t_tet::IndexTupleType IndexTupleType_t_tet;
1195
1196 // fetch the edge-at-tetra and triangle-at-tetra index set
1197 const IndexSetTypeIn_e_tet& index_set_in_e_tet = index_set_holder_in.get_index_set<3,1>();
1198 const IndexSetTypeIn_t_tet& index_set_in_t_tet = index_set_holder_in.get_index_set<3,2>();
1199
1200 // fetch number of coarse mesh tetras
1201 const Index num_tetras = index_set_in_e_tet.get_num_entities();
1202
1203 // loop over all coarse mesh tetras
1204 FEAT_PRAGMA_OMP(parallel for)
1205 for(Index i = 0; i < num_tetras; ++i)
1206 {
1207 // fetch coarse mesh edge-at-tetra and triangle-at-tetra index vector
1208 const IndexTupleType_e_tet& e_tet = index_set_in_e_tet[i];
1209 const IndexTupleType_t_tet& t_tet = index_set_in_t_tet[i];
1210
1211 // fetch fine mesh vertices-at-quad index vectors
1212 IndexTupleType& q_0 = index_set_out[offset + 6*i + 0];
1213 IndexTupleType& q_1 = index_set_out[offset + 6*i + 1];
1214 IndexTupleType& q_2 = index_set_out[offset + 6*i + 2];
1215 IndexTupleType& q_3 = index_set_out[offset + 6*i + 3];
1216 IndexTupleType& q_4 = index_set_out[offset + 6*i + 4];
1217 IndexTupleType& q_5 = index_set_out[offset + 6*i + 5];
1218
1219 // start with X move to the face with lower index
1220 // and end up with edge vertex index
1221 q_0[0] = iotet + i; // X
1222 q_0[1] = iot + t_tet[2]; // t_2
1223 q_0[2] = iot + t_tet[3]; // t_3
1224 q_0[3] = ioe + e_tet[0]; // e_0
1225
1226 q_1[0] = iotet + i; // X
1227 q_1[1] = iot + t_tet[1]; // t_1
1228 q_1[2] = iot + t_tet[3]; // t_3
1229 q_1[3] = ioe + e_tet[1]; // e_1
1230
1231 q_2[0] = iotet + i; // X
1232 q_2[1] = iot + t_tet[1]; // t_1
1233 q_2[2] = iot + t_tet[2]; // t_2
1234 q_2[3] = ioe + e_tet[2]; // e_2
1235
1236 q_3[0] = iotet + i; // X
1237 q_3[1] = iot + t_tet[0]; // t_0
1238 q_3[2] = iot + t_tet[3]; // t_3
1239 q_3[3] = ioe + e_tet[3]; // e_3
1240
1241 q_4[0] = iotet + i; // X
1242 q_4[1] = iot + t_tet[0]; // t_0
1243 q_4[2] = iot + t_tet[2]; // t_2
1244 q_4[3] = ioe + e_tet[4]; // e_4
1245
1246 q_5[0] = iotet + i; // X
1247 q_5[1] = iot + t_tet[0]; // t_0
1248 q_5[2] = iot + t_tet[1]; // t_1
1249 q_5[3] = ioe + e_tet[5]; // e_5
1250
1251 }
1252
1253 // return quad count
1254 return 6*num_tetras;
1255 }
1256 };
1257
1264 template<>
1265 struct ShapeConvertIndex<Shape::Hypercube<3>, 3>
1266 {
1267 static constexpr int shape_dim = 3;
1268 static constexpr int cell_dim = 3;
1269 typedef Shape::Hypercube<shape_dim> ShapeType;
1270 typedef Shape::Simplex<shape_dim> OtherShapeType;
1271 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1272 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
1273 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
1274
1275 static Index refine(
1276 IndexSetType& index_set_out,
1277 const Index offset,
1278 const Index index_offsets[],
1279 const IndexSetHolderType& index_set_holder_in)
1280 {
1281 // fetch vertex index offsets
1282 const Index iov = index_offsets[0];
1283 const Index ioe = index_offsets[1];
1284 const Index iot = index_offsets[2];
1285 const Index iotet = index_offsets[3];
1286
1287 // typedef index set vector reference for output index set
1288 typedef IndexSetType::IndexTupleType IndexTupleType;
1289
1290 // typedef index set types
1291 typedef IndexSet<4> IndexSetTypeIn_v_tet;
1292 typedef IndexSet<6> IndexSetTypeIn_e_tet;
1293 typedef IndexSet<4> IndexSetTypeIn_t_tet;
1294
1295 // typedef index vector type
1296 typedef IndexSetTypeIn_v_tet::IndexTupleType IndexTupleType_v_tet;
1297 typedef IndexSetTypeIn_e_tet::IndexTupleType IndexTupleType_e_tet;
1298 typedef IndexSetTypeIn_t_tet::IndexTupleType IndexTupleType_t_tet;
1299
1300 // fetch the vertices-at-tetra, edge-at-tetra and triangle-at-tetra index set
1301 const IndexSetTypeIn_v_tet& index_set_in_v_tet = index_set_holder_in.get_index_set<3,0>();
1302 const IndexSetTypeIn_e_tet& index_set_in_e_tet = index_set_holder_in.get_index_set<3,1>();
1303 const IndexSetTypeIn_t_tet& index_set_in_t_tet = index_set_holder_in.get_index_set<3,2>();
1304
1305 // fetch number of coarse mesh tetras
1306 const Index num_tetras = index_set_in_e_tet.get_num_entities();
1307
1308 // loop over all coarse mesh tetras
1309 FEAT_PRAGMA_OMP(parallel for)
1310 for(Index i = 0; i < num_tetras; ++i)
1311 {
1312 // fetch coarse mesh vertices-at-tetra, edge-at-tetra and triangle-at-tetra index vector
1313 const IndexTupleType_v_tet& v_tet = index_set_in_v_tet[i];
1314 const IndexTupleType_e_tet& e_tet = index_set_in_e_tet[i];
1315 const IndexTupleType_t_tet& t_tet = index_set_in_t_tet[i];
1316
1317 // fetch fine mesh vertices-at-hexa index vectors
1318 IndexTupleType& c_0 = index_set_out[offset + 4*i + 0];
1319 IndexTupleType& c_1 = index_set_out[offset + 4*i + 1];
1320 IndexTupleType& c_2 = index_set_out[offset + 4*i + 2];
1321 IndexTupleType& c_3 = index_set_out[offset + 4*i + 3];
1322
1323 // start with X. move the face with lowest index,
1324 // move to the next lowest face index and proceed according to
1325 // the conventions (wiki) and end with the corresponding corner
1326 c_0[0] = iotet + i; // X
1327 c_0[1] = iot + t_tet[1]; // t_1
1328 c_0[2] = iot + t_tet[2]; // t_2
1329 c_0[3] = ioe + e_tet[2]; // e_2
1330 c_0[4] = iot + t_tet[3]; // t_3
1331 c_0[5] = ioe + e_tet[1]; // e_1
1332 c_0[6] = ioe + e_tet[0]; // e_0
1333 c_0[7] = iov + v_tet[0]; // v_0
1334
1335 c_1[0] = iotet + i; // X
1336 c_1[1] = iot + t_tet[0]; // t_0
1337 c_1[2] = iot + t_tet[2]; // t_2
1338 c_1[3] = ioe + e_tet[4]; // e_4
1339 c_1[4] = iot + t_tet[3]; // t_3
1340 c_1[5] = ioe + e_tet[3]; // e_3
1341 c_1[6] = ioe + e_tet[0]; // e_0
1342 c_1[7] = iov + v_tet[1]; // v_1
1343
1344 c_2[0] = iotet + i; // X
1345 c_2[1] = iot + t_tet[0]; // t_0
1346 c_2[2] = iot + t_tet[1]; // t_1
1347 c_2[3] = ioe + e_tet[5]; // e_5
1348 c_2[4] = iot + t_tet[3]; // t_3
1349 c_2[5] = ioe + e_tet[3]; // e_3
1350 c_2[6] = ioe + e_tet[1]; // e_1
1351 c_2[7] = iov + v_tet[2]; // v_2
1352
1353 c_3[0] = iotet + i; // X
1354 c_3[1] = iot + t_tet[0]; // t_0
1355 c_3[2] = iot + t_tet[1]; // t_1
1356 c_3[3] = ioe + e_tet[5]; // e_5
1357 c_3[4] = iot + t_tet[2]; // t_2
1358 c_3[5] = ioe + e_tet[4]; // e_4
1359 c_3[6] = ioe + e_tet[2]; // e_2
1360 c_3[7] = iov + v_tet[3]; // v_3
1361 }
1362
1363 // return hexa count
1364 return 4*num_tetras;
1365 }
1366 };
1367 /* *************************************************************************************** */
1368 /* *************************************************************************************** */
1369 /* *************************************************************************************** */
1370 /* *************************************************************************************** */
1371 /* *************************************************************************************** */
1372
1373 template<
1374 typename Shape_,
1375 int cell_dim_,
1376 // The following "dummy" argument is necessary for partial specialization;
1377 // shape_dim_ *must* coincide with Shape_::dimension !!!
1378 int shape_dim_ = Shape_::dimension>
1379 struct ShapeConvertIndexShapeWrapper
1380 {
1381 static_assert(shape_dim_ == Shape_::dimension, "invalid shape dimension");
1382 // the case shape_dim_ = cell_dim_ is specialized below
1383 static_assert(shape_dim_ > cell_dim_, "invalid cell dimension");
1384 static_assert(cell_dim_ > 0, "invalid cell dimension");
1385
1386 typedef Shape_ ShapeType;
1387 typedef typename OtherShape<Shape_>::Type OtherShapeType;
1388 typedef typename Shape::FaceTraits<ShapeType, cell_dim_>::ShapeType CellType;
1389 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
1390 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
1391
1392 static String name()
1393 {
1394 return "ShapeConvertIndexShapeWrapper<" + Shape_::name() + "," + stringify(cell_dim_) + ","
1395 + stringify(shape_dim_) + ">";
1396 }
1397
1398 static Index refine(
1399 IndexSetType& index_set_out,
1400 const Index index_offsets[],
1401 const IndexSetHolderType& index_set_holder_in)
1402 {
1403 typedef typename Shape::FaceTraits<ShapeType, shape_dim_ - 1>::ShapeType FacetType;
1404
1405 // refine facets
1406 Index offset = ShapeConvertIndexShapeWrapper<FacetType, cell_dim_>
1407 ::refine(index_set_out, index_offsets, index_set_holder_in);
1408
1409 // call index refiner and return new offset
1410 Index num_faces = ShapeConvertIndex<Shape_, cell_dim_>
1411 ::refine(index_set_out, offset, index_offsets, index_set_holder_in);
1412
1413 // validate number of created indices
1414 Index num_shapes = index_set_holder_in.template get_index_set<shape_dim_, 0>().get_num_entities();
1415 XASSERTM(num_faces == (ShapeConvertTraits<OtherShapeType, cell_dim_>::count * num_shapes),
1416 "ShapeConvertIndex output does not match ShapeConvertTraits prediction");
1417
1418 // calculate new offset and return
1419 return offset + num_faces;
1420 }
1421 };
1422
1423 template<
1424 typename Shape_,
1425 int cell_dim_>
1426 struct ShapeConvertIndexShapeWrapper<Shape_, cell_dim_, cell_dim_>
1427 {
1428 static_assert(cell_dim_ == Shape_::dimension, "invalid shape dimension");
1429 static_assert(cell_dim_ > 0, "invalid cell dimension");
1430
1431 typedef Shape_ ShapeType;
1432 typedef typename OtherShape<Shape_>::Type OtherShapeType;
1433 typedef typename Shape::FaceTraits<ShapeType, cell_dim_>::ShapeType CellType;
1434 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
1435 typedef IndexSetHolder<OtherShapeType> IndexSetHolderType;
1436
1437 static String name()
1438 {
1439 return "ShapeConvertIndexShapeWrapper<" + Shape_::name() + "," + stringify(cell_dim_) + ","
1440 + "," + stringify(cell_dim_) + ">";
1441 }
1442
1443 static Index refine(
1444 IndexSetType& index_set_out,
1445 const Index index_offsets[],
1446 const IndexSetHolderType& index_set_holder_in)
1447 {
1448 // call index refiner
1449 Index num_faces = ShapeConvertIndex<Shape_, cell_dim_>
1450 ::refine(index_set_out, 0, index_offsets, index_set_holder_in);
1451
1452 // validate number of created indices
1453 Index num_shapes = index_set_holder_in.template get_index_set<cell_dim_, 0>().get_num_entities();
1454 XASSERTM(num_faces == (ShapeConvertTraits<OtherShapeType, cell_dim_>::count * num_shapes),
1455 "ShapeConvertIndex output does not match ShapeConvertTraits prediction");
1456
1457 // return offset
1458 return num_faces;
1459 }
1460 };
1461
1462 /* ********************************************************************************************************* */
1463
1464 template<
1465 typename Shape_,
1466 int cell_dim_ = Shape_::dimension>
1467 struct ShapeConvertIndexWrapper
1468 {
1469 // the case cell_dim_ = 1 is handled by the partial specialization below
1470 static_assert(cell_dim_ > 1, "invalid cell dimension");
1471 static_assert(cell_dim_ <= Shape_::dimension, "invalid cell dimension");
1472
1473 typedef Shape_ ShapeType;
1474 typedef typename OtherShape<Shape_>::Type OtherShapeType;
1475 typedef typename Shape::FaceTraits<ShapeType, cell_dim_>::ShapeType CellType;
1476 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
1477 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
1478 typedef IndexSetHolder<OtherShapeType> IndexSetHolderTypeIn;
1479
1480 static String name()
1481 {
1482 return "ShapeConvertIndexWrapper<" + Shape_::name() + "," + stringify(cell_dim_) + ">";
1483 }
1484
1485 static void refine(
1486 IndexSetHolderType& index_set_holder_out,
1487 const Index num_entities[],
1488 const IndexSetHolderTypeIn& index_set_holder_in)
1489 {
1490 // recursive call of IndexRefineWrapper
1491 ShapeConvertIndexWrapper<Shape_, cell_dim_ - 1>
1492 ::refine(index_set_holder_out, num_entities, index_set_holder_in);
1493
1494 // fetch output index set
1495 IndexSetType& index_set_out = index_set_holder_out.template get_index_set<cell_dim_, 0>();
1496
1497 // calculate index offsets
1498 Index index_offsets[Shape_::dimension+1];
1499 EntityCounter<ShapeConvertTraits,OtherShapeType, 0>::offset(index_offsets, num_entities);
1500
1501 // call face wrapper
1502 ShapeConvertIndexShapeWrapper<Shape_, cell_dim_>
1503 ::refine(index_set_out, index_offsets, index_set_holder_in);
1504 }
1505 };
1506
1507 template<typename Shape_>
1508 struct ShapeConvertIndexWrapper<Shape_, 1>
1509 {
1510 static_assert(1 <= Shape_::dimension, "invalid shape dimension");
1511
1512 typedef Shape_ ShapeType;
1513 typedef typename OtherShape<Shape_>::Type OtherShapeType;
1514 typedef typename Shape::FaceTraits<ShapeType, 1>::ShapeType CellType;
1515 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
1516 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
1517 typedef IndexSetHolder<OtherShapeType> IndexSetHolderTypeIn;
1518
1519 static String name()
1520 {
1521 return "ShapeConvertIndexWrapper<" + Shape_::name() + ",1>";
1522 }
1523
1524 static void refine(
1525 IndexSetHolderType& index_set_holder_out,
1526 const Index num_entities[],
1527 const IndexSetHolderTypeIn& index_set_holder_in)
1528 {
1529 // fetch output index set
1530 IndexSetType& index_set_out = index_set_holder_out.template get_index_set<1, 0>();
1531
1532 // calculate index offsets
1533 Index index_offsets[Shape_::dimension+1];
1534 EntityCounter<ShapeConvertTraits,OtherShapeType, 0>::offset(index_offsets, num_entities);
1535
1536 // call face wrapper
1537 ShapeConvertIndexShapeWrapper<Shape_, 1>::refine(index_set_out, index_offsets, index_set_holder_in);
1538 }
1539 };
1540 } // namespace Intern
1542 } // namespace Geometry
1543} // namespace FEAT
XASSERTM
#define XASSERTM(expr, msg)
Assertion macro definition with custom message.
Definition: assertion.hpp:263
FEAT
FEAT namespace.
Definition: adjactor.hpp:12
FEAT::stringify
String stringify(const T_ &item)
Converts an item into a String.
Definition: string.hpp:993
FEAT::Index
std::uint64_t Index
Index data type.
Definition: base_header.hpp:122