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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
#include <assert.h>
#include <sstream>
#include <vector>
#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 ManualGridpoints (cGH const * const cctkGH,
gh const & hh,
gh::mregs & regsss)
{
DECLARE_CCTK_PARAMETERS;
if (refinement_levels > 4) {
CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels");
}
assert (refinement_levels >= 1 and refinement_levels <= 4);
// 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);
vector<ivect> ilower(3), iupper(3);
ilower.at(0) = ivect (l1ixmin, l1iymin, l1izmin);
iupper.at(0) = ivect (l1ixmax, l1iymax, l1izmax);
ilower.at(1) = ivect (l2ixmin, l2iymin, l2izmin);
iupper.at(1) = ivect (l2ixmax, l2iymax, l2izmax);
ilower.at(2) = ivect (l3ixmin, l3iymin, l3izmin);
iupper.at(2) = ivect (l3ixmax, l3iymax, l3izmax);
assert (! smart_outer_boundaries);
for (size_t rl=1; rl<regss.size(); ++rl) {
b2vect const ob (false);
b2vect const rb (true);
region_t reg;
reg.map = Carpet::map;
reg.outer_boundaries = ob;
reg.refinement_boundaries = rb;
vector<region_t> regs;
ManualGridpoints_OneLevel
(cctkGH, hh, rl,refinement_levels,
ilower.at(rl-1), iupper.at(rl-1), reg, regs);
// make multiprocessor aware
SplitRegions (cctkGH, regs);
regss.at(rl) = regs;
} // for rl
// make multigrid aware
MakeMultigridBoxes (cctkGH, regss, regsss);
return 1;
}
void ManualGridpoints_OneLevel (const cGH * const cctkGH,
const gh & hh,
const int rl,
const int numrl,
const ivect ilower,
const ivect iupper,
const region_t & reg,
vector<region_t> & regs)
{
const ivect rstr = hh.baseextent.stride();
const ivect rlb = hh.baseextent.lower();
const ivect rub = hh.baseextent.upper();
const ivect levfac = hh.reffacts.at(rl);
assert (all (rstr % levfac == 0));
const ivect str (rstr / levfac);
const ivect lb (ilower);
const ivect ub (iupper);
if (! all(lb>=rlb and ub<=rub)) {
ostringstream buf;
buf << "The refinement region boundaries for refinement level #" << rl << " are not within the main grid. Allowed are the grid point boundaries " << rlb << " - " << rub << "; specified were " << lb << " - " << ub << ends;
CCTK_WARN (0, buf.str().c_str());
}
if (! all(lb<=ub)) {
ostringstream buf;
buf << "The refinement region boundaries for refinement level #" << rl << " have the upper boundary (" << ub << ") less than the lower boundary (" << lb << ")" << ends;
CCTK_WARN (0, buf.str().c_str());
}
if (! all(lb%str==0 and ub%str==0)) {
CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
"The refinement region boundaries for refinement level #%d are not a multiple of the stride for that level", rl);
}
assert (all(lb>=rlb and ub<=rub));
assert (all(lb<=ub));
assert (all(lb%str==0 and ub%str==0));
region_t newreg (reg);
newreg.extent = ibbox(lb, ub, str);
regs.push_back (newreg);
}
} // namespace CarpetRegrid
|
