Dune-Fufem 2.11-git
Loading...
Searching...
No Matches
makering.hh
Go to the documentation of this file.
1// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set ts=8 sw=4 et sts=4:
3
4// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
5// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
6
7#ifndef MAKE_RING_HH
8#define MAKE_RING_HH
9
10#include <memory>
11#include <vector>
12
13#include <dune/grid/common/gridfactory.hh>
14#include <dune/grid/common/boundarysegment.hh>
15
16#include <dune/fufem/arcofcircle.hh>
17#include <dune/fufem/boundarypatch.hh>
18
19template<class GridType>
20void makeRingSegment2D(const Dune::FieldVector<typename GridType::ctype,2>& center,
21 const typename GridType::ctype thickness, const typename GridType::ctype innerRadius,
22 const typename GridType::ctype fromAngle, const typename GridType::ctype toAngle,
23 GridType& grid, size_t nElement, bool closeRing, bool quadrilateral = false,
24 BoundaryPatch<typename GridType::LevelGridView>* innerBoundary = NULL,
25 BoundaryPatch<typename GridType::LevelGridView>* outerBoundary = NULL)
26{
27 DUNE_THROW(Dune::Exception, "This method does not work anymore, please use the one returning a unique_ptr");
28}
29
45template<class GridType>
46std::unique_ptr<GridType> makeRingSegment2D(const Dune::FieldVector<typename GridType::ctype,2>& center,
47 const typename GridType::ctype thickness, const typename GridType::ctype innerRadius,
48 const typename GridType::ctype fromAngle, const typename GridType::ctype toAngle,
49 size_t nElement, bool closeRing, bool quadrilateral = false,
50 BoundaryPatch<typename GridType::LevelGridView>* innerBoundary = NULL,
51 BoundaryPatch<typename GridType::LevelGridView>* outerBoundary = NULL)
52{
53 static_assert(GridType::dimension==2, "makeRingSegment2D can only be called for 2-dimensional grids");
54 typedef typename GridType::ctype ctype;
55
56 ctype eps = 1e-10;
57 if ( (-eps >= fromAngle) || (fromAngle > toAngle) || (toAngle > 2*M_PI) )
58 DUNE_THROW(Dune::Exception, "Parameters do not fulfill 0 <= fromAngle <= toAngle <= 2*M_PI");
59
60 bool fullRing = closeRing && (fromAngle==0.0) && (toAngle == 2*M_PI);
61
62 using namespace Dune;
63
64 Dune::GridFactory<GridType> factory;
65
66 // ///////////////////////
67 // Insert vertices
68 // ///////////////////////
69
70 ArcOfCircle innerArc(center, innerRadius, fromAngle, toAngle);
71
72 ctype step =1.0/((ctype) nElement);
73
74 // insert inner ring vertices
75 for (size_t i=0; i<nElement; i++)
76 factory.insertVertex(innerArc(i*step));
77
78 // if the ring is not closed also make vertices for the right boundary
79 if (!fullRing)
80 factory.insertVertex(innerArc(1.0));
81
82 ArcOfCircle outerArc(center, innerRadius+thickness, fromAngle, toAngle);
83
84 // insert outer ring vertices
85 for (size_t i=0; i<nElement; i++)
86 factory.insertVertex(outerArc(i*step));
87
88 // if the ring is not closed also make vertices for the right boundary
89 if (!fullRing)
90 factory.insertVertex(outerArc(1.0));
91
92 // /////////////////
93 // Insert elements
94 // /////////////////
95
96 std::vector<unsigned int> cornerIDs((quadrilateral) ?4 : 3);
97 GeometryType gt = (quadrilateral) ? GeometryTypes::quadrilateral : GeometryTypes::triangle;
98
99 // if ring is closed last nodes are the first nodes
100 if (!fullRing) {
101
102 if (quadrilateral) {
103
104 for (size_t i=0; i<nElement; i++) {
105 cornerIDs[0] = i+1;
106 cornerIDs[1] = i;
107 cornerIDs[2] = i+nElement+2;
108 cornerIDs[3] = i+nElement+1;
109 factory.insertElement(gt, cornerIDs);
110 }
111 } else {
112
113 for (size_t i=0; i<nElement; i++) {
114 // first triangle
115 cornerIDs[0] = i+1;
116 cornerIDs[1] = i;
117 cornerIDs[2] = i+nElement+2;
118 factory.insertElement(gt, cornerIDs);
119
120 // second
121 cornerIDs[0] = nElement +i +1;
122 cornerIDs[1] = nElement +i +2;
123 cornerIDs[2] = i;
124 factory.insertElement(gt, cornerIDs);
125 }
126 }
127 } else {
128
129
130 if (quadrilateral) {
131
132 for (size_t i=0; i<nElement-1; i++) {
133 cornerIDs[0] = i+1;
134 cornerIDs[1] = i;
135 cornerIDs[2] = i+nElement+1;
136 cornerIDs[3] = i+nElement;
137 factory.insertElement(gt, cornerIDs);
138 }
139 cornerIDs[0] = 0;
140 cornerIDs[1] = nElement-1;
141 cornerIDs[2] = nElement;
142 cornerIDs[3] = 2*nElement -1;
143 factory.insertElement(gt, cornerIDs);
144
145 } else {
146
147 for (size_t i=0; i<nElement-1; i++) {
148 // first triangle
149 cornerIDs[0] = i+1;
150 cornerIDs[1] = i;
151 cornerIDs[2] = i+nElement+1;
152 factory.insertElement(gt, cornerIDs);
153
154 // second
155 cornerIDs[0] = nElement +i;
156 cornerIDs[1] = nElement +i +1;
157 cornerIDs[2] = i;
158 factory.insertElement(gt, cornerIDs);
159 }
160 // first triangle
161 cornerIDs[0] = 0;
162 cornerIDs[1] = nElement-1;
163 cornerIDs[2] = nElement;
164 factory.insertElement(gt, cornerIDs);
165
166 // second triangle
167 cornerIDs[0] = 2*nElement-1;
168 cornerIDs[1] = nElement;
169 cornerIDs[2] = nElement -1;
170 factory.insertElement(gt, cornerIDs);
171
172 }
173 }
174
175 // /////////////////////////////
176 // Create boundary segments
177 // /////////////////////////////
178 std::vector<unsigned int> vertices(2);
179
180 typedef typename std::shared_ptr<Dune::BoundarySegment<2> > SegmentPointer;
181 ctype angleDiff = toAngle-fromAngle;
182
183 for (size_t i=0; i<nElement-1; i++) {
184
185 vertices[0] = i; vertices[1] = i+1;
186 factory.insertBoundarySegment(vertices,
187 SegmentPointer(new ArcOfCircle(center, innerRadius, fromAngle+i*step*angleDiff,
188 fromAngle+(i+1)*step*angleDiff)));
189
190 vertices[0] = nElement+i+(!fullRing); vertices[1] = nElement +i+1+(!fullRing);
191 factory.insertBoundarySegment(vertices,
192 SegmentPointer(new ArcOfCircle(center, innerRadius+thickness, fromAngle+i*step*angleDiff,
193 fromAngle+(i+1)*step*angleDiff)));
194
195 }
196
197 // if ring is closed last nodes are the first nodes
198
199 vertices[0] = nElement-1; vertices[1] = (!fullRing)*nElement;
200 factory.insertBoundarySegment(vertices,
201 SegmentPointer(new ArcOfCircle(center, innerRadius, fromAngle+(nElement-1)*step*angleDiff,
202 toAngle)));
203
204 vertices[0] = 2*nElement-1+(!fullRing); vertices[1] =nElement + (!fullRing)*(nElement+1);
205 factory.insertBoundarySegment(vertices,
206 SegmentPointer(new ArcOfCircle(center, innerRadius+thickness, fromAngle+(nElement-1)*step*angleDiff,
207 toAngle)));
208
209 // //////////////////////////////////////
210 // Finish initialization
211 // //////////////////////////////////////
212
213 auto grid = factory.createGrid();
214
215
216 // create boundarypatches if desired
217 if (innerBoundary) {
218
219 BitSetVector<1> innerVert(2*(nElement+!(fullRing)));
220 size_t half = innerVert.size()/2;
221
222 for (size_t i=0; i<half; i++)
223 innerVert[i] = true;
224
225 innerBoundary->setup(grid->levelGridView(0),innerVert);
226 }
227
228 if (outerBoundary) {
229
230 BitSetVector<1> outerVert(2*(nElement+!(fullRing)));
231 size_t half = outerVert.size()/2;
232
233 for (size_t i=0; i<half; i++)
234 outerVert[half+i] = true;
235
236 outerBoundary->setup(grid->levelGridView(0),outerVert);
237 }
238
239 return grid;
240}
241
242
267template <class GridType>
268std::unique_ptr<GridType> makeTubeSegment3D(
269 const Dune::FieldVector<typename GridType::ctype, 3>& center,
270 const typename GridType::ctype thickness,
271 const typename GridType::ctype length,
272 const typename GridType::ctype innerRadius,
273 const typename GridType::ctype fromAngle,
274 const typename GridType::ctype toAngle, unsigned int nElementRing,
275 unsigned int nElementLength, bool closeTube, unsigned int axis = 0,
276 bool tetrahedra = true, bool parametrizedBoundary = true,
277 BoundaryPatch<typename GridType::LevelGridView> * innerBoundary = NULL,
278 BoundaryPatch<typename GridType::LevelGridView> * outerBoundary = NULL)
279{
280 static_assert(GridType::dimension==3, "makeTubeSegment3D can only be called for 3-dimensional grids");
281 typedef typename GridType::ctype ctype;
282
283 if (fromAngle >= toAngle)
284 DUNE_THROW(Dune::Exception, "Parameters do not fulfill fromAngle < toAngle");
285
286 bool fullTube = closeTube && (fromAngle==0.0) && (toAngle == 2*M_PI);
287
288 using namespace Dune;
289
290 Dune::GridFactory<GridType> factory;
291
292 // ///////////////////////
293 // Insert vertices
294 // ///////////////////////
295
296 unsigned int nRingNodes = nElementRing+!fullTube;
297
298 // project center on the yz-plane
299 FieldVector<ctype,2> localCenter;
300 localCenter[0] = center[(axis+1)%3];
301 localCenter[1] = center[(axis+2)%3];
302
303 ArcOfCircle innerArc(localCenter, innerRadius, fromAngle, toAngle);
304
305 ctype stepRing =1.0/((ctype) nElementRing);
306 ctype stepLength = length/((ctype) nElementLength);
307
308 // insert tube vertices
309 for (size_t j=0; j<=nElementLength; j++) {
310
311 FieldVector<ctype,3> vertex(0);
312
313 // x-coordinate is determined by the center x-coord
314 vertex[axis] = center[axis] + j*stepLength;
315
316 for (size_t i=0; i<nRingNodes; i++) {
317
318 FieldVector<ctype,2> local = innerArc(i*stepRing);
319 vertex[(axis+1)%3] = local[0]; vertex[(axis+2)%3] = local[1];
320
321 factory.insertVertex(vertex);
322 }
323 }
324
325 ArcOfCircle outerArc(localCenter, innerRadius+thickness, fromAngle, toAngle);
326
327 // insert outer tube vertices
328 for (unsigned int j=0; j<=nElementLength; j++) {
329
330 FieldVector<ctype,3> vertex(0);
331
332 // x-coordinate is determined by the center x-coord
333 vertex[axis] = center[axis] + j*stepLength;
334
335 for (size_t i=0; i<nRingNodes; i++) {
336
337 FieldVector<ctype,2> local = outerArc(i*stepRing);
338 vertex[(axis+1)%3] = local[0]; vertex[(axis+2)%3] = local[1];
339
340 factory.insertVertex(vertex);
341 }
342 }
343
344
345 unsigned int innerNodes = (nElementLength+1)*(nRingNodes);
346
347 // /////////////////
348 // Insert elements
349 // /////////////////
350
351 // if ring is closed last nodes are the first nodes
352
353 std::vector<unsigned int> cornerIDs(4+!tetrahedra*2);
354
355 Dune::GeometryType gt = Dune::GeometryTypes::tetrahedron;
356 if (!tetrahedra)
357 gt = Dune::GeometryTypes::prism;
358
359 for (unsigned int i=0; i<nElementLength; i++) {
360
361 unsigned int row = i*nRingNodes;
362
363 for (unsigned int j=0; j<nRingNodes-1; j++) {
364
365 if (tetrahedra) {
366
367 std::vector<unsigned int> hexa = {row+j,row+j+1,
368 row+j+1+nElementRing,row+j+2+nElementRing,
369 row+j+innerNodes,row+j+1+innerNodes,
370 row+j+1+nElementRing+innerNodes,
371 row+j+2+nElementRing+innerNodes};
372 //add 6 tetrahedra
373 cornerIDs={hexa[0],hexa[1],hexa[3],hexa[7]};
374 factory.insertElement(gt,cornerIDs);
375 cornerIDs={hexa[0],hexa[5],hexa[1],hexa[7]};
376 factory.insertElement(gt,cornerIDs);
377 cornerIDs={hexa[0],hexa[4],hexa[5],hexa[7]};
378 factory.insertElement(gt,cornerIDs);
379 cornerIDs={hexa[0],hexa[6],hexa[4],hexa[7]};
380 factory.insertElement(gt,cornerIDs);
381 cornerIDs={hexa[0],hexa[2],hexa[6],hexa[7]};
382 factory.insertElement(gt,cornerIDs);
383 cornerIDs={hexa[0],hexa[3],hexa[2],hexa[7]};
384 factory.insertElement(gt,cornerIDs);
385
386 } else {
387 // Insert a prism
388 cornerIDs[0] = row+j;
389 cornerIDs[1] = row+j+1 ;
390 cornerIDs[2] = row+nElementRing+j+2;
391 cornerIDs[3] = cornerIDs[0]+innerNodes;
392 cornerIDs[4] = cornerIDs[1]+innerNodes;
393 cornerIDs[5] = cornerIDs[2]+innerNodes;
394
395 factory.insertElement(gt, cornerIDs);
396
397 // Insert another prism
398 cornerIDs[0] = row+nElementRing+j+2;
399 cornerIDs[1] = row+nElementRing+j+1;
400 cornerIDs[2] = row+j;
401 cornerIDs[3] = cornerIDs[0]+innerNodes;
402 cornerIDs[4] = cornerIDs[1]+innerNodes;
403 cornerIDs[5] = cornerIDs[2]+innerNodes;
404
405 factory.insertElement(gt, cornerIDs);
406 }
407 }
408
409 if (fullTube) {
410 if (tetrahedra) {
411
412 std::vector<unsigned int> hexa = {row+nElementRing-1,row,
413 row+2*nElementRing-1,row+nElementRing,
414 row+nElementRing-1+innerNodes,
415 row+innerNodes,
416 row+2*nElementRing-1+innerNodes,
417 row+nElementRing+innerNodes};
418 //add 6 tetrahedra
419 cornerIDs={hexa[0],hexa[1],hexa[3],hexa[7]};
420 factory.insertElement(gt,cornerIDs);
421 cornerIDs={hexa[0],hexa[5],hexa[1],hexa[7]};
422 factory.insertElement(gt,cornerIDs);
423 cornerIDs={hexa[0],hexa[4],hexa[5],hexa[7]};
424 factory.insertElement(gt,cornerIDs);
425 cornerIDs={hexa[0],hexa[6],hexa[4],hexa[7]};
426 factory.insertElement(gt,cornerIDs);
427 cornerIDs={hexa[0],hexa[2],hexa[6],hexa[7]};
428 factory.insertElement(gt,cornerIDs);
429 cornerIDs={hexa[0],hexa[3],hexa[2],hexa[7]};
430 factory.insertElement(gt,cornerIDs);
431
432 } else {
433 // Insert a prism
434 cornerIDs[0] = row+nElementRing-1;
435 cornerIDs[1] = row;
436 cornerIDs[2] = row+nElementRing;
437 cornerIDs[3] = cornerIDs[0]+innerNodes;
438 cornerIDs[4] = cornerIDs[1]+innerNodes;
439 cornerIDs[5] = cornerIDs[2]+innerNodes;
440
441 factory.insertElement(gt, cornerIDs);
442
443 // Insert another prism
444 cornerIDs[0] = row+nElementRing;
445 cornerIDs[1] = row+2*nElementRing-1;
446 cornerIDs[2] = row+nElementRing-1;
447 cornerIDs[3] = cornerIDs[0]+innerNodes;
448 cornerIDs[4] = cornerIDs[1]+innerNodes;
449 cornerIDs[5] = cornerIDs[2]+innerNodes;
450
451 factory.insertElement(gt, cornerIDs);
452 }
453 }
454 }
455
456 // /////////////////////////////
457 // Create boundary segments
458 // /////////////////////////////
459
460 if (parametrizedBoundary) {
461 std::vector<unsigned int> vertices(3);
462
463 typedef typename std::shared_ptr<Dune::BoundarySegment<3> > SegmentPointer;
464 ctype angleDiff = toAngle-fromAngle;
465
466 for (size_t i=0; i<nElementLength; i++) {
467
468 unsigned int row = i*nRingNodes;
469
470 FieldVector<ctype,3> localCenter = center;
471 localCenter[0] += i*stepLength;
472
473 // lower triangular faces
474 for (size_t j=0; j<nElementRing-1; j++) {
475
476 vertices[0] = row+j; vertices[1] = row+j+1; vertices[2] = row+nRingNodes+j+1;
477 factory.insertBoundarySegment(vertices,
478 SegmentPointer(new TubeLowerSegment<ctype>(localCenter, innerRadius, stepLength,
479 fromAngle+j*stepRing*angleDiff, fromAngle+(j+1)*stepRing*angleDiff)));
480
481
482 vertices[0] += innerNodes; vertices[1] += innerNodes; vertices[2] +=innerNodes;
483 factory.insertBoundarySegment(vertices,
484 SegmentPointer(new TubeLowerSegment<ctype>(localCenter, innerRadius+thickness, stepLength,
485 fromAngle+j*stepRing*angleDiff, fromAngle+(j+1)*stepRing*angleDiff)));
486 }
487
488 // if ring is closed last nodes are the first nodes
489 vertices[0] = row+nElementRing-1;
490 vertices[1] = row+(!fullTube)*nElementRing;
491 vertices[2] = row+(!fullTube)*(nElementRing+1) + nElementRing;
492
493 factory.insertBoundarySegment(vertices,
494 SegmentPointer(new TubeLowerSegment<ctype>(localCenter, innerRadius, stepLength,
495 fromAngle+(nElementRing-1)*stepRing*angleDiff, toAngle)));
496
497 vertices[0] += innerNodes;
498 vertices[1] += innerNodes;
499 vertices[2] += innerNodes;
500
501 factory.insertBoundarySegment(vertices,
502 SegmentPointer(new TubeLowerSegment<ctype>(localCenter, innerRadius+thickness, stepLength,
503 fromAngle+(nElementRing-1)*stepRing*angleDiff, toAngle)));
504
505
506 // upper triangular faces
507 localCenter[0] += stepLength;
508
509 for (size_t j=0; j<nElementRing-1; j++) {
510
511 vertices[0] = row+nRingNodes+j; vertices[1] = row+nRingNodes+j+1; vertices[2] = row+j;
512 factory.insertBoundarySegment(vertices,
513 SegmentPointer(new TubeUpperSegment<ctype>(localCenter, innerRadius, -stepLength,
514 fromAngle+j*stepRing*angleDiff, fromAngle+(j+1)*stepRing*angleDiff)));
515
516 vertices[0] += innerNodes; vertices[1] += innerNodes; vertices[2] += innerNodes;
517 factory.insertBoundarySegment(vertices,
518 SegmentPointer(new TubeUpperSegment<ctype>(localCenter, innerRadius+thickness, -stepLength,
519 fromAngle+j*stepRing*angleDiff, fromAngle+(j+1)*stepRing*angleDiff)));
520 }
521
522 // if ring is closed last nodes are the first nodes
523 vertices[0] = row+nRingNodes+nElementRing-1;
524 vertices[1] = row+(!fullTube)*nElementRing + nRingNodes;
525 vertices[2] = row + nElementRing -1;
526
527 factory.insertBoundarySegment(vertices,
528 SegmentPointer(new TubeUpperSegment<ctype>(localCenter, innerRadius, -stepLength,
529 fromAngle+(nElementRing-1)*stepRing*angleDiff, toAngle)));
530
531 vertices[0] += innerNodes;
532 vertices[1] += innerNodes;
533 vertices[2] += innerNodes;
534
535 factory.insertBoundarySegment(vertices,
536 SegmentPointer(new TubeUpperSegment<ctype>(localCenter, innerRadius+thickness, -stepLength,
537 fromAngle+(nElementRing-1)*stepRing*angleDiff, toAngle)));
538
539 }
540
541 // add lateral segments at the front and the back
542 std::vector<unsigned int> vert(4);
543
544 // offset to last back nodes
545 int lastRow = nElementLength*nRingNodes;
546 FieldVector<ctype,3> lastCenter = center;
547 lastCenter[0] += length;
548
549 for (size_t j=0; j<nElementRing-1; j++) {
550
551 // front side
552 vert[0] = j; vert[1] = j+1; vert[2] = j + innerNodes; vert[3] = j+innerNodes +1;
553
554 factory.insertBoundarySegment(vert,
555 SegmentPointer(new TubeLateralSegment<ctype>(center, innerRadius, thickness,
556 fromAngle+j*stepRing*angleDiff, fromAngle+(j+1)*stepRing*angleDiff)));
557
558 // back side
559 vert[0] += lastRow; vert[1] += lastRow; vert[2] += lastRow; vert[3] += lastRow;
560
561 factory.insertBoundarySegment(vert,
562 SegmentPointer(new TubeLateralSegment<ctype>(lastCenter, innerRadius, thickness,
563 fromAngle+j*stepRing*angleDiff, fromAngle+(j+1)*stepRing*angleDiff)));
564
565
566 }
567
568 //last segments differ depending on if the tube is closed
569
570 // front side
571 vert[0] = nElementRing-1; vert[1] = (!fullTube)*nElementRing;
572 vert[2] = vert[0] + innerNodes; vert[3] = vert[1] + innerNodes;
573
574 factory.insertBoundarySegment(vert,
575 SegmentPointer(new TubeLateralSegment<ctype>(center, innerRadius, thickness,
576 fromAngle+(nElementRing-1)*stepRing*angleDiff, toAngle)));
577
578 // back side
579 vert[0] += lastRow; vert[1] += lastRow; vert[2] += lastRow; vert[3] += lastRow;
580
581 factory.insertBoundarySegment(vert,
582 SegmentPointer(new TubeLateralSegment<ctype>(lastCenter, innerRadius, thickness,
583 fromAngle+(nElementRing-1)*stepRing*angleDiff, toAngle)));
584
585 }
586 // //////////////////////////////////////
587 // Finish initialization
588 // //////////////////////////////////////
589 std::unique_ptr<GridType> grid = std::unique_ptr<GridType>(factory.createGrid());
590
591 // create boundarypatches if desired
592 if (innerBoundary) {
593
594 BitSetVector<1> innerVert(2*(nElementLength+1)*nRingNodes);
595 size_t half = innerVert.size()/2;
596
597 for (size_t i=0; i<half; i++)
598 innerVert[i] = true;
599
600 innerBoundary->setup(grid->levelGridView(0),innerVert);
601 }
602
603 if (outerBoundary) {
604
605 BitSetVector<1> outerVert(2*(nElementLength+1)*nRingNodes);
606 size_t half = outerVert.size()/2;
607
608 for (size_t i=0; i<half; i++)
609 outerVert[half+i] = true;
610
611 outerBoundary->setup(grid->levelGridView(0),outerVert);
612 }
613
614 return grid;
615}
616
637template <class GridType>
638std::unique_ptr<GridType> makeTorus(
639 const Dune::FieldVector<typename GridType::ctype, 3>& center,
640 const std::array<Dune::FieldVector<typename GridType::ctype, 3>, 3> plane,
641 const typename GridType::ctype thickness,
642 const typename GridType::ctype horizontalRadius,
643 const typename GridType::ctype verticalRadius,
644 unsigned int nElementHorizontal, unsigned int nElementVertical,
645 BoundaryPatch<typename GridType::LevelGridView>* innerBoundary = NULL,
646 BoundaryPatch<typename GridType::LevelGridView> * outerBoundary = NULL)
647{
648 static_assert(GridType::dimension == 3,
649 "makeTorus can only be called for 3-dimensional grids");
650 using ctype = typename GridType::ctype;
651 using FV = Dune::FieldVector<ctype,3>;
652
653 Dune::GridFactory<GridType> factory;
654
655 // ///////////////////////
656 // Insert vertices
657 // ///////////////////////
658
659 // make circle around (0,0)
660 ArcOfCircle hArc(Dune::FieldVector<ctype,2>(0), horizontalRadius+verticalRadius+thickness, 0, 2*M_PI);
661 ArcOfCircle vInnerArc(Dune::FieldVector<ctype,2>(0), verticalRadius, 0, 2*M_PI);
662 ArcOfCircle vOuterArc(Dune::FieldVector<ctype,2>(0), verticalRadius+thickness, 0, 2*M_PI);
663
664 ctype stepH =1.0/nElementHorizontal;
665 ctype stepV = 1.0/nElementVertical;
666
667 // insert torus vertices
668 for (size_t j=0; j<nElementHorizontal; j++) {
669
670 // compute center of the vertical rings
671 FV ringCenter = center;
672 auto corr = hArc(j*stepH);
673 for (int i=0; i<2; i++)
674 ringCenter.axpy(corr[i],plane[i]);
675
676 //auto vPlane = ringCenter - plane[0];
677 //vPlane /= vPlane.two_norm();
678
679 auto vPlane = ringCenter - center;
680 vPlane /= vPlane.two_norm();
681
682 // compute coordinates on the inner torus
683 for (size_t i=0; i<nElementVertical; i++) {
684
685 auto local = vInnerArc(i*stepV);
686 auto coord = ringCenter;
687 coord.axpy(local[0],vPlane);
688 coord.axpy(local[1],plane[2]);
689
690 factory.insertVertex(coord);
691 }
692
693 // compute coordinates on the outer torus
694 for (size_t i=0; i<nElementVertical; i++) {
695
696 auto local = vOuterArc(i*stepV);
697 auto coord = ringCenter;
698 coord.axpy(local[0],vPlane);
699 coord.axpy(local[1],plane[2]);
700
701 factory.insertVertex(coord);
702 }
703 }
704
705 unsigned int innerNodes = nElementHorizontal*nElementVertical;
706
707 // /////////////////
708 // Insert elements
709 // /////////////////
710
711 // if ring is closed last nodes are the first nodes
712
713 std::vector<unsigned int> cornerIDs(4);
714 auto gt = Dune::GeometryTypes::tetrahedron;
715
716 for (unsigned int i=0; i<nElementHorizontal-1; i++) {
717
718 unsigned int row = i*2*nElementVertical;
719
720 for (unsigned int j=0; j<nElementVertical-1; j++) {
721
722
723 std::vector<unsigned int> hexa = {row+j,row+j+1,
724 row+j+2*nElementVertical,row+j+1+2*nElementVertical,
725 row+j+nElementVertical,row+j+1+nElementVertical,
726 row+j+3*nElementVertical,
727 row+j+1+3*nElementVertical};
728
729 //add 6 tetrahedra
730 cornerIDs={hexa[0],hexa[1],hexa[3],hexa[7]};
731 factory.insertElement(gt,cornerIDs);
732 cornerIDs={hexa[0],hexa[5],hexa[1],hexa[7]};
733 factory.insertElement(gt,cornerIDs);
734 cornerIDs={hexa[0],hexa[4],hexa[5],hexa[7]};
735 factory.insertElement(gt,cornerIDs);
736 cornerIDs={hexa[0],hexa[6],hexa[4],hexa[7]};
737 factory.insertElement(gt,cornerIDs);
738 cornerIDs={hexa[0],hexa[2],hexa[6],hexa[7]};
739 factory.insertElement(gt,cornerIDs);
740 cornerIDs={hexa[0],hexa[3],hexa[2],hexa[7]};
741 factory.insertElement(gt,cornerIDs);
742
743 }
744
745 // last element in each vertical circle couples with the first nodes
746 std::vector<unsigned int> hexa = {row+nElementVertical-1,row,
747 row+3*nElementVertical-1,row+2*nElementVertical,
748 row+2*nElementVertical-1,
749 row+nElementVertical,
750 row+4*nElementVertical-1,
751 row+3*nElementVertical};
752 //add 6 tetrahedra
753 cornerIDs={hexa[0],hexa[1],hexa[3],hexa[7]};
754 factory.insertElement(gt,cornerIDs);
755 cornerIDs={hexa[0],hexa[5],hexa[1],hexa[7]};
756 factory.insertElement(gt,cornerIDs);
757 cornerIDs={hexa[0],hexa[4],hexa[5],hexa[7]};
758 factory.insertElement(gt,cornerIDs);
759 cornerIDs={hexa[0],hexa[6],hexa[4],hexa[7]};
760 factory.insertElement(gt,cornerIDs);
761 cornerIDs={hexa[0],hexa[2],hexa[6],hexa[7]};
762 factory.insertElement(gt,cornerIDs);
763 cornerIDs={hexa[0],hexa[3],hexa[2],hexa[7]};
764 factory.insertElement(gt,cornerIDs);
765 }
766
767 // the last row couples with the first
768 unsigned int row = (nElementHorizontal-1)*2*nElementVertical;
769
770 for (unsigned int j=0; j<nElementVertical-1; j++) {
771
772
773 std::vector<unsigned int> hexa = {row+j,row+j+1, j,j+1,
774 row+j+nElementVertical,row+j+1+nElementVertical,
775 j+nElementVertical,
776 j+1+nElementVertical};
777
778 //add 6 tetrahedra
779 cornerIDs={hexa[0],hexa[1],hexa[3],hexa[7]};
780 factory.insertElement(gt,cornerIDs);
781 cornerIDs={hexa[0],hexa[5],hexa[1],hexa[7]};
782 factory.insertElement(gt,cornerIDs);
783 cornerIDs={hexa[0],hexa[4],hexa[5],hexa[7]};
784 factory.insertElement(gt,cornerIDs);
785 cornerIDs={hexa[0],hexa[6],hexa[4],hexa[7]};
786 factory.insertElement(gt,cornerIDs);
787 cornerIDs={hexa[0],hexa[2],hexa[6],hexa[7]};
788 factory.insertElement(gt,cornerIDs);
789 cornerIDs={hexa[0],hexa[3],hexa[2],hexa[7]};
790 factory.insertElement(gt,cornerIDs);
791
792 }
793
794 // last element in each vertical circle couples with the first nodes
795 std::vector<unsigned int> hexa = {row+nElementVertical-1,row,
796 nElementVertical-1,0,
797 row+2*nElementVertical-1,
798 row+nElementVertical,
799 2*nElementVertical-1,
800 nElementVertical};
801 //add 6 tetrahedra
802 cornerIDs={hexa[0],hexa[1],hexa[3],hexa[7]};
803 factory.insertElement(gt,cornerIDs);
804 cornerIDs={hexa[0],hexa[5],hexa[1],hexa[7]};
805 factory.insertElement(gt,cornerIDs);
806 cornerIDs={hexa[0],hexa[4],hexa[5],hexa[7]};
807 factory.insertElement(gt,cornerIDs);
808 cornerIDs={hexa[0],hexa[6],hexa[4],hexa[7]};
809 factory.insertElement(gt,cornerIDs);
810 cornerIDs={hexa[0],hexa[2],hexa[6],hexa[7]};
811 factory.insertElement(gt,cornerIDs);
812 cornerIDs={hexa[0],hexa[3],hexa[2],hexa[7]};
813 factory.insertElement(gt,cornerIDs);
814
815 // //////////////////////////////////////
816 // Finish initialization
817 // //////////////////////////////////////
818 std::unique_ptr<GridType> grid = std::unique_ptr<GridType>(factory.createGrid());
819
820 // create boundarypatches if desired
821 if (innerBoundary) {
822
823 Dune::BitSetVector<1> innerVert(2*innerNodes,false);
824
825 for (unsigned i=0; i<nElementHorizontal; i++) {
826 auto row = i*2*nElementVertical;
827 for (unsigned j=0; j<nElementVertical; j++)
828 innerVert[row+j] = true;
829 }
830
831 innerBoundary->setup(grid->levelGridView(0),innerVert);
832 }
833
834 if (outerBoundary) {
835
836 Dune::BitSetVector<1> outerVert(2*innerNodes,false);
837
838 for (unsigned i=0; i<nElementHorizontal; i++) {
839 auto row = i*2*nElementVertical;
840 for (unsigned j=nElementVertical; j<2*nElementVertical; j++)
841 outerVert[row+j] = true;
842 }
843
844 outerBoundary->setup(grid->levelGridView(0),outerVert);
845 }
846
847 return grid;
848}
849
850#endif
makeTubeSegment3D
std::unique_ptr< GridType > makeTubeSegment3D(const Dune::FieldVector< typename GridType::ctype, 3 > &center, const typename GridType::ctype thickness, const typename GridType::ctype length, const typename GridType::ctype innerRadius, const typename GridType::ctype fromAngle, const typename GridType::ctype toAngle, unsigned int nElementRing, unsigned int nElementLength, bool closeTube, unsigned int axis=0, bool tetrahedra=true, bool parametrizedBoundary=true, BoundaryPatch< typename GridType::LevelGridView > *innerBoundary=NULL, BoundaryPatch< typename GridType::LevelGridView > *outerBoundary=NULL)
Create a 3D tube segment prism grid along a prescribed Euclidean axis.
Definition makering.hh:268
makeTorus
std::unique_ptr< GridType > makeTorus(const Dune::FieldVector< typename GridType::ctype, 3 > &center, const std::array< Dune::FieldVector< typename GridType::ctype, 3 >, 3 > plane, const typename GridType::ctype thickness, const typename GridType::ctype horizontalRadius, const typename GridType::ctype verticalRadius, unsigned int nElementHorizontal, unsigned int nElementVertical, BoundaryPatch< typename GridType::LevelGridView > *innerBoundary=NULL, BoundaryPatch< typename GridType::LevelGridView > *outerBoundary=NULL)
Create a 3D tube segment prism grid along the x-axis.
Definition makering.hh:638
makeRingSegment2D
void makeRingSegment2D(const Dune::FieldVector< typename GridType::ctype, 2 > &center, const typename GridType::ctype thickness, const typename GridType::ctype innerRadius, const typename GridType::ctype fromAngle, const typename GridType::ctype toAngle, GridType &grid, size_t nElement, bool closeRing, bool quadrilateral=false, BoundaryPatch< typename GridType::LevelGridView > *innerBoundary=NULL, BoundaryPatch< typename GridType::LevelGridView > *outerBoundary=NULL)
Definition makering.hh:20
DUNE_THROW
#define DUNE_THROW(E,...)
local
LocalIndex & local()
Dune::BitSetVector::size
size_type size() const
Dune::FieldVector::two_norm
constexpr FieldTraits< value_type >::real_type two_norm() const
Dune::Grid::levelGridView
LevelGridView levelGridView(int level) const
Dune::GridFactory::insertElement
virtual void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices)
Dune::GridFactory::insertVertex
virtual void insertVertex(const FieldVector< ctype, dimworld > &pos)
Dune::GridFactory::insertBoundarySegment
virtual void insertBoundarySegment(const std::vector< unsigned int > &vertices)
Dune::GridFactory::createGrid
virtual std::unique_ptr< GridType > createGrid()
ArcOfCircle
Definition arcofcircle.hh:11
TubeLowerSegment
Boundary segments for a 3D-tube/tube-segment prism grid that is aligned with the x-axis.
Definition arcofcircle.hh:50
TubeUpperSegment
Boundary segments for a 3D-tube/tube-segment prism grid that is aligned with the x-axis.
Definition arcofcircle.hh:91
TubeLateralSegment
Lateral boundary segments for a 3D-tube/tube-segment prism grid that is aligned with the x-axis.
Definition arcofcircle.hh:132
BoundaryPatch
Encapsulate a part of a grid boundary.
Definition boundarypatch.hh:218
std::forward
T forward(T... args)