blob: 667a38643b0a69187a945b1007c8bf68ffbe81fc (
plain)
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
|
#include <cassert>
#include <cmath>
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
#include "mask_surface.hh"
namespace CarpetMask {
using namespace std;
/**
* Set the weight in the excluded regions to zero.
*/
void
CarpetSurfaceSetup (CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
CCTK_REAL * const weight =
static_cast <CCTK_REAL *>
(CCTK_VarDataPtr (cctkGH, 0, "CarpetReduce::weight"));
if (not weight) {
CCTK_WARN (CCTK_WARN_ABORT,
"CarpetReduce is not active, or CarpetReduce::mask does not have storage");
}
for (int n = 0; n < 10; ++ n) {
int const sn = excluded_surface[n];
if (sn >= 0) {
assert (sn < nsurfaces);
CCTK_REAL const shrink_factor = excluded_surface_factor[n];
if (sf_valid[sn] > 0) {
CCTK_REAL const x0 = sf_origin_x[sn];
CCTK_REAL const y0 = sf_origin_y[sn];
CCTK_REAL const z0 = sf_origin_z[sn];
CCTK_REAL const theta0 = sf_origin_theta[sn];
CCTK_REAL const phi0 = sf_origin_phi [sn];
CCTK_REAL const dtheta = sf_delta_theta[sn];
CCTK_REAL const dphi = sf_delta_phi [sn];
int const ntheta = sf_ntheta[sn];
int const nphi = sf_nphi [sn];
for (int k = 0; k < cctk_lsh[2]; ++ k) {
for (int j = 0; j < cctk_lsh[1]; ++ j) {
for (int i = 0; i < cctk_lsh[0]; ++ i) {
int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
CCTK_REAL const dx = x[ind] - x0;
CCTK_REAL const dy = y[ind] - y0;
CCTK_REAL const dz = z[ind] - z0;
CCTK_REAL const rho =
sqrt (pow (dx, 2) + pow (dy, 2) + pow (dz, 2));
if (rho < 1.0e-12) {
// Always excise the surface origin
weight[ind] = 0.0;
} else {
CCTK_REAL const theta = acos (dz / rho);
CCTK_REAL const phi = atan2 (dy, dx);
int const a = floor ((theta - theta0) / dtheta + 0.5);
int const b = floor ((phi - phi0 ) / dphi + 0.5);
assert (a >= 0 and a < ntheta);
assert (b >= 0 and b < nphi );
CCTK_REAL const dr = sf_radius[a + maxntheta * b];
if (rho <= dr * shrink_factor) {
weight[ind] = 0.0;
}
}
}
}
}
} // if surface is valid
} // if r>=0
} // for n
}
} // namespace CarpetMask
|
