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
123
124
125
126
127
128
|
/* $Header:$ */
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
void
CoordBase_SetupMask (CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
CCTK_INT nboundaryzones[6];
CCTK_INT is_internal[6];
CCTK_INT is_staggered[6];
CCTK_INT shiftout[6];
int imin[3], imax[3]; /* boundary extent */
int i, j, k;
int d, f;
int dd;
int ierr;
ierr = GetBoundarySpecification
(6, nboundaryzones, is_internal, is_staggered, shiftout);
if (ierr != 0) {
CCTK_WARN (0, "Could not get boundary specification");
}
/* Loop over all dimensions and faces */
for (d=0; d<3; ++d) {
for (f=0; f<2; ++f) {
/* If this processor has the outer boundary */
if (cctk_bbox[2*d+f]) {
int npoints;
if (is_internal[2*d+f]) {
/* The boundary extends inwards */
npoints = - nboundaryzones[2*d+f] + shiftout[2*d+f];
} else {
/* The boundary extends outwards */
npoints = nboundaryzones[2*d+f] + shiftout[2*d+f] - 1;
}
/* If there are boundary that should be ignored */
if (npoints >= 0) {
if (npoints < 0 || npoints > cctk_lsh[d]) {
CCTK_WARN (0, "Illegal number of boundary points");
}
/* Calculate the extent of the boundary */
for (dd=0; dd<3; ++dd) {
imin[dd] = 0;
imax[dd] = cctk_lsh[dd];
}
if (f==0) {
/* lower face */
imax[d] = imin[d] + npoints;
} else {
/* upper face */
imin[d] = imax[d] - npoints;
}
/* Loop over the boundary */
if (verbose) {
CCTK_VInfo (CCTK_THORNSTRING,
"Setting boundary points in direction %d face %d to weight 0", d, f);
}
for (k=imin[2]; k<imax[2]; ++k) {
for (j=imin[1]; j<imax[1]; ++j) {
for (i=imin[0]; i<imax[0]; ++i) {
int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
weight[ind] = 0.0;
}
}
}
/* When the boundary is not staggered, then give the points
on the boundary the weight 1/2 */
if (! is_staggered[2*d+f]) {
/* Adapt the extent of the boundary */
if (f==0) {
/* lower face */
imin[d] = imax[d];
imax[d] = imin[d] + 1;
} else {
/* upper face */
imax[d] = imin[d];
imin[d] = imax[d] - 1;
}
/* Loop over the points next to boundary */
if (verbose) {
CCTK_VInfo (CCTK_THORNSTRING,
"Setting staggered boundary points in direction %d face %d to weight 1/2", d, f);
}
for (k=imin[2]; k<imax[2]; ++k) {
for (j=imin[1]; j<imax[1]; ++j) {
for (i=imin[0]; i<imax[0]; ++i) {
int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
weight[ind] *= 0.5;
}
}
}
} /* if the boundary is not staggered */
} /* if there are boundary points */
} /* if is outer boundary */
} /* loop over faces */
} /* loop over directions */
}
|