KASKADE 7 development version
gridGeneration.hh
Go to the documentation of this file.
1/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2/* */
3/* This file is part of the library KASKADE 7 */
4/* https://www.zib.de/research/projects/kaskade7-finite-element-toolbox */
5/* */
6/* Copyright (C) 2013-2023 Zuse Institute Berlin */
7/* */
8/* KASKADE 7 is distributed under the terms of the ZIB Academic License. */
9/* see $KASKADE/academic.txt */
10/* */
11/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
12/*
13 * Created on: 25.06.2013
14 * Author: bzflubko
15 */
16
17#ifndef KASKADE_GRID_GENERATION_HH_
18#define KASKADE_GRID_GENERATION_HH_
19
20#include <boost/math/constants/constants.hpp>
21
22#include <algorithm>
23#include <cmath>
24#include <memory>
25#include <type_traits>
26#include <utility>
27#include <vector>
28
29#include <dune/common/fvector.hh>
30// #include "dune/grid/uggrid.hh"
31
32#include "fem/fixdune.hh"
33#include "GridGeneration/cube.hh"
34#include "utilities/detailed_exception.hh"
35#include "utilities/timing.hh"
36
37namespace Kaskade
38{
42 namespace GridGeneration_Detail
43 {
55 template <class Type,
56 int dim = std::is_same_v<Type,Cube<typename Type::Real>>? 3: 2>
57 std::pair<std::vector<Dune::FieldVector<typename Type::Real,dim>>, std::vector<std::vector<unsigned int>>> extractSimplices(std::vector<Type> const& cubes)
58 {
59 using Vec = Dune::FieldVector<typename Type::Real,dim>;
60
61 // A lexical ordering on the vertices (with tolerance).
62 typename Type::Real tol = 1e-6;
63 auto vLess = [&](Vec const& x, Vec const& y)
64 {
65 for (int i=0; i<dim; ++i)
66 if (x[i] < y[i]-tol)
67 return true;
68 else if (x[i] > y[i]+tol)
69 return false;
70
71 return false; // equal, i.e. all components within tolerance, and thus not less
72 };
73
74 auto vEq = [&](Vec const& x, Vec const& y)
75 {
76 return !(vLess(x,y) || vLess(y,x));
77 };
78
79 // First we create a list of unique vertices, eliminating any duplicates,
80 // such that we get unique and stable vertex ids.
81 std::vector<Vec> vertices;
82 for (auto const& cube: cubes)
83 {
84 auto const& localVertices = cube.getVertices();
85 vertices.insert(vertices.end(),localVertices.begin(),localVertices.end());
86 }
87 std::sort(vertices.begin(),vertices.end(),vLess);
88 vertices.erase(std::unique(vertices.begin(),vertices.end(),vEq),vertices.end());
89
90
91
92 std::vector<std::vector<unsigned int>> tetrahedra; // TODO: use std::array here, the number of corners is fixed
93
94 for (Type const& cube: cubes)
95 {
96 auto const& localVertices = cube.getVertices();
97 auto const& localTetrahedra = cube.getSimplices();
98 size_t const nVertices = localVertices.size();
99
100 // maps local vertex ids to global ones (i.e. the index in 'vertices')
101 std::vector<size_t> vertexIdMap(nVertices);
102 for(size_t i=0; i<nVertices; ++i)
103 {
104 auto it = std::lower_bound(vertices.begin(),vertices.end(),localVertices[i],vLess);
105 assert(it != vertices.end());
106 vertexIdMap[i] = it - vertices.begin();
107 }
108
109 // map local connectivity onto global one (wrt. numbering in vector 'vertices')
110 std::vector<unsigned int> newTet;
111 for(std::vector<unsigned int> const& tet: localTetrahedra)
112 {
113 newTet.resize(tet.size());
114 for(size_t i=0; i<newTet.size(); ++i)
115 newTet[i] = vertexIdMap[tet[i]];
116 tetrahedra.push_back(newTet);
117 }
118 }
119
120 return std::make_pair(vertices,tetrahedra);
121 }
122
123
132 template <class Scalar>
133 std::vector<Square<Scalar> > fillRectangle(Dune::FieldVector<Scalar,2> const& x0, Dune::FieldVector<Scalar,2> const& dx, Dune::FieldVector<Scalar,2> const& dh, bool symmetric=true)
134 {
135 std::vector<Square<Scalar> > squares;
136 Dune::FieldVector<Scalar,2> dh2(dh);
137 for(int i=0; i<2; ++i)
138 if(dh[i]>dx[i]) // if dh larger than sidelength, set dh = dx
139 dh2[i] = dx[i];
140
141
142 size_t xi_end = (size_t)round(fabs(dx[0]/dh2[0])),
143 yi_end = (size_t)round(fabs(dx[1]/dh2[1]));
144 Dune::FieldVector<Scalar,2> dsquare(dh2), x(x0);
145
146 for(size_t xi=0; xi<xi_end; ++xi)
147 {
148 x[0] = x0[0] + xi*dh2[0];
149 for(size_t yi=0; yi<yi_end; ++yi)
150 {
151 x[1] = x0[1] + yi*dh2[1];
152 squares.push_back(Square<Scalar>(x,dsquare,symmetric));
153 }
154 }
155
156 return squares;
157 }
158
167 template <class Scalar>
168 std::vector<Square<Scalar> > fillRectangle(Dune::FieldVector<Scalar,2> const& x0, Dune::FieldVector<Scalar,2> const& dx, Scalar dh, bool symmetric=true)
169 {
170 return fillRectangle(x0,dx,Dune::FieldVector<Scalar,2>(dh),symmetric);
171 }
172
181 template <class Scalar>
182 std::vector<Cube<Scalar>> fillCuboid(Dune::FieldVector<Scalar,3> const& x0, Dune::FieldVector<Scalar,3> const& dx,
183 Dune::FieldVector<Scalar,3> const& dh, bool symmetric=true)
184 {
185 std::vector<Cube<Scalar>> cubes;
186
187 Dune::FieldVector<Scalar,3> dh2(dh);
188 for(int i=0; i<3; ++i)
189 if(dh[i]>dx[i]) // if dh larger than sidelength, set dh = dx
190 dh2[i] = dx[i];
191
192 size_t xi_end = (size_t)round(fabs(dx[0]/dh2[0])),
193 yi_end = (size_t)round(fabs(dx[1]/dh2[1])),
194 zi_end = (size_t)round(fabs(dx[2]/dh2[2]));
195
196 if (xi_end==0 || yi_end==0 || zi_end==0)
197 throw GridException("empty cubes list",__FILE__,__LINE__);
198
199 Dune::FieldVector<Scalar,3> x(x0);
200
201 for(size_t xi=0; xi<xi_end; ++xi)
202 {
203 x[0] = x0[0] + xi*dh2[0];
204 for(size_t yi=0; yi<yi_end; ++yi)
205 {
206 x[1] = x0[1] + yi*dh2[1];
207 for(size_t zi=0; zi<zi_end; ++zi)
208 {
209 x[2] = x0[2] + zi*dh2[2];
210 cubes.push_back(Cube<Scalar>(x,dh2,symmetric));
211// std::cout << "insert cube at " << x << " width " << dh2 << "\n";
212 }
213 }
214 }
215
216 return cubes;
217 }
218
219
228 template <class Scalar>
229 std::vector<Cube<Scalar>> fillCuboid(Dune::FieldVector<Scalar,3> const& x0, Dune::FieldVector<Scalar,3> const& dx,
230 Scalar dh, bool symmetric=true)
231 {
232 return fillCuboid(x0,dx,Dune::FieldVector<Scalar,3>(dh),symmetric);
233 }
234
235
236
237 template<class Scalar, int dim>
238 std::vector<Cube<Scalar>> fillTwoLayerCuboid(Dune::FieldVector<Scalar,dim> x0, Dune::FieldVector<Scalar,dim> const& dx1, Dune::FieldVector<Scalar,dim> const& dx2,
239 Dune::FieldVector<Scalar,dim> const& dh1, Scalar dh2, bool symmetric=true)
240 {
241 auto cubes = fillCuboid(x0,dx1,dh1,symmetric);
242
243 x0[dim-1] += dx1[dim-1];
244 auto cubes2 = fillCuboid(x0,dx2,dh2,symmetric);
245 for( Cube<Scalar> const& cube : cubes2 ) cubes.push_back(cube);
246
247 return cubes;
248 }
249
250 template<class Scalar, int dim>
251 std::vector<Cube<Scalar> > fillTwoLayerCuboid(Dune::FieldVector<Scalar,dim> x0, Dune::FieldVector<Scalar,dim> const& dx1, Dune::FieldVector<Scalar,dim> const& dx2,
252 Dune::FieldVector<Scalar,dim> const& dh1, Dune::FieldVector<Scalar,dim> const& dh2, bool symmetric=true)
253 {
254 auto cubes = fillCuboid(x0,dx1,dh1,symmetric);
255
256 x0[dim-1] += dx1[dim-1];
257 auto cubes2 = fillCuboid(x0,dx2,dh2,symmetric);
258 for( Cube<Scalar> const& cube : cubes2 ) cubes.push_back(cube);
259
260 return cubes;
261 }
262
263 template <int dim>
264 struct FillCuboid;
265
266 template <>
267 struct FillCuboid<3>
268 {
269 template <class Scalar, class EdgeLength>
270 static std::vector<Cube<Scalar> > apply(Dune::FieldVector<Scalar,3> const& x0, Dune::FieldVector<Scalar,3> const& dx,
271 EdgeLength const& dh, bool symmetric=true)
272 {
273 return fillCuboid(x0, dx, dh, symmetric);
274 }
275 };
276
277 template <>
278 struct FillCuboid<2>
279 {
280 template <class Scalar, class EdgeLength>
281 static std::vector<Square<Scalar> > apply(Dune::FieldVector<Scalar,2> const& x0, Dune::FieldVector<Scalar,2> const& dx, EdgeLength const& dh, bool symmetric=true)
282 {
283 return fillRectangle(x0, dx, dh, symmetric);
284 }
285 };
286
287
298 template <class Scalar, int dim>
299 std::vector< typename std::conditional<dim==3,Cube<Scalar>,Square<Scalar> >::type >
300 fillLShape(Dune::FieldVector<Scalar,dim> x0, Scalar dx, Scalar dy, Dune::FieldVector<Scalar,dim> const& thickness, Scalar dh, bool symmetric=true)
301 {
302 // fill horizontal line
303 Dune::FieldVector<Scalar,dim> d0(thickness); d0[0] = dx;
304 auto cubes = FillCuboid<dim>::apply(x0,d0,dh,symmetric);
305
306 // fill vertical line
307 d0 = thickness; d0[1] = dy;
308 x0[1] += thickness[1];
309 auto cubes2 = FillCuboid<dim>::apply(x0,d0,dh,symmetric);
310
311 // merge
312 for(auto const& cube : cubes2)
313 cubes.push_back(cube);
314
315 return cubes;
316 }
317
328 template <class Scalar, int dim>
329 std::vector< typename std::conditional<dim==3,Cube<Scalar>,Square<Scalar>>::type>
330 fillTShape(Dune::FieldVector<Scalar,dim> x0, Scalar dx, Scalar dy, Dune::FieldVector<Scalar,dim> const& thickness, Scalar dh, bool symmetric=true)
331 {
332 // fill vertical line
333 x0[0] -= 0.5*thickness[0];
334 Dune::FieldVector<Scalar,dim> d0(thickness); d0[1] = dy - thickness[1];
335 auto cubes = FillCuboid<dim>::apply(x0,d0,dh,symmetric);
336
337 // fill vertical line
338 x0[0] += 0.5*thickness[0] - 0.5*dx;
339 x0[1] += d0[1];
340 d0[1] = thickness[1]; d0[0] = dx;
341 auto cubes2 = FillCuboid<dim>::apply(x0,d0,dh,symmetric);
342
343 // merge
344 for(auto const& cube : cubes2) cubes.push_back(cube);
345
346 return cubes;
347 }
348
349
350 } // end of namespace GridGeneration_Detail
370 template <class Grid, class Type>
371 std::unique_ptr<Grid> createGrid(std::vector<Type> const& cubes)
372 {
373 Timings& timer = Timings::instance();
374 timer.start("extract simplices");
375 auto simplices = GridGeneration_Detail::extractSimplices(cubes);
376 timer.stop("extract simplices");
377 Dune::GridFactory<Grid> factory;
378 Dune::GeometryType gt(Dune::GeometryType::simplex,Type::dim);
379
380 timer.start("insert vertex");
381 for (Dune::FieldVector<double,Type::dim> const& vertex: simplices.first)
382 factory.insertVertex(vertex);
383 timer.stop("insert vertex");
384 timer.start("insert element");
385 for (std::vector<unsigned int> const& elem: simplices.second)
386 factory.insertElement(gt,elem);
387 timer.stop("insert element");
388
389 return std::unique_ptr<Grid>(factory.createGrid());
390 }
391
413 template <class Grid, class Scalar, int dim, class EdgeLength>
414 std::unique_ptr<Grid> createCuboid(Dune::FieldVector<Scalar,dim> const& x0, Dune::FieldVector<Scalar,dim> const& dx,
415 EdgeLength const& dh, bool symmetric=true)
416 {
417 Timings& timer = Timings::instance();
418 timer.start("create cubes");
419 auto const& cubes = GridGeneration_Detail::FillCuboid<dim>::apply(x0,dx,dh,symmetric);
420 timer.stop("create cubes");
421 return createGrid<Grid>(cubes);
422 }
423
430 template <class Grid>
431 std::unique_ptr<Grid> createUnitSquare(double dh = 1., bool symmetric=true)
432 {
433 using Vec = Dune::FieldVector<double,2>;
434 return createCuboid<Grid>(Vec(0.), Vec(1.), dh, symmetric);
435 }
436
442 template <class Grid>
443 std::unique_ptr<Grid> createUnitCube(double dh = 1., bool symmetric=true)
444 {
445 typedef Dune::FieldVector<double,3> Vec;
446 return createCuboid<Grid>(Vec(0.), Vec(1.), dh, symmetric);
447 }
448
449
459 template <class Grid, class Scalar>
460 std::unique_ptr<Grid> createRectangle(Dune::FieldVector<Scalar,2> const& x0, Dune::FieldVector<Scalar,2> const& dx, Scalar dh, bool symmetric=true)
461 {
462 return createGrid<Grid>(GridGeneration_Detail::fillRectangle(x0,dx,dh,symmetric));
463 }
464
476 template <class Grid, class Scalar, int dim>
477 std::unique_ptr<Grid> createLShape(Dune::FieldVector<Scalar,dim> const& x0, Scalar dx, Scalar dy,
478 Dune::FieldVector<Scalar,dim> const& thickness, Scalar dh,
479 bool symmetric=true)
480 {
481 return createGrid<Grid>(GridGeneration_Detail::fillLShape(x0,dx,dy,thickness,dh,symmetric));
482 }
483
495 template <class Grid, class Scalar, int dim>
496 std::unique_ptr<Grid> createTShape(Dune::FieldVector<Scalar,dim> const& x0, Scalar dx, Scalar dy, Dune::FieldVector<Scalar,dim> const& thickness, Scalar dh, bool symmetric=true)
497 {
498 return createGrid<Grid>(GridGeneration_Detail::fillTShape(x0,dx,dy,thickness,dh,symmetric));
499 }
500
504 template <class Grid, class Scalar, int dim>
505 std::unique_ptr<Grid> createTwoLayerCuboid(Dune::FieldVector<Scalar,dim> const& x0, Dune::FieldVector<Scalar,dim> const& dx1, Dune::FieldVector<Scalar,dim> const& dx2,
506 Dune::FieldVector<Scalar,dim> const& dh1, Scalar dh2, bool symmetric=true)
507 {
508 return createGrid<Grid>(GridGeneration_Detail::fillTwoLayerCuboid(x0,dx1,dx2,dh1,dh2,symmetric));
509 }
510
514 template <class Grid, class Scalar, int dim>
515 std::unique_ptr<Grid> createTwoLayerCuboid(Dune::FieldVector<Scalar,dim> const& x0, Dune::FieldVector<Scalar,dim> const& dx1, Dune::FieldVector<Scalar,dim> const& dx2,
516 Dune::FieldVector<Scalar,dim> const& dh1, Dune::FieldVector<Scalar,dim> const& dh2, bool symmetric=true)
517 {
518 return createGrid<Grid>(GridGeneration_Detail::fillTwoLayerCuboid(x0,dx1,dx2,dh1,dh2,symmetric));
519 }
520
521 // --------------------------------------------------------------------------------------------------------------------
522 // --------------------------------------------------------------------------------------------------------------------
523
538 template <class Grid>
539 std::unique_ptr<Grid> createHexahedron()
540 {
541 Dune::GridFactory<Grid> factory;
542 Dune::GeometryType gt(Dune::GeometryType::simplex,3);
543
544 for (int i=0; i<2; ++i)
545 for (int j=0; j<2; ++j)
546 for (int k=0; k<2; ++k)
547 {
548 Dune::FieldVector<double,3> x;
549 x[0] = i;
550 x[1] = j;
551 x[2] = k;
552 factory.insertVertex(x);
553 }
554
555 factory.insertElement(gt,std::vector<unsigned int>{7,1,4,2});
556 factory.insertElement(gt,std::vector<unsigned int>{0,4,1,2});
557 factory.insertElement(gt,std::vector<unsigned int>{5,1,4,7});
558 factory.insertElement(gt,std::vector<unsigned int>{6,2,4,7});
559 factory.insertElement(gt,std::vector<unsigned int>{3,1,7,2});
560
561 return std::unique_ptr<Grid>(factory.createGrid());
562 }
563
573 template <class Grid>
574 std::unique_ptr<Grid> createOctahedron(Dune::FieldVector<typename Grid::ctype,3> extent = Dune::FieldVector<typename Grid::ctype,3>(1.0))
575 {
576 Dune::GridFactory<Grid> factory;
577 Dune::GeometryType gt(Dune::GeometryType::simplex,3);
578
579 for (int i=0; i<3; ++i)
580 {
581 Dune::FieldVector<double,3> x;
582 x[i] = extent[i];
583 factory.insertVertex(x);
584 factory.insertVertex(-x);
585 }
586 factory.insertVertex(Dune::FieldVector<double,3>());
587
588 factory.insertElement(gt,std::vector<unsigned int>{0,2,4,6});
589 factory.insertElement(gt,std::vector<unsigned int>{0,3,4,6});
590 factory.insertElement(gt,std::vector<unsigned int>{0,2,5,6});
591 factory.insertElement(gt,std::vector<unsigned int>{0,3,5,6});
592 factory.insertElement(gt,std::vector<unsigned int>{1,2,4,6});
593 factory.insertElement(gt,std::vector<unsigned int>{1,3,4,6});
594 factory.insertElement(gt,std::vector<unsigned int>{1,2,5,6});
595 factory.insertElement(gt,std::vector<unsigned int>{1,3,5,6});
596
597 return std::unique_ptr<Grid>(factory.createGrid());
598 }
599
608 template <class Grid>
609 std::unique_ptr<Grid> createIcosahedron()
610 {
611 Dune::GridFactory<Grid> factory;
612 Dune::GeometryType gt(Dune::GeometryType::simplex,3);
613
614 auto vector = [=](double x, double y, double z)
615 {
616 Dune::FieldVector<double,3> p;
617 p[0] = x; p[1] = y; p[2] = z;
618 return p;
619 };
620
621 double s = (1+std::sqrt(5))/2;
622
623 factory.insertVertex(vector( 0, 1, s)); // 0
624 factory.insertVertex(vector( 0, 1,-s)); // 1
625 factory.insertVertex(vector( 0,-1, s)); // 2
626 factory.insertVertex(vector( 0,-1,-s)); // 3
627 factory.insertVertex(vector( 1, s, 0)); // 4
628 factory.insertVertex(vector( 1,-s, 0)); // 5
629 factory.insertVertex(vector(-1, s, 0)); // 6
630 factory.insertVertex(vector(-1,-s, 0)); // 7
631 factory.insertVertex(vector( s, 0, 1)); // 8
632 factory.insertVertex(vector( s, 0,-1)); // 9
633 factory.insertVertex(vector(-s, 0, 1)); // 10
634 factory.insertVertex(vector(-s, 0,-1)); // 11
635 factory.insertVertex(vector( 0, 0, 0)); // 12
636
637 factory.insertElement(gt,std::vector<unsigned int>{0,2, 8,12});
638 factory.insertElement(gt,std::vector<unsigned int>{0,2,10,12});
639 factory.insertElement(gt,std::vector<unsigned int>{1,3, 9,12});
640 factory.insertElement(gt,std::vector<unsigned int>{1,3,11,12});
641
642 factory.insertElement(gt,std::vector<unsigned int>{4,6, 0,12});
643 factory.insertElement(gt,std::vector<unsigned int>{4,6, 1,12});
644 factory.insertElement(gt,std::vector<unsigned int>{5,7, 2,12});
645 factory.insertElement(gt,std::vector<unsigned int>{5,7, 3,12});
646
647 factory.insertElement(gt,std::vector<unsigned int>{8,9, 4,12});
648 factory.insertElement(gt,std::vector<unsigned int>{8,9, 5,12});
649 factory.insertElement(gt,std::vector<unsigned int>{10,11,6,12});
650 factory.insertElement(gt,std::vector<unsigned int>{10,11,7,12});
651
652 factory.insertElement(gt,std::vector<unsigned int>{0,4, 8,12});
653 factory.insertElement(gt,std::vector<unsigned int>{0,6,10,12});
654 factory.insertElement(gt,std::vector<unsigned int>{3,7,11,12});
655 factory.insertElement(gt,std::vector<unsigned int>{3,5,9,12});
656 factory.insertElement(gt,std::vector<unsigned int>{1,4,9,12});
657 factory.insertElement(gt,std::vector<unsigned int>{1,6,11,12});
658 factory.insertElement(gt,std::vector<unsigned int>{2,7,10,12});
659 factory.insertElement(gt,std::vector<unsigned int>{2,5,8,12});
660
661 return std::unique_ptr<Grid>(factory.createGrid());
662 }
663
664 //-----------------------------------------------------------------------------------------------
665
682 template <class Grid>
683 std::unique_ptr<Grid> createCylinder(double radius, double height, int layers, int circumEdges=6)
684 {
685 return createCylinder<Grid>(radius,height,layers,circumEdges,[](auto const& x) { return x; });
686 }
687
704 template <class Grid, class Deformation>
705 std::unique_ptr<Grid> createCylinder(double radius, double height, int layers, int circumEdges,
706 Deformation const& deformation)
707 {
708 using boost::math::double_constants::pi;
709
710 Dune::GridFactory<Grid> factory;
711 Dune::GeometryType gt(Dune::GeometryType::simplex,3);
712
713 auto vector = [](double x, double y, double z)
714 {
715 Dune::FieldVector<double,3> p;
716 p[0] = x; p[1] = y; p[2] = z;
717 return p;
718 };
719
720 double h = height/layers;
721 // Create vertices. These are circumEdges+1 on each of the layers+1 circular disks
722 for (int i=0; i<=layers; ++i)
723 {
724 double z = i*h;
725 factory.insertVertex(deformation(vector(0,0,z)));
726 for (int j=0; j<circumEdges; ++j)
727 {
728// double angle = 2*pi*(j-i/6.0)/circumEdges;
729 double angle = 2*pi*j/circumEdges;
730 auto v = vector(radius*std::cos(angle),radius*std::sin(angle),z);
731 factory.insertVertex(deformation(v));
732 }
733 }
734
735 // Create the cells. Each prismatic disc consists of circumEdges prisms
736 // connecting the central axis with a surface quadrilateral. These prisms are split
737 // into three tetrahedra.
738 for (int i=0; i<layers; ++i)
739 {
740 // The offset of the vertex indices: circular disks below current anti-prismatic disk
741 // times the number of vertices in each circular disk.
742 int n = circumEdges+1;
743 int idxOffset = i*n;
744
745 unsigned int c0 = idxOffset; // lower center vertex
746 unsigned int c1 = idxOffset+n; // top center vertex
747
748 for (int j=0; j<circumEdges; ++j)
749 {
750 unsigned int a0 = c0+1+j; // bottom left corner of boundary quadrilateral
751 unsigned int b0 = c0+1+(j+1)%circumEdges; // bottom right corner
752 unsigned int a1 = c1+1+j; // top left corner
753 unsigned int b1 = c1+1+(j+1)%circumEdges; // top right corner
754 factory.insertElement(gt,std::vector<unsigned int>{a0,b0,c0,c1});
755 factory.insertElement(gt,std::vector<unsigned int>{a0,b1,a1,c1});
756 factory.insertElement(gt,std::vector<unsigned int>{a0,b0,b1,c1});
757 }
758 }
759
760 return std::unique_ptr<Grid>(factory.createGrid());
761 }
762
763 //-----------------------------------------------------------------------------------------------
764
774 template <class Grid>
775 std::unique_ptr<Grid> createPentagon(double radius=1)
776 {
777
778 Dune::GridFactory<Grid> factory;
779 Dune::GeometryType gt(Dune::GeometryType::simplex,2);
780
781 auto vector = [=](double x, double y)
782 {
783 Dune::FieldVector<double,2> p;
784 p[0] = x; p[1] = y;
785 return p;
786 };
787
788 double c1, c2, s1, s2;
789 c1 = 0.25*(std::sqrt(5)-1);
790 c2 = 0.25*(std::sqrt(5)+1);
791 s1 = 0.25*(std::sqrt(10+2*std::sqrt(5)));
792 s2 = 0.25*(std::sqrt(10-2*std::sqrt(5)));
793
794 factory.insertVertex(radius*vector( 1., 0.)); // 0
795 factory.insertVertex(radius*vector( c1, s1)); // 1
796 factory.insertVertex(radius*vector(-1.*c2,s2)); // 2
797 factory.insertVertex(radius*vector(-1.*c2,-1.*s2)); // 3
798 factory.insertVertex(radius*vector( c1, -1.0*s1)); // 4
799 factory.insertVertex(radius*vector( 0., 0.)); // 5
800
801
802 factory.insertElement(gt,std::vector<unsigned int>{0,1,5});
803 factory.insertElement(gt,std::vector<unsigned int>{1,2,5});
804 factory.insertElement(gt,std::vector<unsigned int>{2,3,5});
805 factory.insertElement(gt,std::vector<unsigned int>{3,4,5});
806 factory.insertElement(gt,std::vector<unsigned int>{4,0,5});
807
808 return std::unique_ptr<Grid>(factory.createGrid());
809 }
810
811
829 template <class Grid>
830 std::unique_ptr<Grid> createUshape(double l, double d, double eps)
831 {
832 assert(eps > 0);
833 assert(d > 0);
834 assert(l > d);
835
836 Dune::GridFactory<Grid> factory;
837 Dune::GeometryType gt(Dune::GeometryType::simplex,2);
838
839 auto vector = [=](double x, double y)
840 {
841 Dune::FieldVector<double,2> p;
842 p[0] = x; p[1] = y;
843 return p;
844 };
845
846 factory.insertVertex(vector( 0, 0)); // 0
847 factory.insertVertex(vector( (l+d)/2, 0)); // 1
848 factory.insertVertex(vector( l, 0)); // 2
849 factory.insertVertex(vector( l, d)); // 3
850 factory.insertVertex(vector( (l+d)/2, d)); // 4
851 factory.insertVertex(vector( d, d)); // 5
852 factory.insertVertex(vector( d, d+eps)); // 6
853 factory.insertVertex(vector( (l+d)/2, d+eps)); // 7
854 factory.insertVertex(vector( l, d+eps)); // 8
855 factory.insertVertex(vector( l, 2*d+eps)); // 9
856 factory.insertVertex(vector( (l+d)/2, 2*d+eps)); // 10
857 factory.insertVertex(vector( 0, 2*d+eps)); // 11
858 factory.insertVertex(vector( 0, d+eps/2)); // 12
859
860 factory.insertElement(gt,std::vector<unsigned int>{0,1,5});
861 factory.insertElement(gt,std::vector<unsigned int>{1,4,5});
862 factory.insertElement(gt,std::vector<unsigned int>{1,2,4});
863 factory.insertElement(gt,std::vector<unsigned int>{2,3,4});
864 factory.insertElement(gt,std::vector<unsigned int>{5,6,12});
865 factory.insertElement(gt,std::vector<unsigned int>{6,7,10});
866 factory.insertElement(gt,std::vector<unsigned int>{7,8,9});
867 factory.insertElement(gt,std::vector<unsigned int>{7,9,10});
868 factory.insertElement(gt,std::vector<unsigned int>{6,10,11});
869 factory.insertElement(gt,std::vector<unsigned int>{6,11,12});
870 factory.insertElement(gt,std::vector<unsigned int>{0,5,12});
871
872 return std::unique_ptr<Grid>(factory.createGrid());
873 }
874
875} // end of namespace Kaskade
876
877#endif /* GRIDGENERATION_HH_ */
Kaskade::Timings
Supports gathering and reporting execution times information for nested program parts.
Definition: timing.hh:64
Kaskade::Timings::instance
static Timings & instance()
Returns a reference to a single default instance.
Kaskade::Timings::stop
void stop(std::string const &name)
Stops the timing of given section.
Kaskade::Timings::start
struct Times const * start(std::string const &name)
Starts or continues the timing of given section.
detailed_exception.hh
fixdune.hh
This file contains various utility functions that augment the basic functionality of Dune.
Kaskade::createLShape
std::unique_ptr< Grid > createLShape(Dune::FieldVector< Scalar, dim > const &x0, Scalar dx, Scalar dy, Dune::FieldVector< Scalar, dim > const &thickness, Scalar dh, bool symmetric=true)
Fill L-shaped domain with cubes.
Definition: gridGeneration.hh:477
Kaskade::createUnitSquare
std::unique_ptr< Grid > createUnitSquare(double dh=1., bool symmetric=true)
Creates a unit square, filled regularly.
Definition: gridGeneration.hh:431
Kaskade::createPentagon
std::unique_ptr< Grid > createPentagon(double radius=1)
Creates an regular pentagon consisting of five triangles, centered at the origin and with radius r.
Definition: gridGeneration.hh:775
Kaskade::createUnitCube
std::unique_ptr< Grid > createUnitCube(double dh=1., bool symmetric=true)
Creates a unit cube.
Definition: gridGeneration.hh:443
Kaskade::createCuboid
std::unique_ptr< Grid > createCuboid(Dune::FieldVector< Scalar, dim > const &x0, Dune::FieldVector< Scalar, dim > const &dx, EdgeLength const &dh, bool symmetric=true)
Creates a uniform simplicial grid on a rectangular or cuboid domain.
Definition: gridGeneration.hh:414
Kaskade::createUshape
std::unique_ptr< Grid > createUshape(double l, double d, double eps)
Creates an U-shaped domain.
Definition: gridGeneration.hh:830
Kaskade::createCylinder
std::unique_ptr< Grid > createCylinder(double radius, double height, int layers, int circumEdges=6)
Creates a cylinder.
Definition: gridGeneration.hh:683
Kaskade::createTwoLayerCuboid
std::unique_ptr< Grid > createTwoLayerCuboid(Dune::FieldVector< Scalar, dim > const &x0, Dune::FieldVector< Scalar, dim > const &dx1, Dune::FieldVector< Scalar, dim > const &dx2, Dune::FieldVector< Scalar, dim > const &dh1, Scalar dh2, bool symmetric=true)
Definition: gridGeneration.hh:505
Kaskade::createRectangle
std::unique_ptr< Grid > createRectangle(Dune::FieldVector< Scalar, 2 > const &x0, Dune::FieldVector< Scalar, 2 > const &dx, Scalar dh, bool symmetric=true)
Fill rectangle with squares.
Definition: gridGeneration.hh:460
Kaskade::createTShape
std::unique_ptr< Grid > createTShape(Dune::FieldVector< Scalar, dim > const &x0, Scalar dx, Scalar dy, Dune::FieldVector< Scalar, dim > const &thickness, Scalar dh, bool symmetric=true)
Fill L-shaped domain with cubes.
Definition: gridGeneration.hh:496
Kaskade::createGrid
std::unique_ptr< Grid > createGrid(std::vector< Type > const &cubes)
Extract simplicial grid from list of cubes.
Definition: gridGeneration.hh:371
Kaskade::createIcosahedron
std::unique_ptr< Grid > createIcosahedron()
Creates an icosahedron consisting of 20 tetrahedra.
Definition: gridGeneration.hh:609
Kaskade::createOctahedron
std::unique_ptr< Grid > createOctahedron(Dune::FieldVector< typename Grid::ctype, 3 > extent=Dune::FieldVector< typename Grid::ctype, 3 >(1.0))
Creates an octahedron consisting of eight tetrahedra.
Definition: gridGeneration.hh:574
Kaskade::createHexahedron
std::unique_ptr< Grid > createHexahedron()
Creates a hexahedron consisting of five tetrahedra.
Definition: gridGeneration.hh:539
Generated on Fri Sep 6 2024 13:13:00 for KASKADE 7 development version by doxygen 1.9.4