1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
|
#include <assert.h>
#include "cctk.h"
#include "cctk_Parameters.h"
#include "gh.hh"
#include "carpet.hh"
#include "regrid.hh"
namespace CarpetRegrid {
using namespace std;
using namespace Carpet;
int Centre (cGH const * const cctkGH,
gh const & hh,
gh::mregs & regsss)
{
DECLARE_CCTK_PARAMETERS;
assert (refinement_levels >= 1);
// do nothing if the levels already exist
if (reflevel == refinement_levels) return 0;
assert (regsss.size() >= 1);
vector<vector<region_t> > regss = regsss.at(0);
regss.resize (refinement_levels);
bvect const symmetric (symmetry_x, symmetry_y, symmetry_z);
ivect const zero(0), one(1), two(2);
ivect rstr = hh.baseextent.stride();
ivect rlb = hh.baseextent.lower();
ivect rub = hh.baseextent.upper();
assert (! smart_outer_boundaries);
for (size_t rl=1; rl<regss.size(); ++rl) {
// save old values
ivect const oldrlb = rlb;
ivect const oldrub = rub;
// refined boxes have smaller stride
assert (all(rstr%(hh.reffacts.at(rl)/hh.reffacts.at(rl-1)) == 0));
rstr /= hh.reffacts.at(rl)/hh.reffacts.at(rl-1);
// calculate new extent
ivect const quarter = (rub - rlb) / 4 / rstr * rstr;
ivect const half = (rub - rlb) / 2 / rstr * rstr;
rlb = oldrlb + symmetric.ifthen(zero, quarter);
rub = oldrub - symmetric.ifthen(half, quarter);
assert (all(rlb >= oldrlb and rub <= oldrub));
vector<region_t> regs (1);
ibbox const ext (rlb, rub, rstr);
regs.at(0).extent = ext;
b2vect const ob (false);
regs.at(0).outer_boundaries = ob;
b2vect const rb (true);
regs.at(0).refinement_boundaries = rb;
regs.at(0).map = Carpet::map;
// make multiprocessor aware
SplitRegions (cctkGH, regs);
regss.at(rl) = regs;
} // for rl
// make multigrid aware
MakeMultigridBoxes (cctkGH, regss, regsss);
return 1;
}
} // namespace CarpetRegrid
|
