-
Notifications
You must be signed in to change notification settings - Fork 4
Expand file tree
/
Copy pathapplyAGive.cpp
More file actions
304 lines (274 loc) · 15.3 KB
/
applyAGive.cpp
File metadata and controls
304 lines (274 loc) · 15.3 KB
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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
#include "../../../include/Residual/ResidualGive/residualGive.h"
static inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, const PolarGrid& grid,
const LevelCache& level_cache, bool DirBC_Interior, Vector<double>& result,
ConstVector<double>& x)
{
/* ---------------------------------------- */
/* Compute or retrieve stencil coefficients */
/* ---------------------------------------- */
const int center = grid.index(i_r, i_theta);
double coeff_beta, arr, att, art, detDF;
level_cache.obtainValues(i_r, i_theta, center, r, theta, coeff_beta, arr, att, art, detDF);
/* -------------------- */
/* Node in the interior */
/* -------------------- */
if (i_r > 1 && i_r < grid.nr() - 2) {
double h1 = grid.radialSpacing(i_r - 1);
double h2 = grid.radialSpacing(i_r);
double k1 = grid.angularSpacing(i_theta - 1);
double k2 = grid.angularSpacing(i_theta);
double coeff1 = 0.5 * (k1 + k2) / h1;
double coeff2 = 0.5 * (k1 + k2) / h2;
double coeff3 = 0.5 * (h1 + h2) / k1;
double coeff4 = 0.5 * (h1 + h2) / k2;
const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
const int left = grid.index(i_r - 1, i_theta);
const int right = grid.index(i_r + 1, i_theta);
const int bottom = grid.index(i_r, i_theta_M1);
const int top = grid.index(i_r, i_theta_P1);
/* Fill result(i,j) */
result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */
- coeff1 * arr * x[left] /* Left */
- coeff2 * arr * x[right] /* Right */
- coeff3 * att * x[bottom] /* Bottom */
- coeff4 * att * x[top] /* Top */
+ ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) *
x[center]) /* Center: (Left, Right, Bottom, Top) */;
/* Fill result(i-1,j) */
result[left] += (-coeff1 * arr * x[center] /* Right */
+ coeff1 * arr * x[left] /* Center: (Right) */
- 0.25 * art * x[top] /* Top Right */
+ 0.25 * art * x[bottom]); /* Bottom Right */
/* Fill result(i+1,j) */
result[right] += (-coeff2 * arr * x[center] /* Left */
+ coeff2 * arr * x[right] /* Center: (Left) */
+ 0.25 * art * x[top] /* Top Left */
- 0.25 * art * x[bottom]); /* Bottom Left */
/* Fill result(i,j-1) */
result[bottom] += (-coeff3 * att * x[center] /* Top */
+ coeff3 * att * x[bottom] /* Center: (Top) */
- 0.25 * art * x[right] /* Top Right */
+ 0.25 * art * x[left]); /* Top Left */
/* Fill result(i,j+1) */
result[top] += (-coeff4 * att * x[center] /* Bottom */
+ coeff4 * att * x[top] /* Center: (Bottom) */
+ 0.25 * art * x[right] /* Bottom Right */
- 0.25 * art * x[left]); /* Bottom Left */
}
/* -------------------------- */
/* Node on the inner boundary */
/* -------------------------- */
else if (i_r == 0) {
/* ------------------------------------------------ */
/* Case 1: Dirichlet boundary on the inner boundary */
/* ------------------------------------------------ */
if (DirBC_Interior) {
/* Fill result(i,j) */
result[center] += x[center];
/* Give value to the interior nodes! */
double h2 = grid.radialSpacing(i_r);
double k1 = grid.angularSpacing(i_theta - 1);
double k2 = grid.angularSpacing(i_theta);
double coeff2 = 0.5 * (k1 + k2) / h2;
const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
const int right = grid.index(i_r + 1, i_theta);
const int bottom = grid.index(i_r, i_theta_M1);
const int top = grid.index(i_r, i_theta_P1);
/* Fill result(i+1,j) */
result[right] += (-coeff2 * arr * x[center] /* Left */
+ coeff2 * arr * x[right] /* Center: (Left) */
+ 0.25 * art * x[top] /* Top Left */
- 0.25 * art * x[bottom]); /* Bottom Left */
}
else {
/* ------------------------------------------------------------- */
/* Case 2: Across origin discretization on the interior boundary */
/* ------------------------------------------------------------- */
/* h1 gets replaced with 2 * R0. */
/* (i_r-1,i_theta) gets replaced with (i_r, i_theta + (grid.ntheta()/2)). */
/* Some more adjustments from the changing the 9-point stencil to the artifical 7-point stencil. */
double h1 = 2.0 * grid.radius(0);
double h2 = grid.radialSpacing(i_r);
double k1 = grid.angularSpacing(i_theta - 1);
double k2 = grid.angularSpacing(i_theta);
double coeff1 = 0.5 * (k1 + k2) / h1;
double coeff2 = 0.5 * (k1 + k2) / h2;
double coeff3 = 0.5 * (h1 + h2) / k1;
double coeff4 = 0.5 * (h1 + h2) / k2;
const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
const int i_theta_Across = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2);
const int left = grid.index(i_r, i_theta_Across);
const int right = grid.index(i_r + 1, i_theta);
const int bottom = grid.index(i_r, i_theta_M1);
const int top = grid.index(i_r, i_theta_P1);
/* Fill result(i,j) */
result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */
- coeff1 * arr * x[left] /* Left */
- coeff2 * arr * x[right] /* Right */
- coeff3 * att * x[bottom] /* Bottom */
- coeff4 * att * x[top] /* Top */
+ ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) *
x[center]); /* Center: (Left, Right, Bottom, Top) */
/* Fill result(i-1,j) */
/* From view the view of the across origin node, the directions are roatated by 180 degrees in the stencil! */
result[left] += (-coeff1 * arr * x[center] /* Right -> Left */
+ coeff1 * arr * x[left]); /* Center: (Right) -> Center: (Left)*/
/* + 0.25 * art * x[top]; // Top Right -> Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
/* - 0.25 * art * x[bottom]; // Bottom Right -> Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
/* Fill result(i+1,j) */
result[right] += (-coeff2 * arr * x[center] /* Left */
+ coeff2 * arr * x[right] /* Center: (Left) */
+ 0.25 * art * x[top] /* Top Left */
- 0.25 * art * x[bottom]); /* Bottom Left */
/* Fill result(i,j-1) */
result[bottom] += (-coeff3 * att * x[center] /* Top */
+ coeff3 * att * x[bottom] /* Center: (Top) */
- 0.25 * art * x[right]); /* Top Right */
/* + 0.25 * art * x[left]; // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
/* Fill result(i,j+1) */
result[top] += (-coeff4 * att * x[center] /* Bottom */
+ coeff4 * att * x[top] /* Center: (Bottom) */
+ 0.25 * art * x[right]); /* Bottom Right */
/* - 0.25 * art * x[left]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
}
}
/* ------------------------------- */
/* Node next to the inner boundary */
/* ------------------------------- */
else if (i_r == 1) {
double h1 = grid.radialSpacing(i_r - 1);
double h2 = grid.radialSpacing(i_r);
double k1 = grid.angularSpacing(i_theta - 1);
double k2 = grid.angularSpacing(i_theta);
double coeff1 = 0.5 * (k1 + k2) / h1;
double coeff2 = 0.5 * (k1 + k2) / h2;
double coeff3 = 0.5 * (h1 + h2) / k1;
double coeff4 = 0.5 * (h1 + h2) / k2;
const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
const int left = grid.index(i_r - 1, i_theta);
const int right = grid.index(i_r + 1, i_theta);
const int bottom = grid.index(i_r, i_theta_M1);
const int top = grid.index(i_r, i_theta_P1);
/* Fill result(i,j) */
result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */
- coeff1 * arr * x[left] /* Left */
- coeff2 * arr * x[right] /* Right */
- coeff3 * att * x[bottom] /* Bottom */
- coeff4 * att * x[top] /* Top */
+ ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) *
x[center]); /* Center: (Left, Right, Bottom, Top) */
/* Fill result(i-1,j) */
if (!DirBC_Interior) { /* Don't give to the inner dirichlet boundary! */
result[left] += (-coeff1 * arr * x[center] /* Right */
+ coeff1 * arr * x[left] /* Center: (Right) */
- 0.25 * art * x[top] /* Top Right */
+ 0.25 * art * x[bottom]); /* Bottom Right */
}
/* Fill result(i+1,j) */
result[right] += (-coeff2 * arr * x[center] /* Left */
+ coeff2 * arr * x[right] /* Center: (Left) */
+ 0.25 * art * x[top] /* Top Left */
- 0.25 * art * x[bottom]); /* Bottom Left */
/* Fill result(i,j-1) */
result[bottom] += (-coeff3 * att * x[center] /* Top */
+ coeff3 * att * x[bottom] /* Center: (Top) */
- 0.25 * art * x[right] /* Top Right */
+ 0.25 * art * x[left]); /* Top Left */
/* Fill result(i,j+1) */
result[top] += (-coeff4 * att * x[center] /* Bottom */
+ coeff4 * att * x[top] /* Center: (Bottom) */
+ 0.25 * art * x[right] /* Bottom Right */
- 0.25 * art * x[left]); /* Bottom Left */
}
/* ------------------------------- */
/* Node next to the outer boundary */
/* ------------------------------- */
else if (i_r == grid.nr() - 2) {
double h1 = grid.radialSpacing(i_r - 1);
double h2 = grid.radialSpacing(i_r);
double k1 = grid.angularSpacing(i_theta - 1);
double k2 = grid.angularSpacing(i_theta);
double coeff1 = 0.5 * (k1 + k2) / h1;
double coeff2 = 0.5 * (k1 + k2) / h2;
double coeff3 = 0.5 * (h1 + h2) / k1;
double coeff4 = 0.5 * (h1 + h2) / k2;
const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
const int left = grid.index(i_r - 1, i_theta);
const int right = grid.index(i_r + 1, i_theta);
const int bottom = grid.index(i_r, i_theta_M1);
const int top = grid.index(i_r, i_theta_P1);
/* Fill result(i,j) */
result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */
- coeff1 * arr * x[left] /* Left */
- coeff2 * arr * x[right] /* Right */
- coeff3 * att * x[bottom] /* Bottom */
- coeff4 * att * x[top] /* Top */
+ ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) *
x[center]); /* Center: (Left, Right, Bottom, Top) */
/* Fill result(i-1,j) */
result[left] += (-coeff1 * arr * x[center] /* Right */
+ coeff1 * arr * x[left] /* Center: (Right) */
- 0.25 * art * x[top] /* Top Right */
+ 0.25 * art * x[bottom]); /* Bottom Right */
/* Don't give to the outer dirichlet boundary! */
/* Fill result(i+1,j) */
/* result[right] -= ( */
/* - coeff2 * arr * x[center] // Left */
/* + coeff2 * arr * x[right] // Center: (Left) */
/* + 0.25 * art * x[top] // Top Left */
/* - 0.25 * art * x[bottom]); // Bottom Left */
/* Fill result(i,j-1) */
result[bottom] += (-coeff3 * att * x[center] /* Top */
+ coeff3 * att * x[bottom] /* Center: (Top) */
- 0.25 * art * x[right] /* Top Right */
+ 0.25 * art * x[left]); /* Top Left */
/* Fill result(i,j+1) */
result[top] += (-coeff4 * att * x[center] /* Bottom */
+ coeff4 * att * x[top] /* Center: (Bottom) */
+ 0.25 * art * x[right] /* Bottom Right */
- 0.25 * art * x[left]); /* Bottom Left */
}
/* ----------------------------- */
/* Node on to the outer boundary */
/* ----------------------------- */
else if (i_r == grid.nr() - 1) {
/* Fill result of (i,j) */
result[center] += x[center];
/* Give value to the interior nodes! */
double h1 = grid.radialSpacing(i_r - 1);
double k1 = grid.angularSpacing(i_theta - 1);
double k2 = grid.angularSpacing(i_theta);
double coeff1 = 0.5 * (k1 + k2) / h1;
const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
const int left = grid.index(i_r - 1, i_theta);
const int bottom = grid.index(i_r, i_theta_M1);
const int top = grid.index(i_r, i_theta_P1);
/* Fill result(i-1,j) */
result[left] += (-coeff1 * arr * x[center] /* Right */
+ coeff1 * arr * x[left] /* Center: (Right) */
- 0.25 * art * x[top] /* Top Right */
+ 0.25 * art * x[bottom]); /* Bottom Right */
}
}
void ResidualGive::applyCircleSection(const int i_r, Vector<double> result, ConstVector<double> x) const
{
const double r = grid_.radius(i_r);
for (int i_theta = 0; i_theta < grid_.ntheta(); i_theta++) {
const double theta = grid_.theta(i_theta);
node_apply_a_give(i_r, i_theta, r, theta, grid_, level_cache_, DirBC_Interior_, result, x);
}
}
void ResidualGive::applyRadialSection(const int i_theta, Vector<double> result, ConstVector<double> x) const
{
const double theta = grid_.theta(i_theta);
for (int i_r = grid_.numberSmootherCircles(); i_r < grid_.nr(); i_r++) {
const double r = grid_.radius(i_r);
node_apply_a_give(i_r, i_theta, r, theta, grid_, level_cache_, DirBC_Interior_, result, x);
}
}