1
+ #include " mockgrid.h"
2
+ void create_grid (gmgpolar& test_p)
3
+ {
4
+ time_t seed = 5000 ;
5
+ std::default_random_engine gen (seed); // deterministic seed to reproduce
6
+ std::uniform_real_distribution<double > dis (gyro::dcntl[Param::R0], gyro::dcntl[Param::R]);
7
+ std::uniform_real_distribution<double > theta_distribution (0 , 2 * PI);
8
+
9
+ std::cout << " creating randomized grid with std::default_random_engine seed " + std::to_string (seed) << std::endl;
10
+ level* new_level = new level (0 );
11
+ new_level->nr = pow (2 , gyro::icntl[Param::nr_exp]);
12
+ new_level->r = std::vector<double >(new_level->nr + 1 );
13
+
14
+ std::vector<double > rands (new_level->nr - 1 );
15
+ for (int j = 0 ; j < new_level->nr - 1 ; ++j) {
16
+ rands[j] = dis (gen);
17
+ }
18
+ std::sort (rands.begin (), rands.end ());
19
+ new_level->r [0 ] = gyro::dcntl[Param::R0];
20
+ for (int i = 1 ; i < new_level->nr ; i++) {
21
+ new_level->r [i] = rands[i - 1 ]; // random grid on coarser level
22
+ }
23
+ new_level->r [new_level->nr ] = gyro::dcntl[Param::R];
24
+ new_level->nr ++;
25
+ int ntmp = pow (2 , ceil (log2 (new_level->nr )));
26
+ new_level->ntheta = gyro::icntl[Param::periodic] ? ntmp : ntmp + 1 ;
27
+
28
+ new_level->theta = std::vector<double >(new_level->ntheta );
29
+
30
+ std::vector<double > randst (ntmp - 1 );
31
+ for (int k = 0 ; k < ntmp - 1 ; ++k) {
32
+ randst[k] = theta_distribution (gen);
33
+ }
34
+ std::sort (randst.begin (), randst.end ());
35
+ new_level->theta [0 ] = 0 ;
36
+ for (int i = 1 ; i < ntmp; i++) {
37
+ new_level->theta [i] = randst[i - 1 ];
38
+ }
39
+ if (!gyro::icntl[Param::periodic]) {
40
+ new_level->theta [ntmp] = 2 * PI;
41
+ }
42
+
43
+ new_level->ntheta_int = gyro::icntl[Param::periodic] ? new_level->ntheta : new_level->ntheta - 1 ;
44
+ new_level->nr_int = new_level->nr - 1 ;
45
+
46
+ new_level->thetaplus = std::vector<double >(new_level->ntheta_int );
47
+ for (int k = 0 ; k < new_level->ntheta_int - 1 ; k++) {
48
+ new_level->thetaplus [k] = new_level->theta [k + 1 ] - new_level->theta [k];
49
+ }
50
+ new_level->thetaplus [new_level->ntheta_int - 1 ] = 2 * PI - new_level->theta [new_level->ntheta_int - 1 ];
51
+
52
+ new_level->hplus = std::vector<double >(new_level->nr_int );
53
+ for (int k = 0 ; k < new_level->nr_int ; k++) {
54
+ new_level->hplus [k] = new_level->r [k + 1 ] - new_level->r [k];
55
+ }
56
+
57
+ level* coarser_level = new level (1 );
58
+ coarser_level->nr = pow (2 , gyro::icntl[Param::nr_exp] - 1 ) + 1 ;
59
+ ntmp = pow (2 , ceil (log2 (coarser_level->nr )));
60
+ coarser_level->ntheta = gyro::icntl[Param::periodic] ? ntmp : ntmp + 1 ;
61
+
62
+ coarser_level->ntheta_int = gyro::icntl[Param::periodic] ? coarser_level->ntheta : coarser_level->ntheta - 1 ;
63
+ coarser_level->nr_int = coarser_level->nr - 1 ;
64
+ test_p.v_level .push_back (new_level);
65
+ test_p.v_level .push_back (coarser_level);
66
+ test_p.levels_orig = 2 ;
67
+ }
0 commit comments