FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
unit_cube_patch_generator.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#include <kernel/geometry/conformal_mesh.hpp>
9#include <kernel/geometry/common_factories.hpp>
10#include <kernel/geometry/mesh_node.hpp>
11#include <kernel/geometry/intern/face_index_mapping.hpp>
12
13#include <vector>
14#include <stdexcept>
15
16namespace FEAT
17{
18 namespace Geometry
19 {
21 template<typename Mesh_>
22 class UnitCubePatchGenerator;
23
25 template<typename Coord_>
26 class UnitCubePatchGenerator<ConformalMesh<Shape::Hypercube<1>, 1, Coord_>>
27 {
28 public:
29 typedef Shape::Hypercube<1> ShapeType;
30 typedef Coord_ CoordType;
31 typedef ConformalMesh<ShapeType, 1, Coord_> MeshType;
32 typedef MeshPart<MeshType> PartType;
33 typedef RootMeshNode<MeshType> MeshNodeType;
34
35 private:
36 class RootMeshFactory :
37 public Geometry::UnitCubeFactory<MeshType>
38 {
39 public:
40 typedef typename MeshType::VertexSetType VertexSetType;
41
42 private:
43 const int _i, _n;
44
45 public:
46 RootMeshFactory(int i, int n) :
47 _i(i), _n(n)
48 {
49 }
50
51 virtual void fill_vertex_set(VertexSetType& vertex_set) override
52 {
53 vertex_set[0][0] = CoordType(_i ) / CoordType(_n);
54 vertex_set[1][0] = CoordType(_i+1) / CoordType(_n);
55 }
56 };
57
58 class Halo0Factory :
59 public Geometry::Factory<PartType>
60 {
61 public:
62 typedef Geometry::Factory<PartType> BaseClass;
63 typedef typename BaseClass::AttributeSetContainer AttributeSetContainer;
64 typedef typename BaseClass::IndexSetHolderType IndexSetHolderType;
65 typedef typename BaseClass::TargetSetHolderType TargetSetHolderType;
66
67 private:
68 const int _k;
69
70 public:
71 explicit Halo0Factory(int k) :
72 _k(k)
73 {
74 }
75
76 virtual Index get_num_entities(int dim) override
77 {
78 return Index(dim == 0 ? 1 : 0);
79 }
80
81 virtual void fill_attribute_sets(AttributeSetContainer&) override
82 {
83 }
84
85 virtual void fill_index_sets(std::unique_ptr<IndexSetHolderType>&) override
86 {
87 }
88
89 virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
90 {
91 target_set_holder.template get_target_set<0>()[0] = Index(_k);
92 }
93
94 };
95
96 static std::unique_ptr<PartType> create_halo0(int k)
97 {
98 Halo0Factory factory(k);
99 return factory.make_unique();
100 }
101
102 public:
103 static Index create_unique(int rank, int nprocs, std::unique_ptr<MeshNodeType>& node, std::vector<int>& ranks)
104 {
105 MeshNodeType* nnode = nullptr;
106 std::vector<Index> cranks, ctags;
107 Index level = Index(create(rank, nprocs, nnode, cranks, ctags));
108 for(auto it = cranks.begin(); it != cranks.end(); ++it)
109 ranks.push_back(int(*it));
110 node.reset(nnode);
111 return level;
112 }
113
114 // deprecated: use make_unique instead
115 static int create(int rank, int nprocs, MeshNodeType*& node, std::vector<Index>& ranks, std::vector<Index>& ctags)
116 {
117 XASSERT(nprocs > 0);
118 XASSERT((rank >= 0) && (rank < nprocs));
119
120 // determine slice count and refinement level
121 int level(0);
122 int n(1);
123 for(; n < nprocs; n *= 2, ++level) {}
124 XASSERTM(n == nprocs, "number of processes must be a power of 2");
125
126 // decompose rank to (i)
127 const int ii(rank);
128
129 // create root mesh node
130 {
131 RootMeshFactory factory(ii, n);
132 node = new MeshNodeType(factory.make_unique());
133 }
134
135 // left neighbor
136 if(ii > 0)
137 {
138 ranks.push_back(Index(ii-1));
139 ctags.push_back(Index(ii-1));
140 node->add_halo(int(ranks.back()), create_halo0(0));
141 node->add_mesh_part("bnd:0", nullptr);
142 }
143 else
144 {
145 node->add_mesh_part("bnd:0", create_halo0(0));
146 }
147 // right neighbor
148 if(ii+1 < n)
149 {
150 ranks.push_back(Index(ii+1));
151 ctags.push_back(Index(ii));
152 node->add_halo(int(ranks.back()), create_halo0(1));
153 node->add_mesh_part("bnd:1", nullptr);
154 }
155 else
156 {
157 node->add_mesh_part("bnd:1", create_halo0(1));
158 }
159
160 // return refinement level
161 return level;
162 }
163 }; // UnitCubePatchGenerator<Hypercube<1>>
164
165 template<typename Coord_>
166 class UnitCubePatchGenerator<ConformalMesh<Shape::Hypercube<2>, 2, Coord_>>
167 {
168 public:
169 typedef Shape::Hypercube<2> ShapeType;
170 typedef Coord_ CoordType;
171 typedef ConformalMesh<ShapeType, 2, Coord_> MeshType;
172 typedef MeshPart<MeshType> PartType;
173 typedef RootMeshNode<MeshType> MeshNodeType;
174
175 private:
176 class RootMeshFactory :
177 public Geometry::UnitCubeFactory<MeshType>
178 {
179 public:
180 typedef typename MeshType::VertexSetType VertexSetType;
181
182 private:
183 const int _i, _j, _n;
184
185 public:
186 RootMeshFactory(int i, int j, int n) :
187 _i(i), _j(j), _n(n)
188 {
189 }
190
191 virtual void fill_vertex_set(VertexSetType& vertex_set) override
192 {
193 vertex_set[0][0] = CoordType(_i ) / CoordType(_n);
194 vertex_set[0][1] = CoordType(_j ) / CoordType(_n);
195 vertex_set[1][0] = CoordType(_i+1) / CoordType(_n);
196 vertex_set[1][1] = CoordType(_j ) / CoordType(_n);
197 vertex_set[2][0] = CoordType(_i ) / CoordType(_n);
198 vertex_set[2][1] = CoordType(_j+1) / CoordType(_n);
199 vertex_set[3][0] = CoordType(_i+1) / CoordType(_n);
200 vertex_set[3][1] = CoordType(_j+1) / CoordType(_n);
201 }
202 };
203
204 class Halo0Factory :
205 public Geometry::Factory<PartType>
206 {
207 public:
208 typedef Geometry::Factory<PartType> BaseClass;
209 typedef typename BaseClass::AttributeSetContainer AttributeSetContainer;
210 typedef typename BaseClass::IndexSetHolderType IndexSetHolderType;
211 typedef typename BaseClass::TargetSetHolderType TargetSetHolderType;
212
213 private:
214 const int _k;
215
216 public:
217 explicit Halo0Factory(int k) :
218 _k(k)
219 {
220 }
221
222 virtual Index get_num_entities(int dim) override
223 {
224 return Index(dim == 0 ? 1 : 0);
225 }
226
227 virtual void fill_attribute_sets(AttributeSetContainer&) override
228 {
229 }
230
231 virtual void fill_index_sets(std::unique_ptr<IndexSetHolderType>&) override
232 {
233 }
234
235 virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
236 {
237 target_set_holder.template get_target_set<0>()[0] = Index(_k);
238 }
239 };
240
241 class Halo1Factory :
242 public Geometry::Factory<PartType>
243 {
244 public:
245 typedef Geometry::Factory<PartType> BaseClass;
246 typedef typename BaseClass::AttributeSetContainer AttributeSetContainer;
247 typedef typename BaseClass::IndexSetHolderType IndexSetHolderType;
248 typedef typename BaseClass::TargetSetHolderType TargetSetHolderType;
249
250 private:
251 const int _k;
252
253 public:
254 explicit Halo1Factory(int k) :
255 _k(k)
256 {
257 }
258
259 virtual Index get_num_entities(int dim) override
260 {
261 return Index(dim == 0 ? 2 : (dim == 1 ? 1 : 0));
262 }
263
264 virtual void fill_attribute_sets(AttributeSetContainer&) override
265 {
266 }
267
268 virtual void fill_index_sets(std::unique_ptr<IndexSetHolderType>&) override
269 {
270 }
271
272 virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
273 {
274 typedef Geometry::Intern::FaceIndexMapping<Shape::Quadrilateral, 1, 0> Fim;
275 target_set_holder.template get_target_set<0>()[0] = Index(Fim::map(_k, 0));
276 target_set_holder.template get_target_set<0>()[1] = Index(Fim::map(_k, 1));
277 target_set_holder.template get_target_set<1>()[0] = Index(_k);
278 }
279 };
280
281 static std::unique_ptr<PartType> create_halo0(int k)
282 {
283 Halo0Factory factory(k);
284 return factory.make_unique();
285 }
286
287 static std::unique_ptr<PartType> create_halo1(int k)
288 {
289 Halo1Factory factory(k);
290 return factory.make_unique();
291 }
292
293 public:
294 static Index create_unique(int rank, int nprocs, std::unique_ptr<MeshNodeType>& node, std::vector<int>& ranks)
295 {
296 MeshNodeType* nnode = nullptr;
297 std::vector<Index> cranks, ctags;
298 Index level = Index(create(rank, nprocs, nnode, cranks, ctags));
299 for(auto it = cranks.begin(); it != cranks.end(); ++it)
300 ranks.push_back(int(*it));
301 node.reset(nnode);
302 return level;
303 }
304
305 /******************WIP********************************************
306 * For now an additional wrapper...
307 */
308 static Index create(int rank, int nprocs, std::unique_ptr<MeshNodeType>& node, std::vector<int>& ranks, std::vector<int>& tags)
309 {
310 MeshNodeType* nnode = nullptr;
311 std::vector<Index> cranks, ctags;
312 int lvl = create(rank, nprocs, nnode, cranks, ctags);
313 node = std::unique_ptr<MeshNodeType>(nnode);
314 for(auto it = cranks.begin(); it != cranks.end(); ++it)
315 ranks.push_back(int(*it));
316 for(auto it = ctags.begin(); it != ctags.end(); ++it)
317 tags.push_back(int(*it));
318 return Index(lvl);
319 }
320 /*************************************************************************************/
321
322 static int create(int rank, int nprocs, MeshNodeType*& node, std::vector<Index>& ranks, std::vector<Index>& ctags)
323 {
324 XASSERT(nprocs > 0);
325 XASSERT((rank >= 0) && (rank < nprocs));
326
327 // determine slice count and refinement level
328 int level(0);
329 int n(1);
330 for(; n*n < nprocs; n *= 2, ++level) {}
331 XASSERTM(n*n == nprocs, "number of processes must be a power of 4");
332
333 // decompose rank to (i,j)
334 const int ii(rank % n);
335 const int jj(rank / n);
336
337 // comm tag offsets
338 const Index cto_x0 = Index(0);
339 const Index cto_x1 = cto_x0 + Index((n-1)*(n-1));
340 const Index cto_h = cto_x1 + Index((n-1)*(n-1));
341 const Index cto_v = cto_h + Index(n*(n-1));
342
343 // create root mesh node
344 {
345 RootMeshFactory factory(ii, jj, n);
346 node = new MeshNodeType(factory.make_unique());
347 }
348
349 // lower left neighbor
350 if((ii > 0) && (jj > 0))
351 {
352 ranks.push_back(Index(n*(jj-1) + ii-1));
353 ctags.push_back(cto_x0 + Index((n-1)*(jj-1) + ii-1));
354 node->add_halo(int(ranks.back()), create_halo0(0));
355 }
356 // lower right neighbor
357 if((ii+1 < n) && (jj > 0))
358 {
359 ranks.push_back(Index(n*(jj-1) + ii+1));
360 ctags.push_back(cto_x1 + Index((n-1)*(jj-1) + ii));
361 node->add_halo(int(ranks.back()), create_halo0(1));
362 }
363 // upper left neighbor
364 if((ii > 0) && (jj+1 < n))
365 {
366 ranks.push_back(Index(n*(jj+1) + ii-1));
367 ctags.push_back(cto_x1 + Index((n-1)*jj + ii-1));
368 node->add_halo(int(ranks.back()), create_halo0(2));
369 }
370 // upper right neighbor
371 if((ii+1 < n) && (jj+1 < n))
372 {
373 ranks.push_back(Index(n*(jj+1) + ii+1));
374 ctags.push_back(cto_x0 + Index((n-1)*jj + ii));
375 node->add_halo(int(ranks.back()), create_halo0(3));
376 }
377
378 // bottom neighbor
379 if(jj > 0)
380 {
381 ranks.push_back(Index(n*(jj-1) + ii));
382 ctags.push_back(cto_v + Index(n*(jj-1) + ii));
383 node->add_halo(int(ranks.back()), create_halo1(0));
384 node->add_mesh_part("bnd:0", nullptr);
385 }
386 else
387 {
388 node->add_mesh_part("bnd:0", create_halo1(0));
389 }
390 // top neighbor
391 if(jj+1 < n)
392 {
393 ranks.push_back(Index(n*(jj+1) + ii));
394 ctags.push_back(cto_v + Index(n*jj + ii));
395 node->add_halo(int(ranks.back()), create_halo1(1));
396 node->add_mesh_part("bnd:1", nullptr);
397 }
398 else
399 {
400 node->add_mesh_part("bnd:1", create_halo1(1));
401 }
402 // left neighbor
403 if(ii > 0)
404 {
405 ranks.push_back(Index(n*(jj) + ii-1));
406 ctags.push_back(cto_h + Index((n-1)*jj + ii - 1));
407 node->add_halo(int(ranks.back()), create_halo1(2));
408 node->add_mesh_part("bnd:2", nullptr);
409 }
410 else
411 {
412 node->add_mesh_part("bnd:2", create_halo1(2));
413 }
414 // right neighbor
415 if(ii+1 < n)
416 {
417 ranks.push_back(Index(n*(jj) + ii+1));
418 ctags.push_back(cto_h + Index((n-1)*jj + ii));
419 node->add_halo(int(ranks.back()), create_halo1(3));
420 node->add_mesh_part("bnd:3", nullptr);
421 }
422 else
423 {
424 node->add_mesh_part("bnd:3", create_halo1(3));
425 }
426
427 // return refinement level
428 return level;
429 }
430 }; // UnitCubePatchGenerator<Hypercube<2>>
431
432
433 template<typename Coord_>
434 class UnitCubePatchGenerator<ConformalMesh<Shape::Hypercube<3>, 3, Coord_>>
435 {
436 public:
437 typedef Shape::Hypercube<3> ShapeType;
438 typedef Coord_ CoordType;
439 typedef ConformalMesh<ShapeType, 3, Coord_> MeshType;
440 typedef MeshPart<MeshType> PartType;
441 typedef RootMeshNode<MeshType> MeshNodeType;
442
443 private:
444 class RootMeshFactory :
445 public Geometry::UnitCubeFactory<MeshType>
446 {
447 public:
448 typedef typename MeshType::VertexSetType VertexSetType;
449
450 private:
451 const int _i, _j, _k, _n;
452
453 public:
454 RootMeshFactory(int i, int j, int k, int n) :
455 _i(i), _j(j), _k(k), _n(n)
456 {
457 }
458
459 virtual void fill_vertex_set(VertexSetType& vertex_set) override
460 {
461 vertex_set[0][0] = CoordType(_i ) / CoordType(_n);
462 vertex_set[0][1] = CoordType(_j ) / CoordType(_n);
463 vertex_set[0][2] = CoordType(_k ) / CoordType(_n);
464
465 vertex_set[1][0] = CoordType(_i+1) / CoordType(_n);
466 vertex_set[1][1] = CoordType(_j ) / CoordType(_n);
467 vertex_set[1][2] = CoordType(_k ) / CoordType(_n);
468
469 vertex_set[2][0] = CoordType(_i ) / CoordType(_n);
470 vertex_set[2][1] = CoordType(_j+1) / CoordType(_n);
471 vertex_set[2][2] = CoordType(_k ) / CoordType(_n);
472
473 vertex_set[3][0] = CoordType(_i+1) / CoordType(_n);
474 vertex_set[3][1] = CoordType(_j+1) / CoordType(_n);
475 vertex_set[3][2] = CoordType(_k ) / CoordType(_n);
476
477 vertex_set[4][0] = CoordType(_i ) / CoordType(_n);
478 vertex_set[4][1] = CoordType(_j ) / CoordType(_n);
479 vertex_set[4][2] = CoordType(_k+1) / CoordType(_n);
480
481 vertex_set[5][0] = CoordType(_i+1) / CoordType(_n);
482 vertex_set[5][1] = CoordType(_j ) / CoordType(_n);
483 vertex_set[5][2] = CoordType(_k+1) / CoordType(_n);
484
485 vertex_set[6][0] = CoordType(_i ) / CoordType(_n);
486 vertex_set[6][1] = CoordType(_j+1) / CoordType(_n);
487 vertex_set[6][2] = CoordType(_k+1) / CoordType(_n);
488
489 vertex_set[7][0] = CoordType(_i+1) / CoordType(_n);
490 vertex_set[7][1] = CoordType(_j+1) / CoordType(_n);
491 vertex_set[7][2] = CoordType(_k+1) / CoordType(_n);
492 }
493 };
494
495 class Halo0Factory :
496 public Geometry::Factory<PartType>
497 {
498 public:
499 typedef Geometry::Factory<PartType> BaseClass;
500 typedef typename BaseClass::AttributeSetContainer AttributeSetContainer;
501 typedef typename BaseClass::IndexSetHolderType IndexSetHolderType;
502 typedef typename BaseClass::TargetSetHolderType TargetSetHolderType;
503
504 private:
505 const int _p;
506
507 public:
508 explicit Halo0Factory(int p) :
509 _p(p)
510 {
511 }
512
513 virtual Index get_num_entities(int dim) override
514 {
515 return Index(dim == 0 ? 1 : 0);
516 }
517
518 virtual void fill_attribute_sets(AttributeSetContainer&) override
519 {
520 }
521
522 virtual void fill_index_sets(std::unique_ptr<IndexSetHolderType>&) override
523 {
524 }
525
526 virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
527 {
528 target_set_holder.template get_target_set<0>()[0] = Index(_p);
529 }
530 };
531
532 class Halo1Factory :
533 public Geometry::Factory<PartType>
534 {
535 public:
536 typedef Geometry::Factory<PartType> BaseClass;
537 typedef typename BaseClass::AttributeSetContainer AttributeSetContainer;
538 typedef typename BaseClass::IndexSetHolderType IndexSetHolderType;
539 typedef typename BaseClass::TargetSetHolderType TargetSetHolderType;
540
541 private:
542 const int _p;
543
544 public:
545 explicit Halo1Factory(int p) :
546 _p(p)
547 {
548 }
549
550 virtual Index get_num_entities(int dim) override
551 {
552 return Index(dim == 0 ? 2 : (dim == 1 ? 1 : 0));
553 }
554
555 virtual void fill_attribute_sets(AttributeSetContainer&) override
556 {
557 }
558
559 virtual void fill_index_sets(std::unique_ptr<IndexSetHolderType>&) override
560 {
561 }
562
563 virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
564 {
565 typedef Geometry::Intern::FaceIndexMapping<ShapeType, 1, 0> Fim;
566 target_set_holder.template get_target_set<0>()[0] = Index(Fim::map(_p, 0));
567 target_set_holder.template get_target_set<0>()[1] = Index(Fim::map(_p, 1));
568 target_set_holder.template get_target_set<1>()[0] = Index(_p);
569 }
570 };
571
572 class Halo2Factory :
573 public Geometry::Factory<PartType>
574 {
575 public:
576 typedef Geometry::Factory<PartType> BaseClass;
577 typedef typename BaseClass::AttributeSetContainer AttributeSetContainer;
578 typedef typename BaseClass::IndexSetHolderType IndexSetHolderType;
579 typedef typename BaseClass::TargetSetHolderType TargetSetHolderType;
580
581 private:
582 const int _p;
583
584 public:
585 explicit Halo2Factory(int p) :
586 _p(p)
587 {
588 }
589
590 virtual Index get_num_entities(int dim) override
591 {
592 return Index(dim == 0 ? 4 : (dim == 1 ? 4 : (dim == 2 ? 1 : 0)));
593 }
594
595 virtual void fill_attribute_sets(AttributeSetContainer&) override
596 {
597 }
598
599 virtual void fill_index_sets(std::unique_ptr<IndexSetHolderType>&) override
600 {
601 }
602
603 virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
604 {
605 typedef Geometry::Intern::FaceIndexMapping<ShapeType, 2, 0> Fim0;
606 target_set_holder.template get_target_set<0>()[0] = Index(Fim0::map(_p, 0));
607 target_set_holder.template get_target_set<0>()[1] = Index(Fim0::map(_p, 1));
608 target_set_holder.template get_target_set<0>()[2] = Index(Fim0::map(_p, 2));
609 target_set_holder.template get_target_set<0>()[3] = Index(Fim0::map(_p, 3));
610 typedef Geometry::Intern::FaceIndexMapping<ShapeType, 2, 1> Fim1;
611 target_set_holder.template get_target_set<1>()[0] = Index(Fim1::map(_p, 0));
612 target_set_holder.template get_target_set<1>()[1] = Index(Fim1::map(_p, 1));
613 target_set_holder.template get_target_set<1>()[2] = Index(Fim1::map(_p, 2));
614 target_set_holder.template get_target_set<1>()[3] = Index(Fim1::map(_p, 3));
615 target_set_holder.template get_target_set<2>()[0] = Index(_p);
616 }
617 };
618
619 static std::unique_ptr<PartType> create_halo0(int k)
620 {
621 Halo0Factory factory(k);
622 return factory.make_unique();
623 }
624
625 static std::unique_ptr<PartType> create_halo1(int k)
626 {
627 Halo1Factory factory(k);
628 return factory.make_unique();
629 }
630
631 static std::unique_ptr<PartType> create_halo2(int k)
632 {
633 Halo2Factory factory(k);
634 return factory.make_unique();
635 }
636
637
638 public:
639 static Index create_unique(int rank, int nprocs, std::unique_ptr<MeshNodeType>& node, std::vector<int>& ranks)
640 {
641 MeshNodeType* nnode = nullptr;
642 std::vector<Index> cranks, ctags;
643 Index level = Index(create(rank, nprocs, nnode, cranks, ctags));
644 for(auto it = cranks.begin(); it != cranks.end(); ++it)
645 ranks.push_back(int(*it));
646 node.reset(nnode);
647 return level;
648 }
649
650
651 static int create(int rank, int nprocs, MeshNodeType*& node, std::vector<Index>& ranks, std::vector<Index>& ctags)
652 {
653 XASSERT(nprocs > 0);
654 XASSERT((rank >= 0) && (rank < nprocs));
655
656 // determine slice count and refinement level
657 int level(0);
658 int n(1);
659 for(; n*n*n < nprocs; n *= 2, ++level) {}
660 XASSERTM(n*n*n == nprocs, "number of processes must be a power of 8");
661
662 // decompose rank to (i,j,k)
663 const int ii(rank % n);
664 const int jj((rank / n) % n);
665 const int kk(rank / (n*n));
666
667 // comm tag offsets
668 const Index cto_x0 = Index(0);
669 const Index cto_x1 = cto_x0 + Index((n-1)*(n-1)*(n-1));
670 const Index cto_x2 = cto_x1 + Index((n-1)*(n-1)*(n-1));
671 const Index cto_x3 = cto_x2 + Index((n-1)*(n-1)*(n-1));
672
673 const Index cto_h0 = cto_x3 + Index((n-1)*(n-1)*(n-1));
674 const Index cto_h1 = cto_h0 + Index(n*(n-1)*(n-1));
675
676 const Index cto_v0 = cto_h1 + Index(n*(n-1)*(n-1));
677 const Index cto_v1 = cto_v0 + Index(n*(n-1)*(n-1));
678
679 const Index cto_z0 = cto_v1 + Index(n*(n-1)*(n-1));
680 const Index cto_z1 = cto_z0 + Index(n*(n-1)*(n-1));
681
682 const Index cto_f0 = cto_z1 + Index(n*(n-1)*(n-1));
683 const Index cto_f1 = cto_f0 + Index((n-1)*n*n);
684 const Index cto_f2 = cto_f1 + Index((n-1)*n*n);
685
686 // create root mesh
687 {
688 RootMeshFactory factory(ii, jj, kk, n);
689 node = new MeshNodeType(factory.make_unique());
690 }
691
692 // lower/upper -> kk
693 // south/north -> jj
694 // west /east -> ii
695 // lower south west neighbor
696 if((ii > 0) && (jj > 0) && (kk > 0))
697 {
698 ranks.push_back(Index(n*n*(kk-1) + n*(jj-1) + ii-1));
699 ctags.push_back(cto_x0 + Index((n-1)*(n-1)*(kk-1) + (n-1)*(jj-1) + ii-1));
700 node->add_halo(int(ranks.back()), create_halo0(0));
701 }
702 // lower south east neighbor
703 if((ii+1 < n) && (jj > 0) && (kk > 0))
704 {
705 ranks.push_back(Index(n*n*(kk-1) + n*(jj-1) + ii+1));
706 ctags.push_back(cto_x1 + Index((n-1)*(n-1)*(kk-1) + (n-1)*(jj-1) + ii));
707 node->add_halo(int(ranks.back()), create_halo0(1));
708 }
709 // lower north west neighbor
710 if((ii > 0) && (jj+1 < n) && (kk > 0))
711 {
712 ranks.push_back(Index(n*n*(kk-1) + n*(jj+1) + ii-1));
713 ctags.push_back(cto_x2 + Index((n-1)*(n-1)*(kk-1) + (n-1)*jj + ii-1));
714 node->add_halo(int(ranks.back()), create_halo0(2));
715 }
716 // lower north east neighbor
717 if((ii+1 < n) && (jj+1 < n) && (kk > 0))
718 {
719 ranks.push_back(Index(n*n*(kk-1) + n*(jj+1) + ii+1));
720 ctags.push_back(cto_x3 + Index((n-1)*(n-1)*(kk-1) + (n-1)*jj + ii));
721 node->add_halo(int(ranks.back()), create_halo0(3));
722 }
723
724 // upper south west
725 if((ii > 0) && (jj > 0) && (kk+1 < n))
726 {
727 ranks.push_back(Index(n*n*(kk+1) + n*(jj-1) + ii-1));
728 ctags.push_back(cto_x3 + Index((n-1)*(n-1)*(kk) + (n-1)*(jj-1) + ii-1));
729 node->add_halo(int(ranks.back()), create_halo0(4));
730 }
731 // upper south east neighbor
732 if((ii+1 < n) && (jj > 0) && (kk+1 < n))
733 {
734 ranks.push_back(Index(n*n*(kk+1) + n*(jj-1) + ii+1));
735 ctags.push_back(cto_x2 + Index((n-1)*(n-1)*(kk) + (n-1)*(jj-1) + ii));
736 node->add_halo(int(ranks.back()), create_halo0(5));
737 }
738 // upper north west neighbor
739 if((ii > 0) && (jj+1 < n) && (kk+1 < n))
740 {
741 ranks.push_back(Index(n*n*(kk+1) + n*(jj+1) + ii-1));
742 ctags.push_back(cto_x1 + Index((n-1)*(n-1)*(kk) + (n-1)*jj + ii-1));
743 node->add_halo(int(ranks.back()), create_halo0(6));
744 }
745 // upper north east neighbor
746 if((ii+1 < n) && (jj+1 < n) && (kk+1 < n))
747 {
748 ranks.push_back(Index(n*n*(kk+1) + n*(jj+1) + ii+1));
749 ctags.push_back(cto_x0 + Index((n-1)*(n-1)*(kk) + (n-1)*jj + ii));
750 node->add_halo(int(ranks.back()), create_halo0(7));
751 }
752
753 // lower south neighbor
754 if((jj > 0) && (kk > 0))
755 {
756 ranks.push_back(Index(n*n*(kk-1) + n*(jj-1) + ii));
757 ctags.push_back(cto_h0 + Index((n-1)*n*(kk-1) + n*(jj-1) + ii));
758 node->add_halo(int(ranks.back()), create_halo1(0));
759 }
760 // lower north neighbor
761 if((jj + 1 < n) && (kk > 0))
762 {
763 ranks.push_back(Index(n*n*(kk-1) + n*(jj+1) + ii));
764 ctags.push_back(cto_h1 + Index((n-1)*n*(kk-1) + n*jj + ii));
765 node->add_halo(int(ranks.back()), create_halo1(1));
766 }
767 // upper south neighbor
768 if((jj > 0) && (kk + 1 < n))
769 {
770 ranks.push_back(Index(n*n*(kk+1) + n*(jj-1) + ii));
771 ctags.push_back(cto_h1 + Index((n-1)*n*(kk) + n*(jj-1) + ii));
772 node->add_halo(int(ranks.back()), create_halo1(2));
773 }
774 // upper north neighbor
775 if((jj + 1 < n) && (kk + 1 < n))
776 {
777 ranks.push_back(Index(n*n*(kk+1) + n*(jj+1) + ii));
778 ctags.push_back(cto_h0 + Index((n-1)*n*(kk) + n*jj + ii));
779 node->add_halo(int(ranks.back()), create_halo1(3));
780 }
781
782 // lower west neighbor
783 if((ii > 0) && (kk > 0))
784 {
785 ranks.push_back(Index(n*n*(kk-1) + n*jj + ii-1));
786 ctags.push_back(cto_v0 + Index((n-1)*n*(kk-1) + n*jj + ii-1));
787 node->add_halo(int(ranks.back()), create_halo1(4));
788 }
789 // lower east neighbor
790 if((ii + 1 < n) && (kk > 0))
791 {
792 ranks.push_back(Index(n*n*(kk-1) + n*jj + ii+1));
793 ctags.push_back(cto_v1 + Index((n-1)*n*(kk-1) + n*jj + ii));
794 node->add_halo(int(ranks.back()), create_halo1(5));
795 }
796 // upper west neighbor
797 if((ii > 0) && (kk + 1 < n))
798 {
799 ranks.push_back(Index(n*n*(kk+1) + n*jj + ii-1));
800 ctags.push_back(cto_v1 + Index((n-1)*n*(kk) + n*jj + ii-1));
801 node->add_halo(int(ranks.back()), create_halo1(6));
802 }
803 // upper east neighbor
804 if((ii + 1 < n) && (kk + 1 < n))
805 {
806 ranks.push_back(Index(n*n*(kk+1) + n*jj + ii+1));
807 ctags.push_back(cto_v0 + Index((n-1)*n*(kk) + n*jj + ii));
808 node->add_halo(int(ranks.back()), create_halo1(7));
809 }
810
811 // south west neighbor
812 if((ii > 0) && (jj > 0))
813 {
814 ranks.push_back(Index(n*n*kk + n*(jj-1) + ii-1));
815 ctags.push_back(cto_z0 + Index((n-1)*(n-1)*(kk) + (n-1)*(jj-1) + ii-1));
816 node->add_halo(int(ranks.back()), create_halo1(8));
817 }
818 // south east neighbor
819 if((ii + 1 < n) && (jj > 0))
820 {
821 ranks.push_back(Index(n*n*kk + n*(jj-1) + ii+1));
822 ctags.push_back(cto_z1 + Index((n-1)*(n-1)*(kk) + (n-1)*(jj-1) + ii));
823 node->add_halo(int(ranks.back()), create_halo1(9));
824 }
825 // north west neighbor
826 if((ii > 0) && (jj + 1 < n))
827 {
828 ranks.push_back(Index(n*n*kk + n*(jj+1) + ii-1));
829 ctags.push_back(cto_z1 + Index((n-1)*(n-1)*(kk) + (n-1)*jj + ii-1));
830 node->add_halo(int(ranks.back()), create_halo1(10));
831 }
832 // north east neighbor
833 if((ii + 1 < n) && (jj + 1 < n))
834 {
835 ranks.push_back(Index(n*n*kk + n*(jj+1) + ii+1));
836 ctags.push_back(cto_z0 + Index((n-1)*(n-1)*(kk) + (n-1)*jj + ii));
837 node->add_halo(int(ranks.back()), create_halo1(11));
838 }
839
840 // lower neighbor
841 if(kk > 0)
842 {
843 ranks.push_back(Index(n*n*(kk-1) + n*jj + ii));
844 ctags.push_back(cto_f0 + Index(n*n*(kk-1) + jj*n + ii));
845 node->add_halo(int(ranks.back()), create_halo2(0));
846 node->add_mesh_part("bnd:0", nullptr);
847 }
848 else
849 node->add_mesh_part("bnd:0", create_halo2(0));
850 // upper neighbor
851 if(kk + 1 < n)
852 {
853 ranks.push_back(Index(n*n*(kk+1) + n*jj + ii));
854 ctags.push_back(cto_f0 + Index(n*n*kk + jj*n + ii));
855 node->add_halo(int(ranks.back()), create_halo2(1));
856 node->add_mesh_part("bnd:1", nullptr);
857 }
858 else
859 node->add_mesh_part("bnd:1", create_halo2(1));
860 // south neighbor
861 if(jj > 0)
862 {
863 ranks.push_back(Index(n*n*kk + n*(jj-1) + ii));
864 ctags.push_back(cto_f1 + Index(n*n*kk + (jj-1)*n + ii));
865 node->add_halo(int(ranks.back()), create_halo2(2));
866 node->add_mesh_part("bnd:2", nullptr);
867 }
868 else
869 node->add_mesh_part("bnd:2", create_halo2(2));
870 // north neighbor
871 if(jj + 1 < n)
872 {
873 ranks.push_back(Index(n*n*kk + n*(jj+1) + ii));
874 ctags.push_back(cto_f1 + Index(n*n*kk + jj*n + ii));
875 node->add_halo(int(ranks.back()), create_halo2(3));
876 node->add_mesh_part("bnd:3", nullptr);
877 }
878 else
879 node->add_mesh_part("bnd:3", create_halo2(3));
880 // west neighbor
881 if(ii > 0)
882 {
883 ranks.push_back(Index(n*n*kk + n*jj + ii-1));
884 ctags.push_back(cto_f2 + Index(n*n*kk + jj*n + ii-1));
885 node->add_halo(int(ranks.back()), create_halo2(4));
886 node->add_mesh_part("bnd:4", nullptr);
887 }
888 else
889 node->add_mesh_part("bnd:4", create_halo2(4));
890 // east neighbor
891 if(ii + 1 < n)
892 {
893 ranks.push_back(Index(n*n*kk + n*jj + ii+1));
894 ctags.push_back(cto_f2 + Index(n*n*kk + jj*n + ii));
895 node->add_halo(int(ranks.back()), create_halo2(5));
896 node->add_mesh_part("bnd:5", nullptr);
897 }
898 else
899 node->add_mesh_part("bnd:5", create_halo2(5));
900
901 return level;
902 }
903 }; // UnitCubePatchGenerator<Hypercube<3>>
905 } // namespace Geometry
906} // namespace FEAT
XASSERT
#define XASSERT(expr)
Assertion macro definition.
Definition: assertion.hpp:262
XASSERTM
#define XASSERTM(expr, msg)
Assertion macro definition with custom message.
Definition: assertion.hpp:263
FEAT::Geometry::ConformalMesh
Conformal mesh class template.
Definition: conformal_mesh.hpp:43
FEAT::Geometry::MeshPart
Class template for partial meshes.
Definition: mesh_part.hpp:90
FEAT::Geometry::RootMeshNode
Root mesh node class template.
Definition: mesh_node.hpp:748
FEAT
FEAT namespace.
Definition: adjactor.hpp:12
FEAT::Index
std::uint64_t Index
Index data type.
Definition: base_header.hpp:122
FEAT::Shape::Hypercube
Hypercube shape tag struct template.
Definition: shape.hpp:64