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
89
90
|
#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<dim> const & hh,
gh<dim>::rexts & bbsss,
gh<dim>::rbnds & obss,
gh<dim>::rprocs & pss)
{
DECLARE_CCTK_PARAMETERS;
assert (refinement_levels >= 1);
// do nothing if the levels already exist
if (reflevel == refinement_levels) return 0;
assert (bbsss.size() >= 1);
bbsss.resize (refinement_levels);
obss.resize (refinement_levels);
pss.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<bbsss.size(); ++rl) {
// save old values
ivect const oldrlb = rlb;
ivect const oldrub = rub;
// refined boxes have smaller stride
assert (all(rstr%hh.reffact == 0));
rstr /= hh.reffact;
// 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 && rub <= oldrub));
ibbox const bb (rlb, rub, rstr);
vector<ibbox> bbs (1);
bbs.at(0) = bb;
bbvect const ob (false);
gh<dim>::cbnds obs (1);
obs.at(0) = ob;
// make multiprocessor aware
gh<dim>::cprocs ps;
SplitRegions (cctkGH, bbs, obs, ps);
// make multigrid aware
vector<vector<ibbox> > bbss;
MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
bbsss.at(rl) = bbss;
obss.at(rl) = obs;
pss.at(rl) = ps;
} // for rl
return 1;
}
} // namespace CarpetRegrid
|
