FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
standard_index_refiner.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/entity_counter.hpp>
11#include <kernel/geometry/intern/standard_refinement_traits.hpp>
12#include <kernel/geometry/intern/sub_index_mapping.hpp>
13
14namespace FEAT
15{
16 namespace Geometry
17 {
19 namespace Intern
20 {
21 template<
22 typename Shape_,
23 int cell_dim_,
24 int face_dim_>
25 struct StandardIndexRefiner;
26
32 template<>
33 struct StandardIndexRefiner<Shape::Simplex<1>, 1, 0>
34 {
35 static constexpr int shape_dim = 1;
36 static constexpr int cell_dim = 1; // second template parameter
37 static constexpr int face_dim = 0; // third template parameter
38
39 typedef Shape::Simplex<shape_dim> ShapeType;
40 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
41 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
42 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
43
44 static Index refine(
45 IndexSetType& index_set_out,
46 const Index offset,
47 const Index index_offsets[],
48 const IndexSetHolderType& index_set_holder_in)
49 {
50 // fetch vertex index offsets
51 const Index iov = index_offsets[0];
52 const Index ioe = index_offsets[1];
53
54 // typedef index set vector reference for output index set
55 typedef IndexSetType::IndexTupleType IndexTupleType;
56
57 // typedef index set type; an edge has 2 vertices, so we need an IndexSet<2>
58 typedef IndexSet<2> IndexSetTypeEV;
59
60 // typedef index vector reference type
61 typedef IndexSetTypeEV::IndexTupleType IndexTupleTypeEV;
62
63 // fetch the vertices-at-edge index set
64 const IndexSetTypeEV& index_set_e_v = index_set_holder_in.get_index_set<1,0>();
65
66 // fetch number of coarse mesh edges
67 const Index num_edges = index_set_e_v.get_num_entities();
68
69 // loop over all coarse mesh edges
70 FEAT_PRAGMA_OMP(parallel for)
71 for(Index i = 0; i < num_edges; ++i)
72 {
73 // fetch coarse mesh vertices-at-edge index vector
74 const IndexTupleTypeEV& e_v = index_set_e_v[i];
75
76 // fetch fine mesh vertices-at-edge index vectors
77 IndexTupleType& idx0 = index_set_out[offset + 2*i + 0];
78 IndexTupleType& idx1 = index_set_out[offset + 2*i + 1];
79
80 // calculate fine edge-vertex indices
81 idx0[0] = iov + e_v[0];
82 idx0[1] = ioe + i;
83 idx1[0] = ioe + i;
84 idx1[1] = iov + e_v[1];
85 }
86
87 // return fine edge count
88 return 2*num_edges;
89 }
90 }; // StandardIndexRefiner<Simplex<1>,1,0>
91
92 /* ************************************************************************************* */
93 /* ************************************************************************************* */
94 /* ************************************************************************************* */
95
101 template<>
102 struct StandardIndexRefiner<Shape::Simplex<2>, 1, 0>
103 {
104 static constexpr int shape_dim = 2;
105 static constexpr int cell_dim = 1; // second template parameter
106 static constexpr int face_dim = 0; // third template parameter
107
108 typedef Shape::Simplex<shape_dim> ShapeType;
109 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
110 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
111 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
112
113 static Index refine(
114 IndexSetType& index_set_out,
115 const Index offset,
116 const Index index_offsets[],
117 const IndexSetHolderType& index_set_holder_in)
118 {
119 // fetch vertex index offsets
120 const Index ioe = index_offsets[1];
121
122 // typedef index set vector reference for output index set
123 typedef IndexSetType::IndexTupleType IndexTupleType;
124
125 // typedef index set type; a triangle has 3 edges, so we need an IndexSet<3>
126 typedef IndexSet<3> IndexSetTypeSE;
127
128 // typedef index vector type
129 typedef IndexSetTypeSE::IndexTupleType IndexTupleTypeSE;
130
131 // fetch the edges-at-simplex index set
132 const IndexSetTypeSE& index_set_s_e = index_set_holder_in.get_index_set<2,1>();
133
134 // fetch number of coarse mesh simplices
135 const Index num_simps = index_set_s_e.get_num_entities();
136
137 /*
138 v_2
139 |\
140 | \
141 | \
142 | \
143 | \
144 e_1->-e_0
145 |\ |\
146 | \ | \
147 | ^ v \
148 | \ | \
149 | \| \
150 v_0---e_2----v_1
151 */
152
153 // loop over all coarse mesh simplices
154 FEAT_PRAGMA_OMP(parallel for)
155 for(Index i = 0; i < num_simps; ++i)
156 {
157 // fetch coarse mesh edges-at-triangle index vector
158 const IndexTupleTypeSE& s_e = index_set_s_e[i];
159
160 // fetch fine mesh vertices-at-edge index vectors
161 IndexTupleType& e_0 = index_set_out[offset + 3*i + 0];
162 IndexTupleType& e_1 = index_set_out[offset + 3*i + 1];
163 IndexTupleType& e_2 = index_set_out[offset + 3*i + 2];
164
165 e_0[0] = ioe + s_e[2]; // e_2
166 e_0[1] = ioe + s_e[1]; // e_1
167
168 e_1[0] = ioe + s_e[0]; // e_0
169 e_1[1] = ioe + s_e[2]; // e_2
170
171 e_2[0] = ioe + s_e[1]; // e_1
172 e_2[1] = ioe + s_e[0]; // e_0
173
174 }
175 // return fine edge count
176 return 3*num_simps;
177 }
178 }; // StandardIndexRefiner<Simplex<2>,1,0>
179
185 template<>
186 struct StandardIndexRefiner<Shape::Simplex<2>, 2, 0>
187 {
188 static constexpr int shape_dim = 2;
189 static constexpr int cell_dim = 2; // second template parameter
190 static constexpr int face_dim = 0; // third template parameter
191
192 typedef Shape::Simplex<shape_dim> ShapeType;
193 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
194 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
195 typedef IndexSetHolder<ShapeType> 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 ioe = index_offsets[1];
206
207 // typedef index set vector reference for output index set
208 typedef IndexSetType::IndexTupleType IndexTupleType;
209
210 // typedef index set types; a triangle has 3 vertices and 3 edges, so both need an IndexSet<3>
211 typedef IndexSet<3> IndexSetTypeSV;
212 typedef IndexSet<3> IndexSetTypeSE;
213
214 // typedef index vector type
215 typedef IndexSetTypeSV::IndexTupleType IndexTupleTypeSV;
216 typedef IndexSetTypeSE::IndexTupleType IndexTupleTypeSE;
217
218 // fetch the simplex-vertex index set
219 const IndexSetTypeSV& index_set_s_v = index_set_holder_in.get_index_set<2,0>();
220 const IndexSetTypeSE& index_set_s_e = index_set_holder_in.get_index_set<2,1>();
221
222 // fetch number of coarse mesh simplices
223 const Index num_simps = index_set_s_v.get_num_entities();
224
225 /*
226 v_2
227 |\
228 | \
229 |s \
230 |2 \
231 | \
232 e_1->-e_0
233 |\ s |\
234 | \ 3 | \
235 | ^ v \
236 |s \ | s \
237 |0 \| 1 \
238 v_0---e_2----v_1
239 */
240
241 // loop over all coarse mesh simplices
242 FEAT_PRAGMA_OMP(parallel for)
243 for(Index i = 0; i < num_simps; ++i)
244 {
245 // fetch coarse mesh vertices-at-triangle and edges-at-triangle index vectors
246 const IndexTupleTypeSV& s_v = index_set_s_v[i];
247 const IndexTupleTypeSE& s_e = index_set_s_e[i];
248
249 // fetch fine mesh vertices-at-triangle index vectors
250 IndexTupleType& s_0 = index_set_out[offset + 4*i + 0];
251 IndexTupleType& s_1 = index_set_out[offset + 4*i + 1];
252 IndexTupleType& s_2 = index_set_out[offset + 4*i + 2];
253 IndexTupleType& s_3 = index_set_out[offset + 4*i + 3];
254
255 // calculate fine simplex-vertex indices
256 s_0[0] = iov + s_v[0]; // v_0
257 s_0[1] = ioe + s_e[2]; // e_2
258 s_0[2] = ioe + s_e[1]; // e_1
259
260 s_1[0] = ioe + s_e[2]; // e_2
261 s_1[1] = iov + s_v[1]; // v_1
262 s_1[2] = ioe + s_e[0]; // e_0
263
264 s_2[0] = ioe + s_e[1]; // e_1
265 s_2[1] = ioe + s_e[0]; // e_0
266 s_2[2] = iov + s_v[2]; // s_2
267
268 s_3[0] = ioe + s_e[0]; // e_0
269 s_3[1] = ioe + s_e[1]; // e_1
270 s_3[2] = ioe + s_e[2]; // e_2
271 }
272 // return fine simplex count
273 return 4*num_simps;
274 }
275 }; // StandardIndexRefiner<Simplex<2>,2,0>
276
282 template<>
283 struct StandardIndexRefiner<Shape::Simplex<2>, 2, 1>
284 {
285 static constexpr int shape_dim = 2;
286 static constexpr int cell_dim = 2; // second template parameter
287 static constexpr int face_dim = 1; // third template parameter
288
289 typedef Shape::Simplex<shape_dim> ShapeType;
290 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
291 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
292 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
293
294 static Index refine(
295 IndexSetType& index_set_out,
296 const Index offset,
297 const Index index_offsets[],
298 const IndexSetHolderType& index_set_holder_in)
299 {
300 // fetch edge index offsets
301 const Index ioe = index_offsets[1];
302 const Index ios = index_offsets[2];
303
304 // typedef index set vector reference for output index set
305 typedef IndexSetType::IndexTupleType IndexTupleType;
306
307 // typedef index set types;
308 typedef IndexSet<2> IndexSetTypeEV; // an edge has 2 vertices
309 typedef IndexSet<3> IndexSetTypeSV; // a simplex has 3 vertices
310 typedef IndexSet<3> IndexSetTypeSE; // a simplex has 3 edges
311
312 // typedef index vector type
313 typedef IndexSetTypeSV::IndexTupleType IndexTupleTypeSV;
314 typedef IndexSetTypeSE::IndexTupleType IndexTupleTypeSE;
315
316 // fetch the index sets
317 const IndexSetTypeEV& index_set_e_v = index_set_holder_in.get_index_set<1,0>();
318 const IndexSetTypeSV& index_set_s_v = index_set_holder_in.get_index_set<2,0>();
319 const IndexSetTypeSE& index_set_s_e = index_set_holder_in.get_index_set<2,1>();
320
321 // fetch number of coarse mesh simplices
322 const Index num_simps = index_set_s_v.get_num_entities();
323
324 // typedef the sub-index mapping
325 typedef Intern::SubIndexMapping<ShapeType, face_dim, 0> SubIndexMappingType;
326
327 /*
328 v_2
329 |\
330 | \
331 |s \
332 |2 ^
333 | \
334 e_1->-e_0
335 |\ s |\
336 v \ 3 | \
337 | ^ v \
338 |s \ | s \
339 |0 \| 1 \
340 v_0---e_2-->-v_1
341 */
342
343 // loop over all coarse mesh simplices
344 FEAT_PRAGMA_OMP(parallel for)
345 for(Index i = 0; i < num_simps; ++i)
346 {
347 // fetch coarse mesh index vectors
348 const IndexTupleTypeSV& s_v = index_set_s_v[i];
349 const IndexTupleTypeSE& s_e = index_set_s_e[i];
350
351 // create a sub-index mapping object
352 SubIndexMappingType sim(s_v, s_e, index_set_e_v);
353
354 // fetch fine mesh edges-at-simplex index vectors
355 IndexTupleType& s_0 = index_set_out[offset + 4*i + 0];
356 IndexTupleType& s_1 = index_set_out[offset + 4*i + 1];
357 IndexTupleType& s_2 = index_set_out[offset + 4*i + 2];
358 IndexTupleType& s_3 = index_set_out[offset + 4*i + 3];
359
360 // calculate fine edges-at-simplex indices
361 s_0[0] = ios + 3*i + 0;
362 s_0[1] = ioe + 2*s_e[1] + sim.map(1, 1); // e_1_1
363 s_0[2] = ioe + 2*s_e[2] + sim.map(2, 0); // e_2_0
364
365 s_1[0] = ioe + 2*s_e[0] + sim.map(0, 0); // e_0_0
366 s_1[1] = ios + 3*i + 1;
367 s_1[2] = ioe + 2*s_e[2] + sim.map(2, 1); // e_2_1
368
369 s_2[0] = ioe + 2*s_e[0] + sim.map(0, 1); // e_0_1
370 s_2[1] = ioe + 2*s_e[1] + sim.map(1, 0); // e_1_0
371 s_2[2] = ios + 3*i + 2;
372
373 s_3[0] = ios + 3*i + 0;
374 s_3[1] = ios + 3*i + 1;
375 s_3[2] = ios + 3*i + 2;
376 }
377
378 // return fine triangle count
379 return 4*num_simps;
380 }
381 }; // StandardIndexRefiner<Simplex<2>,2,1>
382
383 /* ************************************************************************************* */
384 /* ************************************************************************************* */
385 /* ************************************************************************************* */
386
392 template<>
393 struct StandardIndexRefiner<Shape::Simplex<3>, 1, 0>
394 {
395 static constexpr int shape_dim = 3;
396 static constexpr int cell_dim = 1; // second template parameter
397 static constexpr int face_dim = 0; // third template parameter
398
399 typedef Shape::Simplex<shape_dim> ShapeType;
400 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
401 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
402 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
403
404 static Index refine(
405 IndexSetType& index_set_out,
406 const Index offset,
407 const Index index_offsets[],
408 const IndexSetHolderType& index_set_holder_in)
409 {
410 // fetch vertex index offsets
411 const Index ioe = index_offsets[1];
412 const Index ios = index_offsets[3];
413
414 // typedef index set vector reference for output index set
415 typedef IndexSetType::IndexTupleType IndexTupleType;
416
417 // typedef index set type (edges at simplices)
418 typedef IndexSet<6> IndexSetTypeSE;
419
420 // typedef index vector type
421 typedef IndexSetTypeSE::IndexTupleType IndexTupleTypeSE;
422
423 // fetch the edge-at-simplex index set
424 const IndexSetTypeSE& index_set_s_e = index_set_holder_in.get_index_set<3,1>();
425
426 // fetch number of coarse mesh simplices
427 const Index num_simps = index_set_s_e.get_num_entities();
428
429 /* v_3
430 /\\
431 / \ \
432 / \ \
433 / \ \
434 / \ \
435 / \ \
436 / \ \
437 / \ x e_5
438 / \ \
439 / \ \
440 / \ \
441 / \ \ v_2
442 e_2 x x e_4 / /
443 / s x \ / /
444 / \ / /
445 / \ / /
446 / / \ /
447 / x \ /
448 / / e_1 \ x e_3
449 / / \ /
450 / / \ /
451 / / \ /
452 / / \ /
453 //_______________________x_____________________\/
454 v_0 e_0 v_1
455
456 orientation:
457 t0: v_1 - v_2 - v_3
458 t1: v_0 - v_2 - v_3
459 t2: v_0 - v_1 - v_3
460 t3: v_0 - v_1 - v_2
461 */
462
463
464 // loop over all coarse mesh simplices
465 FEAT_PRAGMA_OMP(parallel for)
466 for(Index i = 0; i < num_simps; ++i)
467 {
468 // fetch coarse mesh simplex-at-edge index vector
469 const IndexTupleTypeSE& s_e = index_set_s_e[i];
470
471 // fetch fine mesh vertices-at-edge index vectors
472 IndexTupleType& e_0 = index_set_out[offset + 6*i + 0];
473 IndexTupleType& e_1 = index_set_out[offset + 6*i + 1];
474 IndexTupleType& e_2 = index_set_out[offset + 6*i + 2];
475 IndexTupleType& e_3 = index_set_out[offset + 6*i + 3];
476 IndexTupleType& e_4 = index_set_out[offset + 6*i + 4];
477 IndexTupleType& e_5 = index_set_out[offset + 6*i + 5];
478
479 e_0[0] = ioe + s_e[0]; //e_0
480 e_0[1] = ios + i; //s
481
482 e_1[0] = ioe + s_e[1]; //e_1
483 e_1[1] = ios + i; //s
484
485 e_2[0] = ioe + s_e[2]; //e_2
486 e_2[1] = ios + i; //s
487
488 e_3[0] = ioe + s_e[3]; //e_3
489 e_3[1] = ios + i; //s
490
491 e_4[0] = ioe + s_e[4]; //e_4
492 e_4[1] = ios + i; //s
493
494 e_5[0] = ioe + s_e[5]; //e_5
495 e_5[1] = ios + i; //s
496 }
497
498 // return fine edge count
499 return 6*num_simps;
500 }
501 }; // StandardIndexRefiner<Simplex<3>,1,0>
502
508 template<>
509 struct StandardIndexRefiner<Shape::Simplex<3>, 2, 0>
510 {
511 static constexpr int shape_dim = 3;
512 static constexpr int cell_dim = 2; // second template parameter
513 static constexpr int face_dim = 0; // third template parameter
514
515 typedef Shape::Simplex<shape_dim> ShapeType;
516 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
517 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
518 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
519
520 static Index refine(
521 IndexSetType& index_set_out,
522 const Index offset,
523 const Index index_offsets[],
524 const IndexSetHolderType& index_set_holder_in)
525 {
526 // fetch vertex index offsets
527 const Index ioe = index_offsets[1];
528 const Index ios = index_offsets[3];
529
530 // typedef index set vector reference for output index set
531 typedef IndexSetType::IndexTupleType IndexTupleType;
532
533 // typedef index set types
534 typedef IndexSet<4> IndexSetTypeST;
535 typedef IndexSet<6> IndexSetTypeSE;
536
537 // typedef index vector type
538 typedef IndexSetTypeSE::IndexTupleType IndexTupleTypeSE;
539
540 // fetch the index sets
541 const IndexSetTypeST& index_set_s_t = index_set_holder_in.get_index_set<3,2>();
542 const IndexSetTypeSE& index_set_s_e = index_set_holder_in.get_index_set<3,1>();
543
544 // fetch number of coarse mesh simplices
545 const Index num_simps = index_set_s_t.get_num_entities();
546
547 /* v_3
548 /\\
549 / \ \
550 / \ \
551 / \ \
552 / \ \
553 / \ \
554 / \ x e_5
555 / \ \
556 / \ \
557 / \ \
558 / \ \
559 / e_4 \ \ v_2
560 e_2 x\ x / /
561 / \ \ s x \ / /
562 / \ \ \ / /
563 / \ ^ \ / /
564 / \ \ e_1 / \ /
565 / v s_0 x \ /
566 / \ / / \ x e_3
567 / \ / / \ /
568 / / \ ^ \ /
569 / / \ / \ /
570 / / \ / \ /
571 //____________________x________________________\/
572 v_0 e_0 v_1
573
574 orientation:
575 t0: v_1 - v_2 - v_3
576 t1: v_0 - v_2 - v_3
577 t2: v_0 - v_1 - v_3
578 t3: v_0 - v_1 - v_2
579 */
580
581 // loop over all coarse mesh simplices
582 FEAT_PRAGMA_OMP(parallel for)
583 for(Index i = 0; i < num_simps; ++i)
584 {
585 // fetch coarse mesh edges-at-simplex index vectors
586 const IndexTupleTypeSE& s_e = index_set_s_e[i];
587
588 // fetch fine mesh vertices-at-triangle index vectors
589 IndexTupleType& t_0 = index_set_out[offset + 16*i + 0];
590 IndexTupleType& t_1 = index_set_out[offset + 16*i + 1];
591 IndexTupleType& t_2 = index_set_out[offset + 16*i + 2];
592 IndexTupleType& t_3 = index_set_out[offset + 16*i + 3];
593 IndexTupleType& t_4 = index_set_out[offset + 16*i + 4];
594 IndexTupleType& t_5 = index_set_out[offset + 16*i + 5];
595 IndexTupleType& t_6 = index_set_out[offset + 16*i + 6];
596 IndexTupleType& t_7 = index_set_out[offset + 16*i + 7];
597 IndexTupleType& t_8 = index_set_out[offset + 16*i + 8];
598 IndexTupleType& t_9 = index_set_out[offset + 16*i + 9];
599 IndexTupleType& t_10 = index_set_out[offset + 16*i + 10];
600 IndexTupleType& t_11 = index_set_out[offset + 16*i + 11];
601 IndexTupleType& t_12 = index_set_out[offset + 16*i + 12];
602 IndexTupleType& t_13 = index_set_out[offset + 16*i + 13];
603 IndexTupleType& t_14 = index_set_out[offset + 16*i + 14];
604 IndexTupleType& t_15 = index_set_out[offset + 16*i + 15];
605
606 // calculate fine simplex-vertex indices
607
608 // triangles at the corners (e.g. s_0, the orientation is the same as for triangle(local point_index))
609 t_0[0] = ioe + s_e[0]; // e_0
610 t_0[1] = ioe + s_e[1]; // e_1
611 t_0[2] = ioe + s_e[2]; // e_2
612
613 t_1[0] = ioe + s_e[0]; // e_0
614 t_1[1] = ioe + s_e[3]; // e_3
615 t_1[2] = ioe + s_e[4]; // e_4
616
617 t_2[0] = ioe + s_e[1]; // e_1
618 t_2[1] = ioe + s_e[3]; // e_3
619 t_2[2] = ioe + s_e[5]; // e_5
620
621 t_3[0] = ioe + s_e[2]; // e_0
622 t_3[1] = ioe + s_e[4]; // e_1
623 t_3[2] = ioe + s_e[5]; // e_2
624
625 // triangles adjuncted to the center (orientation i < j -> e_i - e_j - s)
626 t_4[0] = ioe + s_e[0]; // e_0
627 t_4[1] = ioe + s_e[1]; // e_1
628 t_4[2] = ios + i; // s
629
630 t_5[0] = ioe + s_e[0]; // e_0
631 t_5[1] = ioe + s_e[2]; // e_2
632 t_5[2] = ios + i; // s
633
634 t_6[0] = ioe + s_e[0]; // e_0
635 t_6[1] = ioe + s_e[3]; // e_3
636 t_6[2] = ios + i; // s
637
638 t_7[0] = ioe + s_e[0]; // e_0
639 t_7[1] = ioe + s_e[4]; // e_4
640 t_7[2] = ios + i; // s
641
642 t_8[0] = ioe + s_e[1]; // e_1
643 t_8[1] = ioe + s_e[2]; // e_2
644 t_8[2] = ios + i; // s
645
646 t_9[0] = ioe + s_e[1]; // e_1
647 t_9[1] = ioe + s_e[3]; // e_3
648 t_9[2] = ios + i; // s
649
650 t_10[0] = ioe + s_e[1]; // e_1
651 t_10[1] = ioe + s_e[5]; // e_5
652 t_10[2] = ios + i; // s
653
654 t_11[0] = ioe + s_e[2]; // e_2
655 t_11[1] = ioe + s_e[4]; // e_4
656 t_11[2] = ios + i; // s
657
658 t_12[0] = ioe + s_e[2]; // e_2
659 t_12[1] = ioe + s_e[5]; // e_5
660 t_12[2] = ios + i; // s
661
662 t_13[0] = ioe + s_e[3]; // e_3
663 t_13[1] = ioe + s_e[4]; // e_4
664 t_13[2] = ios + i; // s
665
666 t_14[0] = ioe + s_e[3]; // e_3
667 t_14[1] = ioe + s_e[5]; // e_5
668 t_14[2] = ios + i; // s
669
670 t_15[0] = ioe + s_e[4]; // e_4
671 t_15[1] = ioe + s_e[5]; // e_5
672 t_15[2] = ios + i; // s
673 }
674
675 // return fine triangle count
676 return 16*num_simps;
677 }
678 }; // StandardIndexRefiner<Simplex<3>,2,0>
679
685 template<>
686 struct StandardIndexRefiner<Shape::Simplex<3>, 3, 0>
687 {
688 static constexpr int shape_dim = 3;
689 static constexpr int cell_dim = 3; // second template parameter
690 static constexpr int face_dim = 0; // third template parameter
691
692 typedef Shape::Simplex<shape_dim> ShapeType;
693 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
694 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
695 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
696
697 static Index refine(
698 IndexSetType& index_set_out,
699 const Index offset,
700 const Index index_offsets[],
701 const IndexSetHolderType& index_set_holder_in)
702 {
703 // fetch vertex index offsets
704 const Index iov = index_offsets[0];
705 const Index ioe = index_offsets[1];
706 const Index ios = index_offsets[3];
707
708 // typedef index set vector reference for output index set
709 typedef IndexSetType::IndexTupleType IndexTupleType;
710
711 // typedef index set types
712 typedef IndexSet<4> IndexSetTypeST;
713 typedef IndexSet<6> IndexSetTypeSE;
714 typedef IndexSet<4> IndexSetTypeSV;
715
716 // typedef index vector type
717 typedef IndexSetTypeSE::IndexTupleType IndexTupleTypeSE;
718 typedef IndexSetTypeSV::IndexTupleType IndexTupleTypeSV;
719
720 // fetch the index sets
721 const IndexSetTypeST& index_set_s_t = index_set_holder_in.get_index_set<3,2>();
722 const IndexSetTypeSE& index_set_s_e = index_set_holder_in.get_index_set<3,1>();
723 const IndexSetTypeSV& index_set_s_v = index_set_holder_in.get_index_set<3,0>();
724
725 // fetch number of coarse mesh simplices
726 const Index num_simps = index_set_s_t.get_num_entities();
727
728 /* v_3
729 /\\
730 / \ \
731 / \ \
732 / \ \
733 / \ \
734 / \ \
735 / \ x e_5
736 / \ \
737 / \ \
738 / \ \
739 / \ \
740 / e_4 \ \ v_2
741 e_2 x\ x / /
742 / \ \ s x \ / /
743 / \ \ \ / /
744 / \ ^ \ / /
745 / \ \ e_1 / \ /
746 / v s_0 x \ /
747 / \ / / \ x e_3
748 / \ / / \ /
749 / / \ ^ \ /
750 / / \ / \ /
751 / / \ / \ /
752 //____________________x________________________\/
753 v_0 e_0 v_1
754
755 orientation:
756 t0: v_1 - v_2 - v_3
757 t1: v_0 - v_2 - v_3
758 t2: v_0 - v_1 - v_3
759 t3: v_0 - v_1 - v_2
760 */
761
762 // loop over all coarse mesh simplices
763 FEAT_PRAGMA_OMP(parallel for)
764 for(Index i = 0; i < num_simps; ++i)
765 {
766 // fetch coarse mesh index vectors
767 const IndexTupleTypeSE& s_e = index_set_s_e[i];
768 const IndexTupleTypeSV& s_v = index_set_s_v[i];
769
770 // fetch fine mesh vertices-at-simplex index vectors
771 IndexTupleType& s_0 = index_set_out[offset + 12*i + 0];
772 IndexTupleType& s_1 = index_set_out[offset + 12*i + 1];
773 IndexTupleType& s_2 = index_set_out[offset + 12*i + 2];
774 IndexTupleType& s_3 = index_set_out[offset + 12*i + 3];
775 IndexTupleType& s_4 = index_set_out[offset + 12*i + 4];
776 IndexTupleType& s_5 = index_set_out[offset + 12*i + 5];
777 IndexTupleType& s_6 = index_set_out[offset + 12*i + 6];
778 IndexTupleType& s_7 = index_set_out[offset + 12*i + 7];
779 IndexTupleType& s_8 = index_set_out[offset + 12*i + 8];
780 IndexTupleType& s_9 = index_set_out[offset + 12*i + 9];
781 IndexTupleType& s_10 = index_set_out[offset + 12*i + 10];
782 IndexTupleType& s_11 = index_set_out[offset + 12*i + 11];
783
784 // calculate fine simplex-vertex indices
785
786 // simplices at the corners
787
788 // v_0
789 s_0[0] = ioe + s_e[0]; // e_0
790 s_0[1] = ioe + s_e[2]; // e_2
791 s_0[2] = ioe + s_e[1]; // e_1
792 s_0[3] = iov + s_v[0]; // v_0
793
794 s_1[0] = ioe + s_e[0]; // e_0
795 s_1[1] = ioe + s_e[1]; // e_1
796 s_1[2] = ioe + s_e[2]; // e_2
797 s_1[3] = ios + i; // s
798
799 // v_1
800 s_2[0] = ioe + s_e[0]; // e_0
801 s_2[1] = ioe + s_e[3]; // e_3
802 s_2[2] = ioe + s_e[4]; // e_4
803 s_2[3] = iov + s_v[1]; // v_0
804
805 s_3[0] = ioe + s_e[0]; // e_0
806 s_3[1] = ioe + s_e[4]; // e_4
807 s_3[2] = ioe + s_e[3]; // e_3
808 s_3[3] = ios + i; // s
809
810 // v_2
811 s_4[0] = ioe + s_e[1]; // e_1
812 s_4[1] = ioe + s_e[5]; // e_5
813 s_4[2] = ioe + s_e[3]; // e_3
814 s_4[3] = iov + s_v[2]; // v_2
815
816 s_5[0] = ioe + s_e[1]; // e_1
817 s_5[1] = ioe + s_e[3]; // e_3
818 s_5[2] = ioe + s_e[5]; // e_5
819 s_5[3] = ios + i; // s
820
821 // v_3
822 s_6[0] = ioe + s_e[2]; // e_2
823 s_6[1] = ioe + s_e[4]; // e_4
824 s_6[2] = ioe + s_e[5]; // e_5
825 s_6[3] = iov + s_v[3]; // v_3
826
827 s_7[0] = ioe + s_e[2]; // e_2
828 s_7[1] = ioe + s_e[5]; // e_5
829 s_7[2] = ioe + s_e[4]; // e_4
830 s_7[3] = ios + i; // s
831
832 // simplices at the coarse mesh faces
833
834 s_8[0] = ioe + s_e[3]; // e_3
835 s_8[1] = ioe + s_e[4]; // e_4
836 s_8[2] = ioe + s_e[5]; // e_5
837 s_8[3] = ios + i; // s
838
839 s_9[0] = ioe + s_e[1]; // e_1
840 s_9[1] = ioe + s_e[5]; // e_5
841 s_9[2] = ioe + s_e[2]; // e_2
842 s_9[3] = ios + i; // s
843
844 s_10[0] = ioe + s_e[0]; // e_0
845 s_10[1] = ioe + s_e[2]; // e_2
846 s_10[2] = ioe + s_e[4]; // e_4
847 s_10[3] = ios + i; // s
848
849 s_11[0] = ioe + s_e[0]; // e_0
850 s_11[1] = ioe + s_e[3]; // e_3
851 s_11[2] = ioe + s_e[1]; // e_1
852 s_11[3] = ios + i; // s
853 }
854 // return fine simplex count
855 return 12*num_simps;
856 }
857 }; // StandardIndexRefiner<Simplex<3>,3,0>
858
864 template<>
865 struct StandardIndexRefiner<Shape::Simplex<3>, 2, 1>
866 {
867 static constexpr int shape_dim = 3;
868 static constexpr int cell_dim = 2; // second template parameter
869 static constexpr int face_dim = 1; // third template parameter
870
871 typedef Shape::Simplex<shape_dim> ShapeType;
872 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
873 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
874 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
875
876 static Index refine(
877 IndexSetType& index_set_out,
878 const Index offset,
879 const Index index_offsets[],
880 const IndexSetHolderType& index_set_holder_in)
881 {
882 // fetch index offsets
883 const Index iot = index_offsets[2];
884 const Index ios = index_offsets[3];
885
886 // typedef index set vector reference for output index set
887 typedef IndexSetType::IndexTupleType IndexTupleType;
888
889 // typedef index set types;
890 typedef IndexSet<3> IndexSetTypeTV;
891 typedef IndexSet<4> IndexSetTypeSV;
892 typedef IndexSet<4> IndexSetTypeST;
893
894 // typedef index vector type
895 typedef IndexSetTypeSV::IndexTupleType IndexTupleTypeSV;
896 typedef IndexSetTypeST::IndexTupleType IndexTupleTypeST;
897
898 // fetch the index sets
899 const IndexSetTypeSV& index_set_s_v = index_set_holder_in.get_index_set<3,0>();
900 const IndexSetTypeST& index_set_s_t = index_set_holder_in.get_index_set<3,2>();
901 const IndexSetTypeTV& index_set_t_v = index_set_holder_in.get_index_set<2,0>();
902
903 // fetch number of coarse mesh simplices
904 const Index num_simps = index_set_s_t.get_num_entities();
905
906 // typedef the sub-index mapping
907 typedef Intern::SubIndexMapping<ShapeType, cell_dim, face_dim> SubIndexMappingType;
908
909 // loop over all coarse mesh simplices
910 FEAT_PRAGMA_OMP(parallel for)
911 for(Index i = 0; i < num_simps; ++i)
912 {
913 // fetch coarse mesh vertices-at-simplex and edges-at-simplex index vectors
914 const IndexTupleTypeSV& s_v = index_set_s_v[i];
915 const IndexTupleTypeST& s_t = index_set_s_t[i];
916
917 // create a sub-index mapping object
918 SubIndexMappingType sim(s_v, s_t, index_set_t_v);
919
920 // fetch fine mesh edges-at-triangle index vectors
921 IndexTupleType& t_0 = index_set_out[offset + 16*i + 0];
922 IndexTupleType& t_1 = index_set_out[offset + 16*i + 1];
923 IndexTupleType& t_2 = index_set_out[offset + 16*i + 2];
924 IndexTupleType& t_3 = index_set_out[offset + 16*i + 3];
925 IndexTupleType& t_4 = index_set_out[offset + 16*i + 4];
926 IndexTupleType& t_5 = index_set_out[offset + 16*i + 5];
927 IndexTupleType& t_6 = index_set_out[offset + 16*i + 6];
928 IndexTupleType& t_7 = index_set_out[offset + 16*i + 7];
929 IndexTupleType& t_8 = index_set_out[offset + 16*i + 8];
930 IndexTupleType& t_9 = index_set_out[offset + 16*i + 9];
931 IndexTupleType& t_10 = index_set_out[offset + 16*i + 10];
932 IndexTupleType& t_11 = index_set_out[offset + 16*i + 11];
933 IndexTupleType& t_12 = index_set_out[offset + 16*i + 12];
934 IndexTupleType& t_13 = index_set_out[offset + 16*i + 13];
935 IndexTupleType& t_14 = index_set_out[offset + 16*i + 14];
936 IndexTupleType& t_15 = index_set_out[offset + 16*i + 15];
937
938 // calculate fine edges-at-triangle indices
939
940 // triangle at the corners
941 t_0[0] = iot + 3*s_t[1] + sim.map(1, 0);
942 t_0[1] = iot + 3*s_t[2] + sim.map(2, 0);
943 t_0[2] = iot + 3*s_t[3] + sim.map(3, 0);
944
945 t_1[0] = iot + 3*s_t[0] + sim.map(0, 0);
946 t_1[1] = iot + 3*s_t[2] + sim.map(2, 1);
947 t_1[2] = iot + 3*s_t[3] + sim.map(3, 1);
948
949 t_2[0] = iot + 3*s_t[0] + sim.map(0, 1);
950 t_2[1] = iot + 3*s_t[1] + sim.map(1, 1);
951 t_2[2] = iot + 3*s_t[3] + sim.map(3, 2);
952
953 t_3[0] = iot + 3*s_t[0] + sim.map(0, 2);
954 t_3[1] = iot + 3*s_t[1] + sim.map(1, 2);
955 t_3[2] = iot + 3*s_t[2] + sim.map(2, 2);
956
957 // triangles at the midpoint
958
959 //e_0 triangles
960 t_4[0] = ios + 6*i + 1;
961 t_4[1] = ios + 6*i + 0;
962 t_4[2] = iot + 3*s_t[3] + sim.map(3, 0);
963
964 t_5[0] = ios + 6*i + 2;
965 t_5[1] = ios + 6*i + 0;
966 t_5[2] = iot + 3*s_t[2] + sim.map(2, 0);
967
968 t_6[0] = ios + 6*i + 3;
969 t_6[1] = ios + 6*i + 0;
970 t_6[2] = iot + 3*s_t[3] + sim.map(3, 1);
971
972 t_7[0] = ios + 6*i + 4;
973 t_7[1] = ios + 6*i + 0;
974 t_7[2] = iot + 3*s_t[2] + sim.map(2, 1);
975
976 //e_1 triangles
977 t_8[0] = ios + 6*i + 2;
978 t_8[1] = ios + 6*i + 1;
979 t_8[2] = iot + 3*s_t[1] + sim.map(1, 0);
980
981 t_9[0] = ios + 6*i + 3;
982 t_9[1] = ios + 6*i + 1;
983 t_9[2] = iot + 3*s_t[3] + sim.map(3, 2);
984
985 t_10[0] = ios + 6*i + 5;
986 t_10[1] = ios + 6*i + 1;
987 t_10[2] = iot + 3*s_t[1] + sim.map(1, 1);
988
989 //e_2 triangles
990 t_11[0] = ios + 6*i + 4;
991 t_11[1] = ios + 6*i + 2;
992 t_11[2] = iot + 3*s_t[2] + sim.map(2, 2);
993
994 t_12[0] = ios + 6*i + 5;
995 t_12[1] = ios + 6*i + 2;
996 t_12[2] = iot + 3*s_t[1] + sim.map(1, 2);
997
998 //e_3 triangles
999 t_13[0] = ios + 6*i + 4;
1000 t_13[1] = ios + 6*i + 3;
1001 t_13[2] = iot + 3*s_t[0] + sim.map(0, 0);
1002
1003 t_14[0] = ios + 6*i + 5;
1004 t_14[1] = ios + 6*i + 3;
1005 t_14[2] = iot + 3*s_t[0] + sim.map(0, 1);
1006
1007 //e_4 triangles
1008 t_15[0] = ios + 6*i + 5;
1009 t_15[1] = ios + 6*i + 4;
1010 t_15[2] = iot + 3*s_t[0] + sim.map(0, 2);
1011
1012 }
1013 // return fine triangle count
1014 return 16*num_simps;
1015 }
1016 }; // StandardIndexRefiner<Simplex<3>,2,1>
1017
1023 template<>
1024 struct StandardIndexRefiner<Shape::Simplex<3>, 3, 1>
1025 {
1026 static constexpr int shape_dim = 3;
1027 static constexpr int cell_dim = 3; // second template parameter
1028 static constexpr int face_dim = 1; // third template parameter
1029
1030 typedef Shape::Simplex<shape_dim> ShapeType;
1031 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1032 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
1033 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
1034
1035 static Index refine(
1036 IndexSetType& index_set_out,
1037 const Index offset,
1038 const Index index_offsets[],
1039 const IndexSetHolderType& index_set_holder_in)
1040 {
1041 // fetch index offsets
1042 const Index ioe = index_offsets[1];
1043 const Index iot = index_offsets[2];
1044 const Index ios = index_offsets[3];
1045
1046 // typedef index set vector reference for output index set
1047 typedef IndexSetType::IndexTupleType IndexTupleType;
1048
1049 // typedef index set types
1050 typedef IndexSet<2> IndexSetTypeEV;
1051 typedef IndexSet<3> IndexSetTypeTV;
1052 typedef IndexSet<4> IndexSetTypeSV;
1053 typedef IndexSet<6> IndexSetTypeSE;
1054 typedef IndexSet<4> IndexSetTypeST;
1055
1056 // typedef index vector type
1057 typedef IndexSetTypeSV::IndexTupleType IndexTupleTypeSV;
1058 typedef IndexSetTypeSE::IndexTupleType IndexTupleTypeSE;
1059 typedef IndexSetTypeST::IndexTupleType IndexTupleTypeST;
1060
1061 // fetch the index sets
1062 const IndexSetTypeEV& index_set_e_v = index_set_holder_in.get_index_set<1,0>();
1063 const IndexSetTypeTV& index_set_t_v = index_set_holder_in.get_index_set<2,0>();
1064 const IndexSetTypeSV& index_set_s_v = index_set_holder_in.get_index_set<3,0>();
1065 const IndexSetTypeSE& index_set_s_e = index_set_holder_in.get_index_set<3,1>();
1066 const IndexSetTypeST& index_set_s_t = index_set_holder_in.get_index_set<3,2>();
1067
1068 // fetch number of coarse mesh simplices
1069 const Index num_simps = index_set_s_t.get_num_entities();
1070
1071 // typedef the sub-index mapping
1072 typedef Intern::SubIndexMapping<ShapeType, 2, 1> SubIndexMappingType;
1073 typedef Intern::SubIndexMapping<ShapeType, face_dim, 0> EdgeIndexMappingType;
1074
1075 // loop over all coarse mesh simplices
1076 FEAT_PRAGMA_OMP(parallel for)
1077 for(Index i = 0; i < num_simps; ++i)
1078 {
1079 // fetch coarse mesh index vectors
1080 const IndexTupleTypeSV& s_v = index_set_s_v[i];
1081 const IndexTupleTypeSE& s_e = index_set_s_e[i];
1082 const IndexTupleTypeST& s_t = index_set_s_t[i];
1083
1084 // create sub-index mapping objects
1085 SubIndexMappingType sim(s_v, s_t, index_set_t_v);
1086 EdgeIndexMappingType edgesim(s_v, s_e, index_set_e_v);
1087
1088 // fetch fine mesh edges-at-simplex index vectors
1089 IndexTupleType& s_0 = index_set_out[offset + 12*i + 0];
1090 IndexTupleType& s_1 = index_set_out[offset + 12*i + 1];
1091 IndexTupleType& s_2 = index_set_out[offset + 12*i + 2];
1092 IndexTupleType& s_3 = index_set_out[offset + 12*i + 3];
1093 IndexTupleType& s_4 = index_set_out[offset + 12*i + 4];
1094 IndexTupleType& s_5 = index_set_out[offset + 12*i + 5];
1095 IndexTupleType& s_6 = index_set_out[offset + 12*i + 6];
1096 IndexTupleType& s_7 = index_set_out[offset + 12*i + 7];
1097 IndexTupleType& s_8 = index_set_out[offset + 12*i + 8];
1098 IndexTupleType& s_9 = index_set_out[offset + 12*i + 9];
1099 IndexTupleType& s_10 = index_set_out[offset + 12*i + 10];
1100 IndexTupleType& s_11 = index_set_out[offset + 12*i + 11];
1101
1102 // calculate fine edges-at-simplex indices
1103
1104 // v_0-simplices
1105 s_0[0] = iot + 3*s_t[2] + sim.map(2, 0);
1106 s_0[1] = iot + 3*s_t[3] + sim.map(3, 0);
1107 s_0[2] = ioe + 2*s_e[0] + edgesim.map(0, 0);
1108 s_0[3] = iot + 3*s_t[1] + sim.map(1, 0);
1109 s_0[4] = ioe + 2*s_e[2] + edgesim.map(2, 0);
1110 s_0[5] = ioe + 2*s_e[1] + edgesim.map(1, 0);
1111
1112 s_1[0] = iot + 3*s_t[3] + sim.map(3, 0);
1113 s_1[1] = iot + 3*s_t[2] + sim.map(2, 0);
1114 s_1[2] = ios + 6*i + 0;
1115 s_1[3] = iot + 3*s_t[1] + sim.map(1, 0);
1116 s_1[4] = ios + 6*i + 1;
1117 s_1[5] = ios + 6*i + 2;
1118
1119 // v_1-simplices
1120 s_2[0] = iot + 3*s_t[3] + sim.map(3, 1);
1121 s_2[1] = iot + 3*s_t[2] + sim.map(2, 1);
1122 s_2[2] = ioe + 2*s_e[0] + edgesim.map(0, 1);
1123 s_2[3] = iot + 3*s_t[0] + sim.map(0, 0);
1124 s_2[4] = ioe + 2*s_e[3] + edgesim.map(3, 0);
1125 s_2[5] = ioe + 2*s_e[4] + edgesim.map(4, 0);
1126
1127 s_3[0] = iot + 3*s_t[2] + sim.map(2, 1);
1128 s_3[1] = iot + 3*s_t[3] + sim.map(3, 1);
1129 s_3[2] = ios + 6*i + 0;
1130 s_3[3] = iot + 3*s_t[0] + sim.map(0, 0);
1131 s_3[5] = ios + 6*i + 3;
1132 s_3[4] = ios + 6*i + 4;
1133
1134 // v_2-simplices
1135 s_4[0] = iot + 3*s_t[1] + sim.map(1, 1);
1136 s_4[1] = iot + 3*s_t[3] + sim.map(3, 2);
1137 s_4[2] = ioe + 2*s_e[1] + edgesim.map(1, 1);
1138 s_4[3] = iot + 3*s_t[0] + sim.map(0, 1);
1139 s_4[4] = ioe + 2*s_e[5] + edgesim.map(5, 0);
1140 s_4[5] = ioe + 2*s_e[3] + edgesim.map(3, 1);
1141
1142 s_5[0] = iot + 3*s_t[3] + sim.map(3, 2);
1143 s_5[1] = iot + 3*s_t[1] + sim.map(1, 1);
1144 s_5[2] = ios + 6*i + 1;
1145 s_5[3] = iot + 3*s_t[0] + sim.map(0, 1);
1146 s_5[4] = ios + 6*i + 3;
1147 s_5[5] = ios + 6*i + 5;
1148
1149 // v_3-simplices
1150 s_6[0] = iot + 3*s_t[2] + sim.map(2, 2);
1151 s_6[1] = iot + 3*s_t[1] + sim.map(1, 2);
1152 s_6[2] = ioe + 2*s_e[2] + edgesim.map(2, 1);
1153 s_6[3] = iot + 3*s_t[0] + sim.map(0, 2);
1154 s_6[4] = ioe + 2*s_e[4] + edgesim.map(4, 1);
1155 s_6[5] = ioe + 2*s_e[5] + edgesim.map(5, 1);
1156
1157 s_7[0] = iot + 3*s_t[1] + sim.map(1, 2);
1158 s_7[1] = iot + 3*s_t[2] + sim.map(2, 2);
1159 s_7[2] = ios + 6*i + 2;
1160 s_7[3] = iot + 3*s_t[0] + sim.map(0, 2);
1161 s_7[4] = ios + 6*i + 5;
1162 s_7[5] = ios + 6*i + 4;
1163
1164 //triangle simplices
1165
1166 s_8[0] = iot + 3*s_t[0] + sim.map(0, 0);
1167 s_8[1] = iot + 3*s_t[0] + sim.map(0, 1);
1168 s_8[2] = ios + 6*i + 3;
1169 s_8[3] = iot + 3*s_t[0] + sim.map(0, 2);
1170 s_8[4] = ios + 6*i + 4;
1171 s_8[5] = ios + 6*i + 5;
1172
1173 s_9[0] = iot + 3*s_t[1] + sim.map(1, 1);
1174 s_9[1] = iot + 3*s_t[1] + sim.map(1, 0);
1175 s_9[2] = ios + 6*i + 1;
1176 s_9[3] = iot + 3*s_t[1] + sim.map(1, 2);
1177 s_9[4] = ios + 6*i + 5;
1178 s_9[5] = ios + 6*i + 2;
1179
1180 s_10[0] = iot + 3*s_t[2] + sim.map(2, 0);
1181 s_10[1] = iot + 3*s_t[2] + sim.map(2, 1);
1182 s_10[2] = ios + 6*i + 0;
1183 s_10[3] = iot + 3*s_t[2] + sim.map(2, 2);
1184 s_10[4] = ios + 6*i + 2;
1185 s_10[5] = ios + 6*i + 4;
1186
1187 s_11[0] = iot + 3*s_t[3] + sim.map(3, 1);
1188 s_11[1] = iot + 3*s_t[3] + sim.map(3, 0);
1189 s_11[2] = ios + 6*i + 0;
1190 s_11[3] = iot + 3*s_t[3] + sim.map(3, 2);
1191 s_11[4] = ios + 6*i + 3;
1192 s_11[5] = ios + 6*i + 1;
1193
1194 }
1195 // return fine simplex count
1196 return 12*num_simps;
1197 }
1198 }; // StandardIndexRefiner<Simplex<3>,3,1>
1199
1205 template<>
1206 struct StandardIndexRefiner<Shape::Simplex<3>, 3, 2>
1207 {
1208 static constexpr int shape_dim = 3;
1209 static constexpr int cell_dim = 3; // second template parameter
1210 static constexpr int face_dim = 2; // third template parameter
1211
1212 typedef Shape::Simplex<shape_dim> ShapeType;
1213 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1214 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
1215 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
1216
1217 static Index refine(
1218 IndexSetType& index_set_out,
1219 const Index offset,
1220 const Index index_offsets[],
1221 const IndexSetHolderType& index_set_holder_in)
1222 {
1223 // fetch index offsets
1224 const Index iot = index_offsets[2];
1225 const Index ios = index_offsets[3];
1226
1227 // typedef index set vector reference for output index set
1228 typedef IndexSetType::IndexTupleType IndexTupleType;
1229
1230 // typedef index set types
1231 typedef IndexSet<3> IndexSetTypeTV;
1232 typedef IndexSet<4> IndexSetTypeSV;
1233 typedef IndexSet<4> IndexSetTypeST;
1234
1235 // typedef index vector type
1236 typedef IndexSetTypeSV::IndexTupleType IndexTupleTypeSV;
1237 typedef IndexSetTypeST::IndexTupleType IndexTupleTypeST;
1238
1239 // fetch the index sets
1240 const IndexSetTypeTV& index_set_t_v = index_set_holder_in.get_index_set<2,0>();
1241 const IndexSetTypeSV& index_set_s_v = index_set_holder_in.get_index_set<3,0>();
1242 const IndexSetTypeST& index_set_s_t = index_set_holder_in.get_index_set<3,2>();
1243
1244 // fetch number of coarse mesh simplices
1245 const Index num_simps = index_set_s_t.get_num_entities();
1246
1247 // typedef the sub-index mapping
1248 typedef Intern::SubIndexMapping<ShapeType, face_dim, 0> SubIndexMappingType;
1249
1250 // loop over all coarse mesh simplices
1251 FEAT_PRAGMA_OMP(parallel for)
1252 for(Index i = 0; i < num_simps; ++i)
1253 {
1254 // fetch coarse mesh vertices-at-simplex and triangle-at-simplex index vectors
1255 const IndexTupleTypeSV& s_v = index_set_s_v[i];
1256 const IndexTupleTypeST& s_t = index_set_s_t[i];
1257
1258 // create a sub-index mapping object
1259 SubIndexMappingType sim(s_v, s_t, index_set_t_v);
1260
1261 // fetch fine mesh triangle-at-simplex index vectors
1262 IndexTupleType& s_0 = index_set_out[offset + 12*i + 0];
1263 IndexTupleType& s_1 = index_set_out[offset + 12*i + 1];
1264 IndexTupleType& s_2 = index_set_out[offset + 12*i + 2];
1265 IndexTupleType& s_3 = index_set_out[offset + 12*i + 3];
1266 IndexTupleType& s_4 = index_set_out[offset + 12*i + 4];
1267 IndexTupleType& s_5 = index_set_out[offset + 12*i + 5];
1268 IndexTupleType& s_6 = index_set_out[offset + 12*i + 6];
1269 IndexTupleType& s_7 = index_set_out[offset + 12*i + 7];
1270 IndexTupleType& s_8 = index_set_out[offset + 12*i + 8];
1271 IndexTupleType& s_9 = index_set_out[offset + 12*i + 9];
1272 IndexTupleType& s_10 = index_set_out[offset + 12*i + 10];
1273 IndexTupleType& s_11 = index_set_out[offset + 12*i + 11];
1274
1275 // calculate fine triangle-at-simplex indices
1276
1277 // v_0-cube
1278 s_0[0] = iot + 4*s_t[1] + sim.map(1, 0);
1279 s_0[2] = iot + 4*s_t[2] + sim.map(2, 0);
1280 s_0[1] = iot + 4*s_t[3] + sim.map(3, 0);
1281 s_0[3] = ios + 16*i + 0;
1282
1283 s_1[0] = ios + 16*i + 8;
1284 s_1[1] = ios + 16*i + 5;
1285 s_1[2] = ios + 16*i + 4;
1286 s_1[3] = ios + 16*i + 0;
1287
1288 // v_1-cube
1289 s_2[0] = iot + 4*s_t[0] + sim.map(0, 0);
1290 s_2[1] = iot + 4*s_t[2] + sim.map(2, 1);
1291 s_2[2] = iot + 4*s_t[3] + sim.map(3, 1);
1292 s_2[3] = ios + 16*i + 1;
1293
1294 s_3[0] = ios + 16*i + 13;
1295 s_3[2] = ios + 16*i + 7;
1296 s_3[1] = ios + 16*i + 6;
1297 s_3[3] = ios + 16*i + 1;
1298
1299 // v_2-cube
1300 s_4[0] = iot + 4*s_t[0] + sim.map(0, 1);
1301 s_4[2] = iot + 4*s_t[1] + sim.map(1, 1);
1302 s_4[1] = iot + 4*s_t[3] + sim.map(3, 2);
1303 s_4[3] = ios + 16*i + 2;
1304
1305 s_5[0] = ios + 16*i + 14;
1306 s_5[1] = ios + 16*i + 10;
1307 s_5[2] = ios + 16*i + 9;
1308 s_5[3] = ios + 16*i + 2;
1309
1310 // v_3-cube
1311 s_6[0] = iot + 4*s_t[0] + sim.map(0, 2);
1312 s_6[1] = iot + 4*s_t[1] + sim.map(1, 2);
1313 s_6[2] = iot + 4*s_t[2] + sim.map(2, 2);
1314 s_6[3] = ios + 16*i + 3;
1315
1316 s_7[0] = ios + 16*i + 15;
1317 s_7[2] = ios + 16*i + 12;
1318 s_7[1] = ios + 16*i + 11;
1319 s_7[3] = ios + 16*i + 3;
1320
1321 //triangle simplices
1322
1323 s_8[0] = ios + 16*i + 15;
1324 s_8[1] = ios + 16*i + 14;
1325 s_8[2] = ios + 16*i + 13;
1326 s_8[3] = iot + 4*s_t[0] + 3;
1327
1328 s_9[0] = ios + 16*i + 12;
1329 s_9[1] = ios + 16*i + 8;
1330 s_9[2] = ios + 16*i + 10;
1331 s_9[3] = iot + 4*s_t[1] + 3;
1332
1333 s_10[0] = ios + 16*i + 11;
1334 s_10[1] = ios + 16*i + 7;
1335 s_10[2] = ios + 16*i + 5;
1336 s_10[3] = iot + 4*s_t[2] + 3;
1337
1338 s_11[0] = ios + 16*i + 9;
1339 s_11[1] = ios + 16*i + 4;
1340 s_11[2] = ios + 16*i + 6;
1341 s_11[3] = iot + 4*s_t[3] + 3;
1342
1343 }
1344 // return fine simplex count
1345 return 12*num_simps;
1346 }
1347 }; // StandardIndexRefiner<Simplex<3>,3,2>
1353 template<>
1354 struct StandardIndexRefiner<Shape::Hypercube<1>, 1, 0>
1355 {
1356 static constexpr int shape_dim = 1;
1357 static constexpr int cell_dim = 1; // second template parameter
1358 static constexpr int face_dim = 0; // third template parameter
1359
1360 typedef Shape::Hypercube<shape_dim> ShapeType;
1361 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1362 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
1363 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
1364
1365 static Index refine(
1366 IndexSetType& index_set_out,
1367 const Index offset,
1368 const Index index_offsets[],
1369 const IndexSetHolderType& index_set_holder_in)
1370 {
1371 // fetch vertex index offsets
1372 const Index iov = index_offsets[0];
1373 const Index ioe = index_offsets[1];
1374
1375 // typedef index set vector reference for output index set
1376 typedef IndexSetType::IndexTupleType IndexTupleType;
1377
1378 // typedef index set type; an edge has 2 vertices, so we need an IndexSet<2>
1379 typedef IndexSet<2> IndexSetTypeEV;
1380
1381 // typedef index vector reference type
1382 typedef IndexSetTypeEV::IndexTupleType IndexTupleTypeEV;
1383
1384 // fetch the vertices-at-edge index set
1385 const IndexSetTypeEV& index_set_e_v = index_set_holder_in.get_index_set<1,0>();
1386
1387 // fetch number of coarse mesh edges
1388 const Index num_edges = index_set_e_v.get_num_entities();
1389
1390 // loop over all coarse mesh edges
1391 FEAT_PRAGMA_OMP(parallel for)
1392 for(Index i = 0; i < num_edges; ++i)
1393 {
1394 // fetch coarse mesh vertices-at-edge index vector
1395 const IndexTupleTypeEV& e_v = index_set_e_v[i];
1396
1397 // fetch fine mesh vertices-at-edge index vectors
1398 IndexTupleType& idx0 = index_set_out[offset + 2*i + 0];
1399 IndexTupleType& idx1 = index_set_out[offset + 2*i + 1];
1400
1401 // calculate fine edge-vertex indices
1402 idx0[0] = iov + e_v[0];
1403 idx0[1] = ioe + i;
1404 idx1[0] = ioe + i;
1405 idx1[1] = iov + e_v[1];
1406 }
1407
1408 // return fine edge count
1409 return 2*num_edges;
1410 }
1411 }; // StandardIndexRefiner<Hypercube<1>,1,0>
1412
1413 /* ************************************************************************************* */
1414 /* ************************************************************************************* */
1415 /* ************************************************************************************* */
1416
1422 template<>
1423 struct StandardIndexRefiner<Shape::Hypercube<2>, 1, 0>
1424 {
1425 static constexpr int shape_dim = 2;
1426 static constexpr int cell_dim = 1; // second template parameter
1427 static constexpr int face_dim = 0; // third template parameter
1428
1429 typedef Shape::Hypercube<shape_dim> ShapeType;
1430 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1431 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
1432 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
1433
1434 static Index refine(
1435 IndexSetType& index_set_out,
1436 const Index offset,
1437 const Index index_offsets[],
1438 const IndexSetHolderType& index_set_holder_in)
1439 {
1440 // fetch vertex index offsets
1441 const Index ioe = index_offsets[1];
1442 const Index ioq = index_offsets[2];
1443
1444 // typedef index set vector reference for output index set
1445 typedef IndexSetType::IndexTupleType IndexTupleType;
1446
1447 // typedef index set type; a quad has 4 edges, so we need an IndexSet<4>
1448 typedef IndexSet<4> IndexSetTypeQE;
1449
1450 // typedef index vector type
1451 typedef IndexSetTypeQE::IndexTupleType IndexTupleTypeQE;
1452
1453 // fetch the edges-at-quad index set
1454 const IndexSetTypeQE& index_set_q_e = index_set_holder_in.get_index_set<2,1>();
1455
1456 // fetch number of coarse mesh quads
1457 const Index num_quads = index_set_q_e.get_num_entities();
1458
1459 // Each coarse mesh quad generates four fine mesh edges (e_i) upon refinement, with
1460 // each one connecting one of the fine mesh vertices generated from edge refinement
1461 // (v_i) with the vertex generated from quad refinement (x).
1462 //
1463 // +------------v_1------------+
1464 // | ^ |
1465 // | | |
1466 // | e_1 |
1467 // | | |
1468 // | | |
1469 // v_2----e_2---->x-----e_3--->v_3
1470 // | ^ |
1471 // | | |
1472 // | e_0 |
1473 // | | |
1474 // | | |
1475 // +------------v_0------------+
1476
1477
1478 // loop over all coarse mesh edges
1479 FEAT_PRAGMA_OMP(parallel for)
1480 for(Index i = 0; i < num_quads; ++i)
1481 {
1482 // fetch coarse mesh edges-at-quad index vector
1483 const IndexTupleTypeQE& q_e = index_set_q_e[i];
1484
1485 // fetch fine mesh vertices-at-edge index vectors
1486 IndexTupleType& e_0 = index_set_out[offset + 4*i + 0];
1487 IndexTupleType& e_1 = index_set_out[offset + 4*i + 1];
1488 IndexTupleType& e_2 = index_set_out[offset + 4*i + 2];
1489 IndexTupleType& e_3 = index_set_out[offset + 4*i + 3];
1490
1491 e_0[0] = ioe + q_e[0]; // v_0
1492 e_0[1] = ioq + i; // x
1493
1494 e_1[0] = ioq + i; // x
1495 e_1[1] = ioe + q_e[1]; // v_1
1496
1497 e_2[0] = ioe + q_e[2]; // v_2
1498 e_2[1] = ioq + i; // x
1499
1500 e_3[0] = ioq + i; // x
1501 e_3[1] = ioe + q_e[3]; // v_3
1502 }
1503
1504 // return fine edge count
1505 return 4*num_quads;
1506 }
1507 }; // StandardIndexRefiner<Hypercube<2>,1,0>
1508
1514 template<>
1515 struct StandardIndexRefiner<Shape::Hypercube<2>, 2, 0>
1516 {
1517 static constexpr int shape_dim = 2;
1518 static constexpr int cell_dim = 2; // second template parameter
1519 static constexpr int face_dim = 0; // third template parameter
1520
1521 typedef Shape::Hypercube<shape_dim> ShapeType;
1522 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1523 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
1524 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
1525
1526 static Index refine(
1527 IndexSetType& index_set_out,
1528 const Index offset,
1529 const Index index_offsets[],
1530 const IndexSetHolderType& index_set_holder_in)
1531 {
1532 // fetch vertex index offsets
1533 const Index iov = index_offsets[0];
1534 const Index ioe = index_offsets[1];
1535 const Index ioq = index_offsets[2];
1536
1537 // typedef index set vector reference for output index set
1538 typedef IndexSetType::IndexTupleType IndexTupleType;
1539
1540 // typedef index set types; a quad has 4 vertices and 4 edges, so both need an IndexSet<4>
1541 typedef IndexSet<4> IndexSetTypeQV;
1542 typedef IndexSet<4> IndexSetTypeQE;
1543
1544 // typedef index vector type
1545 typedef IndexSetTypeQV::IndexTupleType IndexTupleTypeQV;
1546 typedef IndexSetTypeQE::IndexTupleType IndexTupleTypeQE;
1547
1548 // fetch the quad-vertex- and the quad-edge index set
1549 const IndexSetTypeQV& index_set_q_v = index_set_holder_in.get_index_set<2,0>();
1550 const IndexSetTypeQE& index_set_q_e = index_set_holder_in.get_index_set<2,1>();
1551
1552 // fetch number of coarse mesh quads
1553 const Index num_quads = index_set_q_v.get_num_entities();
1554
1555 // Each coarse mesh quad generates four fine mesh quads (q_i) upon refinement, with
1556 // each one connecting one of the coarse mesh vertices (v_i), its two neighbor fine
1557 // mesh vertices generated from edge refinement (e_i) and the vertex generated from
1558 // quad refinement (x):
1559 //
1560 // v_2-----------e_1-----------v_3
1561 // | ^ |
1562 // | | |
1563 // | q_2 | q_3 |
1564 // | | |
1565 // | | |
1566 // e_2----------->x----------->e_3
1567 // | ^ |
1568 // | | |
1569 // | q_0 | q_1 |
1570 // | | |
1571 // | | |
1572 // v_0-----------e_0-----------v_1
1573
1574 // loop over all coarse mesh quads
1575 FEAT_PRAGMA_OMP(parallel for)
1576 for(Index i = 0; i < num_quads; ++i)
1577 {
1578 // fetch coarse mesh vertices-at-quad and edges-at-quad index vectors
1579 const IndexTupleTypeQV& q_v = index_set_q_v[i];
1580 const IndexTupleTypeQE& q_e = index_set_q_e[i];
1581
1582 // fetch fine mesh vertices-at-quad index vectors
1583 IndexTupleType& q_0 = index_set_out[offset + 4*i + 0];
1584 IndexTupleType& q_1 = index_set_out[offset + 4*i + 1];
1585 IndexTupleType& q_2 = index_set_out[offset + 4*i + 2];
1586 IndexTupleType& q_3 = index_set_out[offset + 4*i + 3];
1587
1588 // calculate fine quad-vertex indices
1589 q_0[0] = iov + q_v[0]; // v_0
1590 q_0[1] = ioe + q_e[0]; // e_0
1591 q_0[2] = ioe + q_e[2]; // e_2
1592 q_0[3] = ioq + i; // x
1593
1594 q_1[0] = ioe + q_e[0]; // e_0
1595 q_1[1] = iov + q_v[1]; // v_1
1596 q_1[2] = ioq + i; // x
1597 q_1[3] = ioe + q_e[3]; // e_3
1598
1599 q_2[0] = ioe + q_e[2]; // e_2
1600 q_2[1] = ioq + i; // x
1601 q_2[2] = iov + q_v[2]; // v_2
1602 q_2[3] = ioe + q_e[1]; // e_1
1603
1604 q_3[0] = ioq + i; // x
1605 q_3[1] = ioe + q_e[3]; // e_3
1606 q_3[2] = ioe + q_e[1]; // e_1
1607 q_3[3] = iov + q_v[3]; // v_3
1608 }
1609
1610 // return fine quad count
1611 return 4*num_quads;
1612 }
1613 }; // StandardIndexRefiner<Hypercube<2>,2,0>
1614
1620 template<>
1621 struct StandardIndexRefiner<Shape::Hypercube<2>, 2, 1>
1622 {
1623 static constexpr int shape_dim = 2;
1624 static constexpr int cell_dim = 2; // second template parameter
1625 static constexpr int face_dim = 1; // third template parameter
1626
1627 typedef Shape::Hypercube<shape_dim> ShapeType;
1628 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1629 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
1630 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
1631
1632 static Index refine(
1633 IndexSetType& index_set_out,
1634 const Index offset,
1635 const Index index_offsets[],
1636 const IndexSetHolderType& index_set_holder_in)
1637 {
1638 // fetch edge index offsets
1639 const Index ioe = index_offsets[1];
1640 const Index ioq = index_offsets[2];
1641
1642 // typedef index set vector reference for output index set
1643 typedef IndexSetType::IndexTupleType IndexTupleType;
1644
1645 // typedef index set types;
1646 typedef IndexSet<2> IndexSetTypeEV; // an edge has 2 vertices
1647 typedef IndexSet<4> IndexSetTypeQV; // a quad has 4 vertices
1648 typedef IndexSet<4> IndexSetTypeQE; // a quad has 4 edges
1649
1650 // typedef index vector type
1651 typedef IndexSetTypeQV::IndexTupleType IndexTupleTypeQV;
1652 typedef IndexSetTypeQE::IndexTupleType IndexTupleTypeQE;
1653
1654 // fetch the index sets
1655 const IndexSetTypeEV& index_set_e_v = index_set_holder_in.get_index_set<1,0>();
1656 const IndexSetTypeQV& index_set_q_v = index_set_holder_in.get_index_set<2,0>();
1657 const IndexSetTypeQE& index_set_q_e = index_set_holder_in.get_index_set<2,1>();
1658
1659 // fetch number of coarse mesh quads
1660 const Index num_quads = index_set_q_v.get_num_entities();
1661
1662 // typedef the sub-index mapping
1663 typedef Intern::SubIndexMapping<ShapeType, face_dim, 0> SubIndexMappingType;
1664
1666 //
1667 // +----e_1_0----+----e_1_1----+
1668 // | | |
1669 // | | |
1670 // e_2_1 q_2 f_1 q_3 e_3_1
1671 // | | |
1672 // | | |
1673 // +-----f_2-----+-----f_3-----+
1674 // | | |
1675 // | | |
1676 // e_2_0 q_0 f_0 q_1 e_3_0
1677 // | | |
1678 // | | |
1679 // +----e_0_0----+----e_0_1----+
1680
1681 // loop over all coarse mesh quads
1682 FEAT_PRAGMA_OMP(parallel for)
1683 for(Index i = 0; i < num_quads; ++i)
1684 {
1685 // fetch coarse mesh vertices-at-quad and edges-at-quad index vectors
1686 const IndexTupleTypeQV& q_v = index_set_q_v[i];
1687 const IndexTupleTypeQE& q_e = index_set_q_e[i];
1688
1689 // create a sub-index mapping object
1690 SubIndexMappingType sim(q_v, q_e, index_set_e_v);
1691
1692 // fetch fine mesh edges-at-quad index vectors
1693 IndexTupleType& q_0 = index_set_out[offset + 4*i + 0];
1694 IndexTupleType& q_1 = index_set_out[offset + 4*i + 1];
1695 IndexTupleType& q_2 = index_set_out[offset + 4*i + 2];
1696 IndexTupleType& q_3 = index_set_out[offset + 4*i + 3];
1697
1698 // calculate fine edges-at-quad indices
1699 q_0[0] = ioe + 2*q_e[0] + sim.map(0, 0); // e_0_0
1700 q_0[1] = ioq + 4*i + 2; // f_2
1701 q_0[2] = ioe + 2*q_e[2] + sim.map(2, 0); // e_2_0
1702 q_0[3] = ioq + 4*i + 0; // f_0
1703
1704 q_1[0] = ioe + 2*q_e[0] + sim.map(0, 1); // e_0_1
1705 q_1[1] = ioq + 4*i + 3; // f_3
1706 q_1[2] = ioq + 4*i + 0; // f_0
1707 q_1[3] = ioe + 2*q_e[3] + sim.map(3, 0); // e_3_0
1708
1709 q_2[0] = ioq + 4*i + 2; // f_2
1710 q_2[1] = ioe + 2*q_e[1] + sim.map(1, 0); // e_1_0
1711 q_2[2] = ioe + 2*q_e[2] + sim.map(2, 1); // e_2_1
1712 q_2[3] = ioq + 4*i + 1; // f_1
1713
1714 q_3[0] = ioq + 4*i + 3; // f_3
1715 q_3[1] = ioe + 2*q_e[1] + sim.map(1, 1); // e_1_1
1716 q_3[2] = ioq + 4*i + 1; // f_1
1717 q_3[3] = ioe + 2*q_e[3] + sim.map(3, 1); // e_3_1
1718 }
1719
1720 // return fine quad count
1721 return 4*num_quads;
1722 }
1723 }; // StandardIndexRefiner<Hypercube<2>,2,1>
1724
1725 /* ************************************************************************************* */
1726 /* ************************************************************************************* */
1727 /* ************************************************************************************* */
1728
1734 template<>
1735 struct StandardIndexRefiner<Shape::Hypercube<3>, 1, 0>
1736 {
1737 static constexpr int shape_dim = 3;
1738 static constexpr int cell_dim = 1; // second template parameter
1739 static constexpr int face_dim = 0; // third template parameter
1740
1741 typedef Shape::Hypercube<shape_dim> ShapeType;
1742 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1743 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
1744 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
1745
1746 static Index refine(
1747 IndexSetType& index_set_out,
1748 const Index offset,
1749 const Index index_offsets[],
1750 const IndexSetHolderType& index_set_holder_in)
1751 {
1752 // fetch vertex index offsets
1753 const Index ioq = index_offsets[2];
1754 const Index ioc = index_offsets[3];
1755
1756 // typedef index set vector reference for output index set
1757 typedef IndexSetType::IndexTupleType IndexTupleType;
1758
1759 // typedef index set type
1760 typedef IndexSet<6> IndexSetTypeCQ;
1761
1762 // typedef index vector type
1763 typedef IndexSetTypeCQ::IndexTupleType IndexTupleTypeCQ;
1764
1765 // fetch the quad-at-cube index set
1766 const IndexSetTypeCQ& index_set_c_q = index_set_holder_in.get_index_set<3,2>();
1767
1768 // fetch number of coarse mesh cubes
1769 const Index num_cubes = index_set_c_q.get_num_entities();
1770
1771 //
1772 // q_1
1773 // ^
1774 // | q_3
1775 // e_1 /
1776 // | e_3
1777 // |/
1778 // q_4-----e_4---->x-----e_5---->q_5
1779 // /^
1780 // e_2 |
1781 // / e_0
1782 // q_2 |
1783 // |
1784 // q_0
1785
1786
1787 // loop over all coarse mesh cubes
1788 FEAT_PRAGMA_OMP(parallel for)
1789 for(Index i = 0; i < num_cubes; ++i)
1790 {
1791 // fetch coarse mesh quad-at-cube index vector
1792 const IndexTupleTypeCQ& c_q = index_set_c_q[i];
1793
1794 // fetch fine mesh vertices-at-edge index vectors
1795 IndexTupleType& e_0 = index_set_out[offset + 6*i + 0];
1796 IndexTupleType& e_1 = index_set_out[offset + 6*i + 1];
1797 IndexTupleType& e_2 = index_set_out[offset + 6*i + 2];
1798 IndexTupleType& e_3 = index_set_out[offset + 6*i + 3];
1799 IndexTupleType& e_4 = index_set_out[offset + 6*i + 4];
1800 IndexTupleType& e_5 = index_set_out[offset + 6*i + 5];
1801
1802 e_0[0] = ioq + c_q[0]; // q_0
1803 e_0[1] = ioc + i; // x
1804
1805 e_1[0] = ioc + i; // x
1806 e_1[1] = ioq + c_q[1]; // q_1
1807
1808 e_2[0] = ioq + c_q[2]; // q_2
1809 e_2[1] = ioc + i; // x
1810
1811 e_3[0] = ioc + i; // x
1812 e_3[1] = ioq + c_q[3]; // q_3
1813
1814 e_4[0] = ioq + c_q[4]; // q_4
1815 e_4[1] = ioc + i; // x
1816
1817 e_5[0] = ioc + i; // x
1818 e_5[1] = ioq + c_q[5]; // q_5
1819 }
1820
1821 // return fine edge count
1822 return 6*num_cubes;
1823 }
1824 }; // StandardIndexRefiner<Hypercube<3>,1,0>
1825
1831 template<>
1832 struct StandardIndexRefiner<Shape::Hypercube<3>, 2, 0>
1833 {
1834 static constexpr int shape_dim = 3;
1835 static constexpr int cell_dim = 2; // second template parameter
1836 static constexpr int face_dim = 0; // third template parameter
1837
1838 typedef Shape::Hypercube<shape_dim> ShapeType;
1839 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
1840 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
1841 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
1842
1843 static Index refine(
1844 IndexSetType& index_set_out,
1845 const Index offset,
1846 const Index index_offsets[],
1847 const IndexSetHolderType& index_set_holder_in)
1848 {
1849 // fetch vertex index offsets
1850 const Index ioe = index_offsets[1];
1851 const Index ioq = index_offsets[2];
1852 const Index ioc = index_offsets[3];
1853
1854 // typedef index set vector reference for output index set
1855 typedef IndexSetType::IndexTupleType IndexTupleType;
1856
1857 // typedef index set types
1858 typedef IndexSet<6> IndexSetTypeCQ;
1859 typedef IndexSet<12> IndexSetTypeCE;
1860
1861 // typedef index vector type
1862 typedef IndexSetTypeCQ::IndexTupleType IndexTupleTypeCQ;
1863 typedef IndexSetTypeCE::IndexTupleType IndexTupleTypeCE;
1864
1865 // fetch the index sets
1866 const IndexSetTypeCQ& index_set_c_q = index_set_holder_in.get_index_set<3,2>();
1867 const IndexSetTypeCE& index_set_c_e = index_set_holder_in.get_index_set<3,1>();
1868
1869 // fetch number of coarse mesh cubes
1870 const Index num_cubes = index_set_c_q.get_num_entities();
1871
1872 // x-normal surfaces
1873 // e_3
1874 // /|
1875 // / |
1876 // / |
1877 // q_1 s |
1878 // /| 3 |
1879 // / | q_3
1880 // / | /|
1881 // e_2 s | / |
1882 // | 2 |/ |
1883 // | x s |
1884 // | /| 1 |
1885 // | / | e_1
1886 // |/ | /
1887 // q_2 s | /
1888 // | 0 |/
1889 // | q_0
1890 // | / z
1891 // | / ^ y
1892 // |/ |/
1893 // e_0 ---->x
1894
1895
1896 // y-normal surfaces
1897 //
1898 // e_6-----------q_1-----------e_7
1899 // | | |
1900 // | | |
1901 // | s_6 | s_7 |
1902 // | | |
1903 // | | |
1904 // q_4------------x------------q_5
1905 // | | |
1906 // | | |
1907 // | s_4 | s_5 | z
1908 // | | | ^ y
1909 // | | | |/
1910 // e_4-----------q_0-----------e_5 ---->x
1911
1912
1913 // z-normal surfaces
1914 //
1915 // e_10----------q_3----------e_11
1916 // / / /
1917 // / s_10 / s_11 /
1918 // / / /
1919 // q_4----------- x ---------- q_5
1920 // / / / z
1921 // / s_8 / s_9 / ^ y
1922 // / / / |/
1923 // e_8----------q_2----------e_9 ---->x
1924 //
1925 //
1926
1927 // loop over all coarse mesh cubes
1928 FEAT_PRAGMA_OMP(parallel for)
1929 for(Index i = 0; i < num_cubes; ++i)
1930 {
1931 // fetch coarse mesh quad-at-cube and edges-at-cube index vectors
1932 const IndexTupleTypeCQ& c_q = index_set_c_q[i];
1933 const IndexTupleTypeCE& c_e = index_set_c_e[i];
1934
1935 // fetch fine mesh vertices-at-quad index vectors
1936 IndexTupleType& s_0 = index_set_out[offset + 12*i + 0];
1937 IndexTupleType& s_1 = index_set_out[offset + 12*i + 1];
1938 IndexTupleType& s_2 = index_set_out[offset + 12*i + 2];
1939 IndexTupleType& s_3 = index_set_out[offset + 12*i + 3];
1940 IndexTupleType& s_4 = index_set_out[offset + 12*i + 4];
1941 IndexTupleType& s_5 = index_set_out[offset + 12*i + 5];
1942 IndexTupleType& s_6 = index_set_out[offset + 12*i + 6];
1943 IndexTupleType& s_7 = index_set_out[offset + 12*i + 7];
1944 IndexTupleType& s_8 = index_set_out[offset + 12*i + 8];
1945 IndexTupleType& s_9 = index_set_out[offset + 12*i + 9];
1946 IndexTupleType& s_10 = index_set_out[offset + 12*i + 10];
1947 IndexTupleType& s_11 = index_set_out[offset + 12*i + 11];
1948
1949 // calculate fine quad-vertex indices
1950
1951 // x-normal
1952 s_0[0] = ioe + c_e[0]; // e_0
1953 s_0[1] = ioq + c_q[0]; // q_0
1954 s_0[2] = ioq + c_q[2]; // q_2
1955 s_0[3] = ioc + i; // x
1956
1957 s_1[0] = ioq + c_q[0]; // q_0
1958 s_1[1] = ioe + c_e[1]; // e_1
1959 s_1[2] = ioc + i; // x
1960 s_1[3] = ioq + c_q[3]; // q_3
1961
1962 s_2[0] = ioq + c_q[2]; // q_2
1963 s_2[1] = ioc + i; // x
1964 s_2[2] = ioe + c_e[2]; // e_2
1965 s_2[3] = ioq + c_q[1]; // q_1
1966
1967 s_3[0] = ioc + i; // x
1968 s_3[1] = ioq + c_q[3]; // q_3
1969 s_3[2] = ioq + c_q[1]; // q_1
1970 s_3[3] = ioe + c_e[3]; // e_3
1971
1972 // y-normal
1973 s_4[0] = ioe + c_e[4]; // e_4
1974 s_4[1] = ioq + c_q[0]; // q_0
1975 s_4[2] = ioq + c_q[4]; // q_4
1976 s_4[3] = ioc + i; // x
1977
1978 s_5[0] = ioq + c_q[0]; // q_0
1979 s_5[1] = ioe + c_e[5]; // e_5
1980 s_5[2] = ioc + i; // x
1981 s_5[3] = ioq + c_q[5]; // q_5
1982
1983 s_6[0] = ioq + c_q[4]; // q_4
1984 s_6[1] = ioc + i; // x
1985 s_6[2] = ioe + c_e[6]; // e_6
1986 s_6[3] = ioq + c_q[1]; // q_1
1987
1988 s_7[0] = ioc + i; // x
1989 s_7[1] = ioq + c_q[5]; // q_5
1990 s_7[2] = ioq + c_q[1]; // q_1
1991 s_7[3] = ioe + c_e[7]; // e_7
1992
1993 // z-normal
1994 s_8[0] = ioe + c_e[8]; // e_8
1995 s_8[1] = ioq + c_q[2]; // q_2
1996 s_8[2] = ioq + c_q[4]; // q_4
1997 s_8[3] = ioc + i; // x
1998
1999 s_9[0] = ioq + c_q[2]; // q_2
2000 s_9[1] = ioe + c_e[9]; // e_9
2001 s_9[2] = ioc + i; // x
2002 s_9[3] = ioq + c_q[5]; // q_5
2003
2004 s_10[0] = ioq + c_q[4]; // q_4
2005 s_10[1] = ioc + i; // x
2006 s_10[2] = ioe + c_e[10]; // e_10
2007 s_10[3] = ioq + c_q[3]; // q_3
2008
2009 s_11[0] = ioc + i; // x
2010 s_11[1] = ioq + c_q[5]; // q_5
2011 s_11[2] = ioq + c_q[3]; // q_3
2012 s_11[3] = ioe + c_e[11]; // e_11
2013 }
2014 // return fine quad count
2015 return 12*num_cubes;
2016 }
2017 }; // StandardIndexRefiner<Hypercube<3>,2,0>
2018
2024 template<>
2025 struct StandardIndexRefiner<Shape::Hypercube<3>, 3, 0>
2026 {
2027 static constexpr int shape_dim = 3;
2028 static constexpr int cell_dim = 3; // second template parameter
2029 static constexpr int face_dim = 0; // third template parameter
2030
2031 typedef Shape::Hypercube<shape_dim> ShapeType;
2032 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
2033 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
2034 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
2035
2036 static Index refine(
2037 IndexSetType& index_set_out,
2038 const Index offset,
2039 const Index index_offsets[],
2040 const IndexSetHolderType& index_set_holder_in)
2041 {
2042 // fetch vertex index offsets
2043 const Index iov = index_offsets[0];
2044 const Index ioe = index_offsets[1];
2045 const Index ioq = index_offsets[2];
2046 const Index ioc = index_offsets[3];
2047
2048 // typedef index set vector reference for output index set
2049 typedef IndexSetType::IndexTupleType IndexTupleType;
2050
2051 // typedef index set types
2052 typedef IndexSet<6> IndexSetTypeCQ;
2053 typedef IndexSet<12> IndexSetTypeCE;
2054 typedef IndexSet<8> IndexSetTypeCV;
2055
2056 // typedef index vector type
2057 typedef IndexSetTypeCQ::IndexTupleType IndexTupleTypeCQ;
2058 typedef IndexSetTypeCE::IndexTupleType IndexTupleTypeCE;
2059 typedef IndexSetTypeCV::IndexTupleType IndexTupleTypeCV;
2060
2061 // fetch the index sets
2062 const IndexSetTypeCQ& index_set_c_q = index_set_holder_in.get_index_set<3,2>();
2063 const IndexSetTypeCE& index_set_c_e = index_set_holder_in.get_index_set<3,1>();
2064 const IndexSetTypeCV& index_set_c_v = index_set_holder_in.get_index_set<3,0>();
2065
2066 // fetch number of coarse mesh cubes
2067 const Index num_cubes = index_set_c_q.get_num_entities();
2068
2069 //
2070 // v_6-----------e_3-----------v_7
2071 // /| /| /|
2072 // / | / | / |
2073 // / | / | / |
2074 // e_6------------q_1-----------e_7 |
2075 // /| | /| | /| |
2076 // / | e_10-----/-|--q_3------/-|-e_11
2077 // / | /| / | /| / | /|
2078 // v_4-----------e_2-----------v_5 | / |
2079 // | |/ | | |/ | | |/ |
2080 // | q_4--------|---x---------|--q_5 |
2081 // | /| | | /| | | /| |
2082 // | / | v_2----|-/-|--e_1----|-/-|--v_3
2083 // |/ | / |/ | / |/ | /
2084 // e_8-----------q_2-----------e_9 | /
2085 // | |/ | |/ | |/
2086 // | e_4--------|--q_0--------|--e_5
2087 // | / | / | /
2088 // | / | / | /
2089 // |/ |/ |/
2090 // v_0-----------e_0-----------v_1
2091 //
2092
2093 // loop over all coarse mesh cubes
2094 FEAT_PRAGMA_OMP(parallel for)
2095 for(Index i = 0; i < num_cubes; ++i)
2096 {
2097 // fetch coarse mesh index vectors
2098 const IndexTupleTypeCQ& c_q = index_set_c_q[i];
2099 const IndexTupleTypeCE& c_e = index_set_c_e[i];
2100 const IndexTupleTypeCV& c_v = index_set_c_v[i];
2101
2102 // fetch fine mesh vertices-at-cube index vectors
2103 IndexTupleType& c_0 = index_set_out[offset + 8*i + 0];
2104 IndexTupleType& c_1 = index_set_out[offset + 8*i + 1];
2105 IndexTupleType& c_2 = index_set_out[offset + 8*i + 2];
2106 IndexTupleType& c_3 = index_set_out[offset + 8*i + 3];
2107 IndexTupleType& c_4 = index_set_out[offset + 8*i + 4];
2108 IndexTupleType& c_5 = index_set_out[offset + 8*i + 5];
2109 IndexTupleType& c_6 = index_set_out[offset + 8*i + 6];
2110 IndexTupleType& c_7 = index_set_out[offset + 8*i + 7];
2111
2112 // calculate fine cube-vertex indices
2113
2114 c_0[0] = iov + c_v[0]; // v_0
2115 c_0[1] = ioe + c_e[0]; // e_0
2116 c_0[2] = ioe + c_e[4]; // e_4
2117 c_0[3] = ioq + c_q[0]; // q_0
2118 c_0[4] = ioe + c_e[8]; // e_8
2119 c_0[5] = ioq + c_q[2]; // q_2
2120 c_0[6] = ioq + c_q[4]; // q_4
2121 c_0[7] = ioc + i; // x
2122
2123 c_1[0] = ioe + c_e[0]; // e_0
2124 c_1[1] = iov + c_v[1]; // v_1
2125 c_1[2] = ioq + c_q[0]; // q_0
2126 c_1[3] = ioe + c_e[5]; // e_5
2127 c_1[4] = ioq + c_q[2]; // q_2
2128 c_1[5] = ioe + c_e[9]; // e_9
2129 c_1[6] = ioc + i; // x
2130 c_1[7] = ioq + c_q[5]; // q_5
2131
2132 c_2[0] = ioe + c_e[4]; // e_4
2133 c_2[1] = ioq + c_q[0]; // q_0
2134 c_2[2] = iov + c_v[2]; // v_2
2135 c_2[3] = ioe + c_e[1]; // e_1
2136 c_2[4] = ioq + c_q[4]; // q_4
2137 c_2[5] = ioc + i; // x
2138 c_2[6] = ioe + c_e[10];// e_10
2139 c_2[7] = ioq + c_q[3]; // q_3
2140
2141 c_3[0] = ioq + c_q[0]; // q_0
2142 c_3[1] = ioe + c_e[5]; // e_5
2143 c_3[2] = ioe + c_e[1]; // e_1
2144 c_3[3] = iov + c_v[3]; // v_3
2145 c_3[4] = ioc + i; // x
2146 c_3[5] = ioq + c_q[5]; // q_5
2147 c_3[6] = ioq + c_q[3]; // q_3
2148 c_3[7] = ioe + c_e[11];// e_11
2149
2150 c_4[0] = ioe + c_e[8]; // e_8
2151 c_4[1] = ioq + c_q[2]; // q_2
2152 c_4[2] = ioq + c_q[4]; // q_4
2153 c_4[3] = ioc + i; // x
2154 c_4[4] = iov + c_v[4]; // v_4
2155 c_4[5] = ioe + c_e[2]; // e_2
2156 c_4[6] = ioe + c_e[6]; // e_6
2157 c_4[7] = ioq + c_q[1]; // q_1
2158
2159 c_5[0] = ioq + c_q[2]; // q_2
2160 c_5[1] = ioe + c_e[9]; // e_9
2161 c_5[2] = ioc + i; // x
2162 c_5[3] = ioq + c_q[5]; // q_5
2163 c_5[4] = ioe + c_e[2]; // e_2
2164 c_5[5] = iov + c_v[5]; // v_5
2165 c_5[6] = ioq + c_q[1]; // q_1
2166 c_5[7] = ioe + c_e[7]; // e_7
2167
2168 c_6[0] = ioq + c_q[4]; // q_4
2169 c_6[1] = ioc + i; // x
2170 c_6[2] = ioe + c_e[10];// e_10
2171 c_6[3] = ioq + c_q[3]; // q_3
2172 c_6[4] = ioe + c_e[6]; // e_6
2173 c_6[5] = ioq + c_q[1]; // q_1
2174 c_6[6] = iov + c_v[6]; // v_6
2175 c_6[7] = ioe + c_e[3]; // e_3
2176
2177 c_7[0] = ioc + i; // x
2178 c_7[1] = ioq + c_q[5]; // q_5
2179 c_7[2] = ioq + c_q[3]; // q_3
2180 c_7[3] = ioe + c_e[11];// e_11
2181 c_7[4] = ioq + c_q[1]; // q_1
2182 c_7[5] = ioe + c_e[7]; // e_7
2183 c_7[6] = ioe + c_e[3]; // e_3
2184 c_7[7] = iov + c_v[7]; // v_7
2185 }
2186 // return fine cube count
2187 return 8*num_cubes;
2188 }
2189 }; // StandardIndexRefiner<Hypercube<3>,3,0>
2190
2196 template<>
2197 struct StandardIndexRefiner<Shape::Hypercube<3>, 2, 1>
2198 {
2199 static constexpr int shape_dim = 3;
2200 static constexpr int cell_dim = 2; // second template parameter
2201 static constexpr int face_dim = 1; // third template parameter
2202
2203 typedef Shape::Hypercube<shape_dim> ShapeType;
2204 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
2205 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
2206 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
2207
2208 static Index refine(
2209 IndexSetType& index_set_out,
2210 const Index offset,
2211 const Index index_offsets[],
2212 const IndexSetHolderType& index_set_holder_in)
2213 {
2214 // fetch index offsets
2215 const Index ioq = index_offsets[2];
2216 const Index ioc = index_offsets[3];
2217
2218 // typedef index set vector reference for output index set
2219 typedef IndexSetType::IndexTupleType IndexTupleType;
2220
2221 // typedef index set types;
2222 typedef IndexSet<4> IndexSetTypeQV;
2223 typedef IndexSet<8> IndexSetTypeCV;
2224 typedef IndexSet<6> IndexSetTypeCQ;
2225
2226 // typedef index vector type
2227 typedef IndexSetTypeCV::IndexTupleType IndexTupleTypeCV;
2228 typedef IndexSetTypeCQ::IndexTupleType IndexTupleTypeCQ;
2229
2230 // fetch the index sets
2231 const IndexSetTypeQV& index_set_q_v = index_set_holder_in.get_index_set<2,0>();
2232 const IndexSetTypeCV& index_set_c_v = index_set_holder_in.get_index_set<3,0>();
2233 const IndexSetTypeCQ& index_set_c_q = index_set_holder_in.get_index_set<3,2>();
2234
2235 // fetch number of coarse mesh cubes
2236 const Index num_cubes = index_set_c_q.get_num_entities();
2237
2238 // typedef the sub-index mapping
2239 typedef Intern::SubIndexMapping<ShapeType, cell_dim, face_dim> SubIndexMappingType;
2240
2241 // x-normal surfaces
2242 // e_3
2243 // /|
2244 // / |
2245 // / |
2246 // q_1 s |
2247 // /| 3 |
2248 // / | q_3
2249 // / | /|
2250 // e_2 s | / |
2251 // | 2 |/ |
2252 // | x s |
2253 // | /| 1 |
2254 // | / | e_1
2255 // |/ | /
2256 // q_2 s | /
2257 // | 0 |/
2258 // | q_0
2259 // | / z
2260 // | / ^ y
2261 // |/ |/
2262 // e_0 ---->x
2263
2264
2265 // y-normal surfaces
2266 //
2267 // e_6-----------q_1-----------e_7
2268 // | | |
2269 // | | |
2270 // | s_6 | s_7 |
2271 // | | |
2272 // | | |
2273 // q_4------------x------------q_5
2274 // | | |
2275 // | | |
2276 // | s_4 | s_5 | z
2277 // | | | ^ y
2278 // | | | |/
2279 // e_4-----------q_0-----------e_5 ---->x
2280
2281
2282 // z-normal surfaces
2283 //
2284 // e_10----------q_3----------e_11
2285 // / / /
2286 // / s_10 / s_11 /
2287 // / / /
2288 // q_4----------- x ---------- q_5
2289 // / / / z
2290 // / s_8 / s_9 / ^ y
2291 // / / / |/
2292 // e_8----------q_2----------e_9 ---->x
2293 //
2294 //
2295
2296 // loop over all coarse mesh cubes
2297 FEAT_PRAGMA_OMP(parallel for)
2298 for(Index i = 0; i < num_cubes; ++i)
2299 {
2300 // fetch coarse mesh vertices-at-quad and edges-at-quad index vectors
2301 const IndexTupleTypeCV& c_v = index_set_c_v[i];
2302 const IndexTupleTypeCQ& c_q = index_set_c_q[i];
2303
2304 // create a sub-index mapping object
2305 SubIndexMappingType sim(c_v, c_q, index_set_q_v);
2306
2307 // fetch fine mesh edges-at-quad index vectors
2308 IndexTupleType& q_0 = index_set_out[offset + 12*i + 0];
2309 IndexTupleType& q_1 = index_set_out[offset + 12*i + 1];
2310 IndexTupleType& q_2 = index_set_out[offset + 12*i + 2];
2311 IndexTupleType& q_3 = index_set_out[offset + 12*i + 3];
2312 IndexTupleType& q_4 = index_set_out[offset + 12*i + 4];
2313 IndexTupleType& q_5 = index_set_out[offset + 12*i + 5];
2314 IndexTupleType& q_6 = index_set_out[offset + 12*i + 6];
2315 IndexTupleType& q_7 = index_set_out[offset + 12*i + 7];
2316 IndexTupleType& q_8 = index_set_out[offset + 12*i + 8];
2317 IndexTupleType& q_9 = index_set_out[offset + 12*i + 9];
2318 IndexTupleType& q_10 = index_set_out[offset + 12*i + 10];
2319 IndexTupleType& q_11 = index_set_out[offset + 12*i + 11];
2320
2321 // calculate fine edges-at-quad indices
2322
2323 // x-normal
2324 q_0[0] = ioq + 4*c_q[0] + sim.map(0, 0); // e_0-q_0
2325 q_0[1] = ioc + 6*i + 2; // q_2-x
2326 q_0[2] = ioq + 4*c_q[2] + sim.map(2, 0); // e_0-q_2
2327 q_0[3] = ioc + 6*i + 0; // q_0-x
2328
2329 q_1[0] = ioq + 4*c_q[0] + sim.map(0, 1); // q_0-e_1
2330 q_1[1] = ioc + 6*i + 3; // x-q_3
2331 q_1[2] = ioc + 6*i + 0; // q_0-x
2332 q_1[3] = ioq + 4*c_q[3] + sim.map(3, 0); // e_1-q_3
2333
2334 q_2[0] = ioc + 6*i + 2; // q_2-x
2335 q_2[1] = ioq + 4*c_q[1] + sim.map(1, 0); // e_2-q_1
2336 q_2[2] = ioq + 4*c_q[2] + sim.map(2, 1); // q_2-e_2
2337 q_2[3] = ioc + 6*i + 1; // x-q_1
2338
2339 q_3[0] = ioc + 6*i + 3; // x-q_3
2340 q_3[1] = ioq + 4*c_q[1] + sim.map(1, 1); // q_1-e_3
2341 q_3[2] = ioc + 6*i + 1; // x-q_1
2342 q_3[3] = ioq + 4*c_q[3] + sim.map(3, 1); // q_3-e_3
2343
2344 // y-normal
2345 q_4[0] = ioq + 4*c_q[0] + sim.map(0, 2); // e_4-q_0
2346 q_4[1] = ioc + 6*i + 4; // q_4-x
2347 q_4[2] = ioq + 4*c_q[4] + sim.map(4, 0); // e_4-q_4
2348 q_4[3] = ioc + 6*i + 0; // q_0-x
2349
2350 q_5[0] = ioq + 4*c_q[0] + sim.map(0, 3); // q_0-e_5
2351 q_5[1] = ioc + 6*i + 5; // x-q_5
2352 q_5[2] = ioc + 6*i + 0; // q_0-x
2353 q_5[3] = ioq + 4*c_q[5] + sim.map(5, 0); // e_5-q_5
2354
2355 q_6[0] = ioc + 6*i + 4; // q_4-x
2356 q_6[1] = ioq + 4*c_q[1] + sim.map(1, 2); // e_6-q_1
2357 q_6[2] = ioq + 4*c_q[4] + sim.map(4, 1); // q_4-e_6
2358 q_6[3] = ioc + 6*i + 1; // x-q_1
2359
2360 q_7[0] = ioc + 6*i + 5; // x-q_5
2361 q_7[1] = ioq + 4*c_q[1] + sim.map(1, 3); // q_1-e_7
2362 q_7[2] = ioc + 6*i + 1; // x-q_1
2363 q_7[3] = ioq + 4*c_q[5] + sim.map(5, 1); // q_5-e_7
2364
2365 // z-normal
2366 q_8[0] = ioq + 4*c_q[2] + sim.map(2, 2); // e_8-q_2
2367 q_8[1] = ioc + 6*i + 4; // q_4-x
2368 q_8[2] = ioq + 4*c_q[4] + sim.map(4, 2); // e_8-q_4
2369 q_8[3] = ioc + 6*i + 2; // q_2-x
2370
2371 q_9[0] = ioq + 4*c_q[2] + sim.map(2, 3); // q_2-e_9
2372 q_9[1] = ioc + 6*i + 5; // x-q_5
2373 q_9[2] = ioc + 6*i + 2; // q_2-x
2374 q_9[3] = ioq + 4*c_q[5] + sim.map(5, 2); // e_9-q_5
2375
2376 q_10[0] = ioc + 6*i + 4; // q_4-x
2377 q_10[1] = ioq + 4*c_q[3] + sim.map(3, 2); // e_10-q_3
2378 q_10[2] = ioq + 4*c_q[4] + sim.map(4, 3); // q_4-e_10
2379 q_10[3] = ioc + 6*i + 3; // x-q_3
2380
2381 q_11[0] = ioc + 6*i + 5; // x-q_5
2382 q_11[1] = ioq + 4*c_q[3] + sim.map(3, 3); // q_3-e_11
2383 q_11[2] = ioc + 6*i + 3; // x-q_3
2384 q_11[3] = ioq + 4*c_q[5] + sim.map(5, 3); // q_5-e_11
2385 }
2386 // return fine quad count
2387 return 12*num_cubes;
2388 }
2389 }; // StandardIndexRefiner<Hypercube<3>,2,1>
2390
2396 template<>
2397 struct StandardIndexRefiner<Shape::Hypercube<3>, 3, 1>
2398 {
2399 static constexpr int shape_dim = 3;
2400 static constexpr int cell_dim = 3; // second template parameter
2401 static constexpr int face_dim = 1; // third template parameter
2402
2403 typedef Shape::Hypercube<shape_dim> ShapeType;
2404 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
2405 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
2406 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
2407
2408 static Index refine(
2409 IndexSetType& index_set_out,
2410 const Index offset,
2411 const Index index_offsets[],
2412 const IndexSetHolderType& index_set_holder_in)
2413 {
2414 // fetch edge index offsets
2415 const Index ioe = index_offsets[1];
2416 const Index ioq = index_offsets[2];
2417 const Index ioc = index_offsets[3];
2418
2419 // typedef index set vector reference for output index set
2420 typedef IndexSetType::IndexTupleType IndexTupleType;
2421
2422 // typedef index set types
2423 typedef IndexSet<2> IndexSetTypeEV;
2424 typedef IndexSet<4> IndexSetTypeQV;
2425 typedef IndexSet<8> IndexSetTypeCV;
2426 typedef IndexSet<12> IndexSetTypeCE;
2427 typedef IndexSet<6> IndexSetTypeCQ;
2428
2429 // typedef index vector type
2430 typedef IndexSetTypeCV::IndexTupleType IndexTupleTypeCV;
2431 typedef IndexSetTypeCE::IndexTupleType IndexTupleTypeCE;
2432 typedef IndexSetTypeCQ::IndexTupleType IndexTupleTypeCQ;
2433
2434 // fetch the index sets
2435 const IndexSetTypeEV& index_set_e_v = index_set_holder_in.get_index_set<1,0>();
2436 const IndexSetTypeQV& index_set_q_v = index_set_holder_in.get_index_set<2,0>();
2437 const IndexSetTypeCV& index_set_c_v = index_set_holder_in.get_index_set<3,0>();
2438 const IndexSetTypeCE& index_set_c_e = index_set_holder_in.get_index_set<3,1>();
2439 const IndexSetTypeCQ& index_set_c_q = index_set_holder_in.get_index_set<3,2>();
2440
2441 // fetch number of coarse mesh cubes
2442 const Index num_cubes = index_set_c_q.get_num_entities();
2443
2444 // typedef the sub-index mapping
2445 typedef Intern::SubIndexMapping<ShapeType, 2, 1> SubIndexMappingType;
2446 typedef Intern::SubIndexMapping<ShapeType, face_dim, 0> EdgeIndexMappingType;
2447
2448 //
2449 // v_6-----------e_3-----------v_7
2450 // /| /| /|
2451 // / | / | / |
2452 // / | / | / |
2453 // e_6------------q_1-----------e_7 |
2454 // /| | /| | /| |
2455 // / | e_10-----/-|--q_3------/-|-e_11
2456 // / | /| / | /| / | /|
2457 // v_4-----------e_2-----------v_5 | / |
2458 // | |/ | | |/ | | |/ |
2459 // | q_4--------|---x---------|--q_5 |
2460 // | /| | | /| | | /| |
2461 // | / | v_2----|-/-|--e_1----|-/-|--v_3
2462 // |/ | / |/ | / |/ | /
2463 // e_8-----------q_2-----------e_9 | /
2464 // | |/ | |/ | |/
2465 // | e_4--------|--q_0--------|--e_5
2466 // | / | / | /
2467 // | / | / | /
2468 // |/ |/ |/
2469 // v_0-----------e_0-----------v_1
2470 //
2471
2472 // loop over all coarse mesh cubes
2473 FEAT_PRAGMA_OMP(parallel for)
2474 for(Index i = 0; i < num_cubes; ++i)
2475 {
2476 // fetch coarse mesh index vectors
2477 const IndexTupleTypeCV& c_v = index_set_c_v[i];
2478 const IndexTupleTypeCE& c_e = index_set_c_e[i];
2479 const IndexTupleTypeCQ& c_q = index_set_c_q[i];
2480
2481 // create sub-index mapping objects
2482 SubIndexMappingType sim(c_v, c_q, index_set_q_v);
2483 EdgeIndexMappingType edgesim(c_v, c_e, index_set_e_v);
2484
2485 // fetch fine mesh edges-at-cube index vectors
2486 IndexTupleType& c_0 = index_set_out[offset + 8*i + 0];
2487 IndexTupleType& c_1 = index_set_out[offset + 8*i + 1];
2488 IndexTupleType& c_2 = index_set_out[offset + 8*i + 2];
2489 IndexTupleType& c_3 = index_set_out[offset + 8*i + 3];
2490 IndexTupleType& c_4 = index_set_out[offset + 8*i + 4];
2491 IndexTupleType& c_5 = index_set_out[offset + 8*i + 5];
2492 IndexTupleType& c_6 = index_set_out[offset + 8*i + 6];
2493 IndexTupleType& c_7 = index_set_out[offset + 8*i + 7];
2494
2495 // calculate fine edges-at-cube indices
2496
2497 // v_0-cube
2498 c_0[0] = ioe + 2*c_e[0] + edgesim.map(0, 0); // v_0-e_0
2499 c_0[1] = ioq + 4*c_q[0] + sim.map(0, 2); // e_4-q_0
2500 c_0[2] = ioq + 4*c_q[2] + sim.map(2, 2); // e_8-q_2
2501 c_0[3] = ioc + 6*i + 4; // q_4-x
2502 c_0[4] = ioe + 2*c_e[4] + edgesim.map(4, 0); // v_0-e_4
2503 c_0[5] = ioq + 4*c_q[0] + sim.map(0, 0); // e_0-q_0
2504 c_0[6] = ioq + 4*c_q[4] + sim.map(4, 2); // e_8-q_4
2505 c_0[7] = ioc + 6*i + 2; // q_2-x
2506 c_0[8] = ioe + 2*c_e[8] + edgesim.map(8, 0); // v_0-e_8
2507 c_0[9] = ioq + 4*c_q[2] + sim.map(2, 0); // e_0-q_2
2508 c_0[10] = ioq + 4*c_q[4] + sim.map(4, 0); // e_4-q_4
2509 c_0[11] = ioc + 6*i + 0; // q_0-x
2510
2511 // v_1-cube
2512 c_1[0] = ioe + 2*c_e[0] + edgesim.map(0,1); // e_0-v_1
2513 c_1[1] = ioq + 4*c_q[0] + sim.map(0, 3); // q_0-e_5
2514 c_1[2] = ioq + 4*c_q[2] + sim.map(2,3); // q_2-e_9
2515 c_1[3] = ioc + 6*i + 5; // x-q_5
2516 c_1[4] = ioq + 4*c_q[0] + sim.map(0, 0); // e_0-q_0
2517 c_1[5] = ioe + 2*c_e[5] + edgesim.map(5, 0); // v_1-e_5
2518 c_1[6] = ioc + 6*i + 2; // q_2-x
2519 c_1[7] = ioq + 4*c_q[5] + sim.map(5, 2); // e_9-q_5
2520 c_1[8] = ioq + 4*c_q[2] + sim.map(2,0); // e_0-q_2
2521 c_1[9] = ioe + 2*c_e[9] + edgesim.map(9, 0); // v_1-e_9
2522 c_1[10] = ioc + 6*i + 0; // q_0-x
2523 c_1[11] = ioq + 4*c_q[5] + sim.map(5, 0); // e_5-q_5
2524
2525 // v_2-cube
2526 c_2[0] = ioq + 4*c_q[0] + sim.map(0, 2); // e_4-q_0
2527 c_2[1] = ioe + 2*c_e[1] + edgesim.map(1, 0); // v_2-e_1
2528 c_2[2] = ioc + 6*i + 4; // q_4-x
2529 c_2[3] = ioq + 4*c_q[3] + sim.map(3, 2); // e_10-q_3
2530 c_2[4] = ioe + 2*c_e[4] + edgesim.map(4, 1); // e_4-v_2
2531 c_2[5] = ioq + 4*c_q[0] + sim.map(0, 1); // q_0-e_1
2532 c_2[6] = ioq + 4*c_q[4] + sim.map(4, 3); // q_4-e_10
2533 c_2[7] = ioc + 6*i + 3; // x-q_3
2534 c_2[8] = ioq + 4*c_q[4] + sim.map(4, 0); // e_4-q_4
2535 c_2[9] = ioc + 6*i + 0; // q_0-x
2536 c_2[10] = ioe + 2*c_e[10] + edgesim.map(10, 0); // v_2-e_10
2537 c_2[11] = ioq + 4*c_q[3] + sim.map(3, 0); // e_1-q_3
2538
2539 // v_3-cube
2540 c_3[0] = ioq + 4*c_q[0] + sim.map(0, 3); // q_0-e_5
2541 c_3[1] = ioe + 2*c_e[1] + edgesim.map(1, 1); // e_1-v_3
2542 c_3[2] = ioc + 6*i + 5; // x-q_5
2543 c_3[3] = ioq + 4*c_q[3] + sim.map(3, 3); // q_3-e_11
2544 c_3[4] = ioq + 4*c_q[0] + sim.map(0, 1); // q_0-e_1
2545 c_3[5] = ioe + 2*c_e[5] + edgesim.map(5, 1); // e_5-v_3
2546 c_3[6] = ioc + 6*i + 3; // x-q_3
2547 c_3[7] = ioq + 4*c_q[5] + sim.map(5, 3); // q_5-e_11
2548 c_3[8] = ioc + 6*i + 0; // q_0-x
2549 c_3[9] = ioq + 4*c_q[5] + sim.map(5, 0); // e_5-q_5
2550 c_3[10] = ioq + 4*c_q[3] + sim.map(3, 0); // e_1-q_3
2551 c_3[11] = ioe + 2*c_e[11] + edgesim.map(11, 0); // v_3-e_11
2552
2553 // v_4-cube
2554 c_4[0] = ioq + 4*c_q[2] + sim.map(2, 2); // e_8-q_2
2555 c_4[1] = ioc + 6*i + 4; // q_4-x
2556 c_4[2] = ioe + 2*c_e[2] + edgesim.map(2, 0); // v_4-e_2
2557 c_4[3] = ioq + 4*c_q[1] + sim.map(1, 2); // e_6-q_1
2558 c_4[4] = ioq + 4*c_q[4] + sim.map(4, 2); // e_8-q_4
2559 c_4[5] = ioc + 6*i + 2; // q_2-x
2560 c_4[6] = ioe + 2*c_e[6] + edgesim.map(6, 0); // v_4-e_6
2561 c_4[7] = ioq + 4*c_q[1] + sim.map(1, 0); // e_2-q_1
2562 c_4[8] = ioe + 2*c_e[8] + edgesim.map(8, 1); // e_8-v_4
2563 c_4[9] = ioq + 4*c_q[2] + sim.map(2, 1); // q_2-e_2
2564 c_4[10] = ioq + 4*c_q[4] + sim.map(4, 1); // q_4-e_6
2565 c_4[11] = ioc + 6*i + 1; // x-q_1
2566
2567 // v_5-cube
2568 c_5[0] = ioq + 4*c_q[2] + sim.map(2, 3); // q_2-e_9
2569 c_5[1] = ioc + 6*i + 5; // x-q_5
2570 c_5[2] = ioe + 2*c_e[2] + edgesim.map(2, 1); // e_2-v-5
2571 c_5[3] = ioq + 4*c_q[1] + sim.map(1, 3); // q_1-e_7
2572 c_5[4] = ioc + 6*i + 2; // q_2-x
2573 c_5[5] = ioq + 4*c_q[5] + sim.map(5, 2); // e_9-q_5
2574 c_5[6] = ioq + 4*c_q[1] + sim.map(1, 0); // e_2-q_1
2575 c_5[7] = ioe + 2*c_e[7] + edgesim.map(7, 0); // v_5-e_7
2576 c_5[8] = ioq + 4*c_q[2] + sim.map(2, 1); // q_2-e_2
2577 c_5[9] = ioe + 2*c_e[9] + edgesim.map(9, 1); // e_9-v_5
2578 c_5[10] = ioc + 6*i + 1; // x-q_1
2579 c_5[11] = ioq + 4*c_q[5] + sim.map(5, 1); // q_5-e_7
2580
2581 // v_6-cube
2582 c_6[0] = ioc + 6*i + 4; // q_4-x
2583 c_6[1] = ioq + 4*c_q[3] + sim.map(3, 2); // e_10-q_3
2584 c_6[2] = ioq + 4*c_q[1] + sim.map(1, 2); // e_6-q_1
2585 c_6[3] = ioe + 2*c_e[3] + edgesim.map(3, 0); // v_6-e_3
2586 c_6[4] = ioq + 4*c_q[4] + sim.map(4, 3); // q_4-e_10
2587 c_6[5] = ioc + 6*i + 3; // x-q_3
2588 c_6[6] = ioe + 2*c_e[6] + edgesim.map(6, 1); // e_6-v_6
2589 c_6[7] = ioq + 4*c_q[1] + sim.map(1, 1); // q_1-e_3
2590 c_6[8] = ioq + 4*c_q[4] + sim.map(4, 1); // q_4-e_6
2591 c_6[9] = ioc + 6*i + 1; // x-q_1
2592 c_6[10] = ioe + 2*c_e[10] + edgesim.map(10, 1); // e_10-v_6
2593 c_6[11] = ioq + 4*c_q[3] + sim.map(3, 1); // q_3-e_3
2594
2595 // v_7-cube
2596 c_7[0] = ioc + 6*i + 5; // x-q_5
2597 c_7[1] = ioq + 4*c_q[3] + sim.map(3, 3); // q_3-e_11
2598 c_7[2] = ioq + 4*c_q[1] + sim.map(1, 3); // q_1-e_7
2599 c_7[3] = ioe + 2*c_e[3] + edgesim.map(3, 1); // e_3-v_7
2600 c_7[4] = ioc + 6*i + 3; // x-q_3
2601 c_7[5] = ioq + 4*c_q[5] + sim.map(5, 3); // q_5-e_11
2602 c_7[6] = ioq + 4*c_q[1] + sim.map(1, 1); // q_1-e_3
2603 c_7[7] = ioe + 2*c_e[7] + edgesim.map(7, 1); // e_7-v_7
2604 c_7[8] = ioc + 6*i + 1; // x-q_1
2605 c_7[9] = ioq + 4*c_q[5] + sim.map(5, 1); // q_5-e_7
2606 c_7[10] = ioq + 4*c_q[3] + sim.map(3, 1); // q_3-e_3
2607 c_7[11] = ioe + 2*c_e[11] + edgesim.map(11, 1); // e_11-v_7
2608 }
2609 // return fine cube count
2610 return 8*num_cubes;
2611 }
2612 }; // StandardIndexRefiner<Hypercube<3>,3,1>
2613
2619 template<>
2620 struct StandardIndexRefiner<Shape::Hypercube<3>, 3, 2>
2621 {
2622 static constexpr int shape_dim = 3;
2623 static constexpr int cell_dim = 3; // second template parameter
2624 static constexpr int face_dim = 2; // third template parameter
2625
2626 typedef Shape::Hypercube<shape_dim> ShapeType;
2627 typedef Shape::FaceTraits<ShapeType, cell_dim>::ShapeType CellType;
2628 typedef IndexSet<Shape::FaceTraits<CellType, face_dim>::count> IndexSetType;
2629 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
2630
2631 static Index refine(
2632 IndexSetType& index_set_out,
2633 const Index offset,
2634 const Index index_offsets[],
2635 const IndexSetHolderType& index_set_holder_in)
2636 {
2637 // fetch quad index offsets
2638 const Index ioq = index_offsets[2];
2639 const Index ioc = index_offsets[3];
2640
2641 // typedef index set vector reference for output index set
2642 typedef IndexSetType::IndexTupleType IndexTupleType;
2643
2644 // typedef index set types
2645 typedef IndexSet<4> IndexSetTypeQV;
2646 typedef IndexSet<8> IndexSetTypeCV;
2647 typedef IndexSet<6> IndexSetTypeCQ;
2648
2649 // typedef index vector type
2650 typedef IndexSetTypeCV::IndexTupleType IndexTupleTypeCV;
2651 typedef IndexSetTypeCQ::IndexTupleType IndexTupleTypeCQ;
2652
2653 // fetch the index sets
2654 const IndexSetTypeQV& index_set_q_v = index_set_holder_in.get_index_set<2,0>();
2655 const IndexSetTypeCV& index_set_c_v = index_set_holder_in.get_index_set<3,0>();
2656 const IndexSetTypeCQ& index_set_c_q = index_set_holder_in.get_index_set<3,2>();
2657
2658 // fetch number of coarse mesh cubes
2659 const Index num_cubes = index_set_c_q.get_num_entities();
2660
2661 // typedef the sub-index mapping
2662 typedef Intern::SubIndexMapping<ShapeType, face_dim, 0> SubIndexMappingType;
2663
2664 //
2665 // v_6-----------e_3-----------v_7
2666 // /| /| /|
2667 // / | / | / |
2668 // / | / | / |
2669 // e_6------------q_1-----------e_7 |
2670 // /| | /| | /| |
2671 // / | e_10-----/-|--q_3------/-|-e_11
2672 // / | /| / | /| / | /|
2673 // v_4-----------e_2-----------v_5 | / |
2674 // | |/ | | |/ | | |/ |
2675 // | q_4--------|---x---------|--q_5 |
2676 // | /| | | /| | | /| |
2677 // | / | v_2----|-/-|--e_1----|-/-|--v_3
2678 // |/ | / |/ | / |/ | /
2679 // e_8-----------q_2-----------e_9 | /
2680 // | |/ | |/ | |/
2681 // | e_4--------|--q_0--------|--e_5
2682 // | / | / | /
2683 // | / | / | /
2684 // |/ |/ |/
2685 // v_0-----------e_0-----------v_1
2686 //
2687
2688 // loop over all coarse mesh cubes
2689 FEAT_PRAGMA_OMP(parallel for)
2690 for(Index i = 0; i < num_cubes; ++i)
2691 {
2692 // fetch coarse mesh vertices-at-cube and quads-at-cube index vectors
2693 const IndexTupleTypeCV& c_v = index_set_c_v[i];
2694 const IndexTupleTypeCQ& c_q = index_set_c_q[i];
2695
2696 // create a sub-index mapping object
2697 SubIndexMappingType sim(c_v, c_q, index_set_q_v);
2698
2699 // fetch fine mesh quads-at-cube index vectors
2700 IndexTupleType& c_0 = index_set_out[offset + 8*i + 0];
2701 IndexTupleType& c_1 = index_set_out[offset + 8*i + 1];
2702 IndexTupleType& c_2 = index_set_out[offset + 8*i + 2];
2703 IndexTupleType& c_3 = index_set_out[offset + 8*i + 3];
2704 IndexTupleType& c_4 = index_set_out[offset + 8*i + 4];
2705 IndexTupleType& c_5 = index_set_out[offset + 8*i + 5];
2706 IndexTupleType& c_6 = index_set_out[offset + 8*i + 6];
2707 IndexTupleType& c_7 = index_set_out[offset + 8*i + 7];
2708
2709 // calculate fine quad-at-cube indices
2710
2711 // v_0-cube
2712 c_0[0] = ioq + 4*c_q[0] + sim.map(0, 0); // v_0-e_0-e_4-q_0
2713 c_0[1] = ioc + 12*i + 8; // e_8-q_2-q_4-x
2714 c_0[2] = ioq + 4*c_q[2] + sim.map(2, 0); // v_0-e_0-e_8-q_2
2715 c_0[3] = ioc + 12*i + 4; // e_4-q_0-q_4-x
2716 c_0[4] = ioq + 4*c_q[4] + sim.map(4, 0); // v_0-e_4-e_8-q_4
2717 c_0[5] = ioc + 12*i + 0; // e_0-q_0-q_2-x
2718
2719
2720 // v_1-cube
2721 c_1[0] = ioq + 4*c_q[0] + sim.map(0, 1); // e_0-v_1-q_0-e_5
2722 c_1[1] = ioc + 12*i + 9; // q_2-e_9-x-q_5
2723 c_1[2] = ioq + 4*c_q[2] + sim.map(2, 1); // e_0-v_1-q_2-e_9
2724 c_1[3] = ioc + 12*i + 5; // q_0-e_5-x-q_5
2725 c_1[4] = ioc + 12*i + 0; // e_0-q_0-q_2-x
2726 c_1[5] = ioq + 4*c_q[5] + sim.map(5, 0); // v_1-e_5-e_9-q_5
2727
2728 // v_2-cube
2729 c_2[0] = ioq + 4*c_q[0] + sim.map(0, 2); // e_4-q_0-v_2-e_1
2730 c_2[1] = ioc + 12*i + 10; // q_4-x-e_10-q_3
2731 c_2[2] = ioc + 12*i + 4; // e_4-q_0-q_4-x
2732 c_2[3] = ioq + 4*c_q[3] + sim.map(3, 0); // v_2-e_1-e_10-q_3
2733 c_2[4] = ioq + 4*c_q[4] + sim.map(4, 1); // e_4-v_2-q_4-e_10
2734 c_2[5] = ioc + 12*i + 1; // q_0-e_1-x-q_3
2735
2736 // v_3-cube
2737 c_3[0] = ioq + 4*c_q[0] + sim.map(0, 3); // q_0-e_5-e_1-v_3
2738 c_3[1] = ioc + 12*i + 11; // x-q_5-q_3-e_11
2739 c_3[2] = ioc + 12*i + 5; // q_0-e_5-x-q_5
2740 c_3[3] = ioq + 4*c_q[3] + sim.map(3, 1); // e_1-v_3-q_3-e_11
2741 c_3[4] = ioc + 12*i + 1; // q_0-e_1-x-q_3
2742 c_3[5] = ioq + 4*c_q[5] + sim.map(5, 1); // e_5-v_3-q_5-e_11
2743
2744 // v_4-cube
2745 c_4[0] = ioc + 12*i + 8; // e_8-q_2-q_4-x
2746 c_4[1] = ioq + 4*c_q[1] + sim.map(1, 0); // v_4-e_2-e_6-q_1
2747 c_4[2] = ioq + 4*c_q[2] + sim.map(2, 2); // e_8-q_2-v_4-e_2
2748 c_4[3] = ioc + 12*i + 6; // q_4-x-e_6-q_1
2749 c_4[4] = ioq + 4*c_q[4] + sim.map(4, 2); // e_8-q_4-v_4-e_6
2750 c_4[5] = ioc + 12*i + 2; // q_2-x-e_2-q_1
2751
2752 // v_5-cube
2753 c_5[0] = ioc + 12*i + 9; // q_2-e_9-x-q_5
2754 c_5[1] = ioq + 4*c_q[1] + sim.map(1, 1); // e_2-v_5-q_1-e_7
2755 c_5[2] = ioq + 4*c_q[2] + sim.map(2, 3); // q_2-e_9-e_2-v_5
2756 c_5[3] = ioc + 12*i + 7; // x-q_5-q_1-e_7
2757 c_5[4] = ioc + 12*i + 2; // q_2-x-e_2-q_1
2758 c_5[5] = ioq + 4*c_q[5] + sim.map(5, 2); // e_9-q_5-v_5-e_7
2759
2760 // v_6-cube
2761 c_6[0] = ioc + 12*i + 10; // q_4-x-e_10-q_3
2762 c_6[1] = ioq + 4*c_q[1] + sim.map(1, 2); // e_6-q_1-v_6-e_3
2763 c_6[2] = ioc + 12*i + 6; // q_4-x-e_6-q_1
2764 c_6[3] = ioq + 4*c_q[3] + sim.map(3, 2); // e_10-q_3-v_6-e_3
2765 c_6[4] = ioq + 4*c_q[4] + sim.map(4, 3); // q_4-e_10-e_6-v_6
2766 c_6[5] = ioc + 12*i + 3; // x-q_3-q_1-e_3
2767
2768 // v_7-cube
2769 c_7[0] = ioc + 12*i + 11; // x-q_5-q_3-e_11
2770 c_7[1] = ioq + 4*c_q[1] + sim.map(1, 3); // q_1-e_7-e_3-v_7
2771 c_7[2] = ioc + 12*i + 7; // x-q_5-q_1-e_7
2772 c_7[3] = ioq + 4*c_q[3] + sim.map(3, 3); // q_3-e_11-e_3-v_7
2773 c_7[4] = ioc + 12*i + 3; // x-q_3-q_1-e_3
2774 c_7[5] = ioq + 4*c_q[5] + sim.map(5, 3); // q_5-e_11-e_7-v_7
2775 }
2776 // return fine cube count
2777 return 8*num_cubes;
2778 }
2779 }; // StandardIndexRefiner<Hypercube<3>,3,2>
2780
2781 /* *************************************************************************************** */
2782 /* *************************************************************************************** */
2783 /* *************************************************************************************** */
2784 /* *************************************************************************************** */
2785 /* *************************************************************************************** */
2786
2787 template<
2788 typename Shape_,
2789 int cell_dim_,
2790 int face_dim_,
2791 // The following "dummy" argument is necessary for partial specialization;
2792 // shape_dim_ *must* coincide with Shape_::dimension !!!
2793 int shape_dim_ = Shape_::dimension>
2794 struct IndexRefineShapeWrapper
2795 {
2796 static_assert(shape_dim_ == Shape_::dimension, "invalid shape dimension");
2797 // the case shape_dim_ = cell_dim_ is specialized below
2798 static_assert(shape_dim_ > cell_dim_, "invalid cell dimension");
2799 static_assert(cell_dim_ > face_dim_, "invalid face dimension");
2800 static_assert(face_dim_ >= 0, "invalid face_dimension");
2801
2802 typedef Shape_ ShapeType;
2803 typedef typename Shape::FaceTraits<ShapeType, cell_dim_>::ShapeType CellType;
2804 typedef IndexSet<Shape::FaceTraits<CellType, face_dim_>::count> IndexSetType;
2805 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
2806
2807 static String name()
2808 {
2809 return "IndexRefineShapeWrapper<" + Shape_::name() + "," + stringify(cell_dim_) + ","
2810 + stringify(face_dim_) + "," + stringify(shape_dim_) + ">";
2811 }
2812
2813 static Index refine(
2814 IndexSetType& index_set_out,
2815 const Index index_offsets[],
2816 const IndexSetHolderType& index_set_holder_in)
2817 {
2818 typedef typename Shape::FaceTraits<ShapeType, shape_dim_ - 1>::ShapeType FacetType;
2819
2820 // refine facets
2821 Index offset = IndexRefineShapeWrapper<FacetType, cell_dim_, face_dim_>
2822 ::refine(index_set_out, index_offsets, index_set_holder_in);
2823
2824 // call index refiner and return new offset
2825 Index num_faces = StandardIndexRefiner<Shape_, cell_dim_, face_dim_>
2826 ::refine(index_set_out, offset, index_offsets, index_set_holder_in);
2827
2828 // validate number of created indices
2829 Index num_shapes = index_set_holder_in.template get_index_set<shape_dim_, 0>().get_num_entities();
2830 XASSERTM(num_faces == (StandardRefinementTraits<ShapeType, cell_dim_>::count * num_shapes),
2831 "IndexRefiner output does not match StdRefTraits prediction");
2832
2833 // calculate new offset and return
2834 return offset + num_faces;
2835 }
2836 };
2837
2838 template<
2839 typename Shape_,
2840 int cell_dim_,
2841 int face_dim_>
2842 struct IndexRefineShapeWrapper<Shape_, cell_dim_, face_dim_, cell_dim_>
2843 {
2844 static_assert(cell_dim_ == Shape_::dimension, "invalid shape dimension");
2845 static_assert(cell_dim_ > face_dim_, "invalid face dimension");
2846 static_assert(face_dim_ >= 0, "invalid face_dimension");
2847
2848 typedef Shape_ ShapeType;
2849 typedef typename Shape::FaceTraits<ShapeType, cell_dim_>::ShapeType CellType;
2850 typedef IndexSet<Shape::FaceTraits<CellType, face_dim_>::count> IndexSetType;
2851 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
2852
2853 static String name()
2854 {
2855 return "IndexRefineShapeWrapper<" + Shape_::name() + "," + stringify(cell_dim_) + ","
2856 + stringify(face_dim_) + "," + stringify(cell_dim_) + ">";
2857 }
2858
2859 static Index refine(
2860 IndexSetType& index_set_out,
2861 const Index index_offsets[],
2862 const IndexSetHolderType& index_set_holder_in)
2863 {
2864 // call index refiner
2865 Index num_faces = StandardIndexRefiner<Shape_, cell_dim_, face_dim_>
2866 ::refine(index_set_out, 0, index_offsets, index_set_holder_in);
2867
2868 // validate number of created indices
2869 Index num_shapes = index_set_holder_in.template get_index_set<cell_dim_, 0>().get_num_entities();
2870 XASSERTM(num_faces == (StandardRefinementTraits<ShapeType, cell_dim_>::count * num_shapes),
2871 "IndexRefiner output does not match StdRefTraits prediction");
2872
2873 // return offset
2874 return num_faces;
2875 }
2876 };
2877
2878 /* ********************************************************************************************************* */
2879
2880 template<
2881 typename Shape_,
2882 int cell_dim_,
2883 int face_dim_ = cell_dim_ - 1>
2884 struct IndexRefineFaceWrapper
2885 {
2886 // the case face_dim_ = 0 is handled by the partial specialization below
2887 static_assert(face_dim_ > 0, "invalid face dimension");
2888 static_assert(cell_dim_ > face_dim_, "invalid cell dimension");
2889 static_assert(cell_dim_ <= Shape_::dimension, "invalid cell dimension");
2890
2891 typedef Shape_ ShapeType;
2892 typedef typename Shape::FaceTraits<ShapeType, cell_dim_>::ShapeType CellType;
2893 typedef IndexSet<Shape::FaceTraits<CellType, face_dim_>::count> IndexSetType;
2894 typedef IndexSetWrapper<CellType, face_dim_> IndexSetWrapperType;
2895 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
2896
2897 static String name()
2898 {
2899 return "IndexRefineFaceWrapper<" + Shape_::name() + "," + stringify(cell_dim_) + ","
2900 + stringify(face_dim_) + ">";
2901 }
2902
2903 static void refine(
2904 IndexSetWrapperType& index_set_wrapper_out,
2905 const Index num_entities[],
2906 const IndexSetHolderType& index_set_holder_in)
2907 {
2908 // recursive call of IndexRefineFaceWrapper
2909 IndexRefineFaceWrapper<Shape_, cell_dim_, face_dim_ - 1>
2910 ::refine(index_set_wrapper_out, num_entities, index_set_holder_in);
2911
2912 // fetch output index set
2913 IndexSetType& index_set_out = index_set_wrapper_out.template get_index_set<face_dim_>();
2914
2915 // calculate index offsets
2916 Index index_offsets[Shape_::dimension+1];
2917 EntityCounter<StandardRefinementTraits, Shape_, face_dim_>::offset(index_offsets, num_entities);
2918
2919 // call refinement shape wrapper
2920 IndexRefineShapeWrapper<Shape_, cell_dim_, face_dim_>
2921 ::refine(index_set_out, index_offsets, index_set_holder_in);
2922 }
2923 };
2924
2925 template<
2926 typename Shape_,
2927 int cell_dim_>
2928 struct IndexRefineFaceWrapper<Shape_, cell_dim_, 0>
2929 {
2930 static_assert(cell_dim_ > 0, "invalid cell dimension");
2931 static_assert(cell_dim_ <= Shape_::dimension, "invalid cell dimension");
2932
2933 typedef Shape_ ShapeType;
2934 typedef typename Shape::FaceTraits<ShapeType, cell_dim_>::ShapeType CellType;
2935 typedef IndexSet<Shape::FaceTraits<CellType, 0>::count> IndexSetType;
2936 typedef IndexSetWrapper<CellType, 0> IndexSetWrapperType;
2937 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
2938
2939 static String name()
2940 {
2941 return "IndexRefineFaceWrapper<" + Shape_::name() + "," + stringify(cell_dim_) + ",0>";
2942 }
2943
2944 static void refine(
2945 IndexSetWrapperType& index_set_wrapper_out,
2946 const Index num_entities[],
2947 const IndexSetHolderType& index_set_holder_in)
2948 {
2949 // fetch output index set
2950 IndexSetType& index_set_out = index_set_wrapper_out.template get_index_set<0>();
2951
2952 // calculate index offsets
2953 Index index_offsets[Shape_::dimension+1];
2954 EntityCounter<StandardRefinementTraits, Shape_, 0>::offset(index_offsets, num_entities);
2955
2956 // call refinement shape wrapper
2957 IndexRefineShapeWrapper<Shape_, cell_dim_, 0>
2958 ::refine(index_set_out, index_offsets, index_set_holder_in);
2959 }
2960 };
2961
2962 /* ********************************************************************************************************* */
2963
2964 template<
2965 typename Shape_,
2966 int cell_dim_ = Shape_::dimension>
2967 struct IndexRefineWrapper
2968 {
2969 // the case cell_dim_ = 1 is handled by the partial specialization below
2970 static_assert(cell_dim_ > 1, "invalid cell dimension");
2971 static_assert(cell_dim_ <= Shape_::dimension, "invalid cell dimension");
2972
2973 typedef Shape_ ShapeType;
2974 typedef typename Shape::FaceTraits<ShapeType, cell_dim_>::ShapeType CellType;
2975 typedef IndexSetWrapper<CellType> IndexSetWrapperType;
2976 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
2977
2978 static String name()
2979 {
2980 return "IndexRefineWrapper<" + Shape_::name() + "," + stringify(cell_dim_) + ">";
2981 }
2982
2983 static void refine(
2984 IndexSetHolderType& index_set_holder_out,
2985 const Index num_entities[],
2986 const IndexSetHolderType& index_set_holder_in)
2987 {
2988 // recursive call of IndexRefineWrapper
2989 IndexRefineWrapper<Shape_, cell_dim_ - 1>
2990 ::refine(index_set_holder_out, num_entities, index_set_holder_in);
2991
2992 // fetch output index set wrapper
2993 IndexSetWrapperType& index_set_wrapper_out = index_set_holder_out
2994 .template get_index_set_wrapper<cell_dim_>();
2995
2996 // call face wrapper
2997 IndexRefineFaceWrapper<Shape_, cell_dim_>
2998 ::refine(index_set_wrapper_out, num_entities, index_set_holder_in);
2999 }
3000 };
3001
3002 template<typename Shape_>
3003 struct IndexRefineWrapper<Shape_, 1>
3004 {
3005 static_assert(1 <= Shape_::dimension, "invalid shape dimension");
3006
3007 typedef Shape_ ShapeType;
3008 typedef typename Shape::FaceTraits<ShapeType, 1>::ShapeType CellType;
3009 typedef IndexSetWrapper<CellType> IndexSetWrapperType;
3010 typedef IndexSetHolder<ShapeType> IndexSetHolderType;
3011
3012 static String name()
3013 {
3014 return "IndexRefineWrapper<" + Shape_::name() + ",1>";
3015 }
3016
3017 static void refine(
3018 IndexSetHolderType& index_set_holder_out,
3019 const Index num_entities[],
3020 const IndexSetHolderType& index_set_holder_in)
3021 {
3022 // fetch output index set wrapper
3023 IndexSetWrapperType& index_set_wrapper_out = index_set_holder_out.template get_index_set_wrapper<1>();
3024
3025 // call face wrapper
3026 IndexRefineFaceWrapper<Shape_, 1>::refine(index_set_wrapper_out, num_entities, index_set_holder_in);
3027 }
3028 };
3029 } // namespace Intern
3031 } // namespace Geometry
3032} // 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