APBS 3.0.0
Loading...
Searching...
No Matches
matvecd.c
1
55#include "matvecd.h"
56
57VPUBLIC void Vmatvec(int *nx, int *ny, int *nz,
58 int *ipc, double *rpc,
59 double *ac, double *cc,
60 double *x, double *y) {
61
62 int numdia;
63
64 // Do in one step
65 numdia = VAT(ipc, 11);
66
67 if (numdia == 7) {
68 Vmatvec7(nx, ny, nz,
69 ipc, rpc,
70 ac, cc,
71 x, y);
72 } else if (numdia == 27) {
73 Vmatvec27(nx, ny, nz,
74 ipc, rpc,
75 ac, cc,
76 x, y);
77 } else {
78 Vnm_print(2, "MATVEC: invalid stencil type given...");
79 }
80}
81
82VPUBLIC void Vmatvec7(int *nx, int *ny, int *nz,
83 int *ipc, double *rpc,
84 double *ac, double *cc,
85 double *x, double *y) {
86
87 MAT2(ac, *nx * *ny * *nz, 1);
88
89 Vmatvec7_1s(nx, ny, nz,
90 ipc, rpc,
91 RAT2(ac, 1, 1), cc,
92 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
93 x, y);
94}
95
96
97
98VEXTERNC void Vmatvec7_1s(int *nx, int *ny, int *nz,
99 int *ipc, double *rpc,
100 double *oC, double *cc,
101 double *oE, double *oN, double *uC,
102 double *x, double *y) {
103
104 int i, j, k;
105
106 MAT3(oE, *nx, *ny, *nz);
107 MAT3(oN, *nx, *ny, *nz);
108 MAT3(uC, *nx, *ny, *nz);
109 MAT3(cc, *nx, *ny, *nz);
110 MAT3(oC, *nx, *ny, *nz);
111 MAT3(x, *nx, *ny, *nz);
112 MAT3(y, *nx, *ny, *nz);
113
114 // Do it
115 #pragma omp parallel for private(i, j, k)
116 for (k=2; k<=*nz-1; k++) {
117 for (j=2; j<=*ny-1; j++) {
118 for(i=2; i<=*nx-1; i++) {
119 VAT3(y, i, j, k) =
120 - VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
121 - VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
122 - VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
123 - VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
124 - VAT3( uC, i, j, k-1) * VAT3(x, i, j,k-1)
125 - VAT3( uC, i, j, k) * VAT3(x, i, j,k+1)
126 + (VAT3(oC, i, j, k) + VAT3(cc, i, j, k)) * VAT3(x, i, j, k);
127 }
128 }
129 }
130}
131
132
133
134VPUBLIC void Vmatvec27(int *nx, int *ny, int *nz,
135 int *ipc, double *rpc,
136 double *ac, double *cc,
137 double *x, double *y) {
138
139 MAT2(ac, *nx * *ny * *nz, 1);
140
141 Vmatvec27_1s(nx, ny, nz,
142 ipc, rpc,
143 RAT2(ac, 1, 1), cc,
144 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
145 RAT2(ac, 1, 5), RAT2(ac, 1, 6),
146 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
147 RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
148 x, y);
149}
150
151
152
153VPUBLIC void Vmatvec27_1s(int *nx, int *ny, int *nz,
154 int *ipc, double *rpc,
155 double *oC, double *cc,
156 double *oE, double *oN, double *uC,
157 double *oNE, double *oNW,
158 double *uE, double *uW, double *uN, double *uS,
159 double *uNE, double *uNW, double *uSE, double *uSW,
160 double *x, double *y) {
161
162 int i, j, k;
163
164 double tmpO, tmpU, tmpD;
165
166 MAT3(cc, *nx, *ny, *nz);
167 MAT3(x, *nx, *ny, *nz);
168 MAT3(y, *nx, *ny, *nz);
169
170 MAT3(oC, *nx, *ny, *nz);
171 MAT3(oE, *nx, *ny, *nz);
172 MAT3(oN, *nx, *ny, *nz);
173 MAT3(oNE, *nx, *ny, *nz);
174 MAT3(oNW, *nx, *ny, *nz);
175
176 MAT3(uC, *nx, *ny, *nz);
177 MAT3(uE, *nx, *ny, *nz);
178 MAT3(uW, *nx, *ny, *nz);
179 MAT3(uN, *nx, *ny, *nz);
180 MAT3(uS, *nx, *ny, *nz);
181 MAT3(uNE, *nx, *ny, *nz);
182 MAT3(uNW, *nx, *ny, *nz);
183 MAT3(uSE, *nx, *ny, *nz);
184 MAT3(uSW, *nx, *ny, *nz);
185
186 // Do it
187 #pragma omp parallel for private(i, j, k, tmpO, tmpU, tmpD)
188 for (k=2; k<=*nz-1; k++) {
189 for (j=2; j<=*ny-1; j++) {
190 for(i=2; i<=*nx-1; i++) {
191 tmpO =
192 - VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
193 - VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
194 - VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
195 - VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
196 - VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
197 - VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
198 - VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
199 - VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
200
201 tmpU =
202 - VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
203 - VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
204 - VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
205 - VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
206 - VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
207 - VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
208 - VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
209 - VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
210 - VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
211
212 tmpD =
213 - VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
214 - VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
215 - VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
216 - VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
217 - VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
218 - VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
219 - VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
220 - VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
221 - VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
222
223 VAT3(y, i, j, k) = tmpO + tmpU + tmpD
224 + (VAT3(oC, i, j, k) + VAT3(cc, i, j, k)) * VAT3(x, i, j, k);
225 }
226 }
227 }
228}
229
230
231
232VEXTERNC void Vnmatvec(int *nx, int *ny, int *nz,
233 int *ipc, double *rpc,
234 double *ac, double *cc, double *x, double *y, double *w1) {
235
236 int numdia;
237
238 // Do in one step
239 numdia = VAT(ipc, 11);
240
241 if (numdia == 7) {
242 Vnmatvec7(nx, ny, nz,
243 ipc, rpc,
244 ac, cc,
245 x, y, w1);
246 } else if (numdia == 27) {
247 Vnmatvec27(nx, ny, nz,
248 ipc, rpc,
249 ac, cc,
250 x, y, w1);
251 } else {
252 Vnm_print(2, "MATVEC: invalid stencil type given...");
253 }
254}
255
256
257
258VPUBLIC void Vnmatvec7(int *nx, int *ny, int *nz,
259 int *ipc, double *rpc,
260 double *ac, double *cc,
261 double *x, double *y, double *w1) {
262
263 MAT2(ac, *nx * *ny * *nz, 1);
264
265 WARN_UNTESTED;
266
267 Vnmatvecd7_1s(nx, ny, nz,
268 ipc, rpc,
269 RAT2(ac, 1, 1), cc,
270 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
271 x, y, w1);
272}
273
274
275
276VPUBLIC void Vnmatvecd7_1s(int *nx, int *ny, int *nz,
277 int *ipc, double *rpc,
278 double *oC, double *cc,
279 double *oE, double *oN, double *uC,
280 double *x, double *y, double *w1) {
281
282 int i, j, k;
283 int ipkey;
284
285 MAT3(oE, *nx, *ny, *nz);
286 MAT3(oN, *nx, *ny, *nz);
287 MAT3(uC, *nx, *ny, *nz);
288 MAT3(cc, *nx, *ny, *nz);
289 MAT3(oC, *nx, *ny, *nz);
290 MAT3( x, *nx, *ny, *nz);
291 MAT3( y, *nx, *ny, *nz);
292 MAT3(w1, *nx, *ny, *nz);
293
294 WARN_UNTESTED;
295
296 // first get vector nonlinear term to avoid subroutine calls
297 ipkey = VAT(ipc, 10);
298 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
299
300 // The operator
301 #pragma omp parallel for private(i, j, k)
302 for (k=2; k<=*nz-1; k++)
303 for (j=2; j<=*ny-1; j++)
304 for(i=2; i<=*nx-1; i++)
305 VAT3(y, i, j, k) =
306 - VAT3(oN, i, j, k) * VAT3(x, i, j+1, k)
307 - VAT3(oN, i, j-1, k) * VAT3(x, i, j-1, k)
308 - VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
309 - VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
310 - VAT3(uC, i, j, k-1) * VAT3(x, i, j, k-1)
311 - VAT3(uC, i, j, k) * VAT3(x, i, j, k+1)
312 + VAT3(oC, i, j, k) * VAT3(x, i, j, k)
313 + VAT3(w1, i, j, k);
314}
315
316
317VPUBLIC void Vnmatvec27(int *nx, int *ny, int *nz,
318 int *ipc, double *rpc,
319 double *ac, double *cc,
320 double *x, double *y, double *w1) {
321
322 MAT2(ac, *nx * *ny * *nz, 1);
323
324 WARN_UNTESTED;
325
326 // Do in one step
327 Vnmatvecd27_1s(nx, ny, nz,
328 ipc, rpc,
329 RAT2(ac, 1, 1), cc,
330 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
331 RAT2(ac, 1, 5), RAT2(ac, 1, 6),
332 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
333 RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
334 x, y, w1);
335}
336
337
338
339VPUBLIC void Vnmatvecd27_1s(int *nx, int *ny, int *nz,
340 int *ipc, double *rpc,
341 double *oC, double *cc,
342 double *oE, double *oN, double *uC,
343 double *oNE, double *oNW,
344 double *uE, double *uW, double *uN, double *uS,
345 double *uNE, double *uNW, double *uSE, double *uSW,
346 double *x, double *y, double *w1) {
347
348 int i, j, k;
349 int ipkey;
350
351 double tmpO, tmpU, tmpD;
352
353 MAT3( oE, *nx, *ny, *nz);
354 MAT3( oN, *nx, *ny, *nz);
355 MAT3( uC, *nx, *ny, *nz);
356 MAT3(oNE, *nx, *ny, *nz);
357 MAT3(oNW, *nx, *ny, *nz);
358 MAT3( uE, *nx, *ny, *nz);
359 MAT3( uW, *nx, *ny, *nz);
360 MAT3( uN, *nx, *ny, *nz);
361 MAT3( uS, *nx, *ny, *nz);
362 MAT3(uNE, *nx, *ny, *nz);
363 MAT3(uNW, *nx, *ny, *nz);
364 MAT3(uSE, *nx, *ny, *nz);
365 MAT3(uSW, *nx, *ny, *nz);
366 MAT3( cc, *nx, *ny, *nz);
367 MAT3( oC, *nx, *ny, *nz);
368 MAT3( x, *nx, *ny, *nz);
369 MAT3( y, *nx, *ny, *nz);
370 MAT3( w1, *nx, *ny, *nz);
371
372 WARN_UNTESTED;
373
374 // First get vector noNlinear term to avoid subroutine calls
375 ipkey = VAT(ipc, 10);
376 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
377
378 // The operator
379 #pragma omp parallel for private(i, j, k, tmpO, tmpU, tmpD)
380 for (k=2; k<=*nz-1; k++) {
381 for (j=2; j<=*ny-1; j++) {
382 for(i=2; i<=*nx-1; i++) {
383
384 tmpO =
385 - VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
386 - VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
387 - VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
388 - VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
389 - VAT3(oNE, i, j, k) * VAT3(x, i+1, j+1, k)
390 - VAT3(oNW, i, j, k) * VAT3(x, i-1, j+1, k)
391 - VAT3(oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
392 - VAT3(oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
393
394 tmpU =
395 - VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
396 - VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
397 - VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
398 - VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
399 - VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
400 - VAT3(uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
401 - VAT3(uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
402 - VAT3(uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
403 - VAT3(uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
404
405 tmpD =
406 - VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
407 - VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
408 - VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
409 - VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
410 - VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
411 - VAT3(uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
412 - VAT3(uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
413 - VAT3(uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
414 - VAT3(uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
415
416 VAT3(y, i, j, k) = tmpO + tmpU + tmpD
417 + VAT3(oC, i, j, k) * VAT3(x, i, j, k)
418 + VAT3(w1, i, j, k);
419 }
420 }
421 }
422}
423
424
425
426VPUBLIC void Vmresid(int *nx, int *ny, int *nz,
427 int *ipc, double *rpc,
428 double *ac, double *cc, double *fc,
429 double *x, double *r) {
430
431 int numdia;
432
433 // Do in one step
434 numdia = VAT(ipc, 11);
435 if (numdia == 7) {
436 Vmresid7(nx, ny, nz, ipc, rpc, ac, cc, fc, x, r);
437 } else if (numdia == 27) {
438 Vmresid27(nx, ny, nz, ipc, rpc, ac, cc, fc, x, r);
439 } else {
440 Vnm_print(2, "Vmresid: invalid stencil type given...\n");
441 }
442}
443
444
445
446VPUBLIC void Vmresid7(int *nx, int *ny, int *nz,
447 int *ipc, double *rpc,
448 double *ac, double *cc, double *fc,
449 double *x, double *r) {
450
451 MAT2(ac, *nx * *ny * *nz, 1);
452
453 // Do in one step
454 Vmresid7_1s(nx, ny, nz,
455 ipc, rpc,
456 RAT2(ac, 1,1), cc, fc,
457 RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
458 x,r);
459}
460
461VPUBLIC void Vmresid7_1s(int *nx, int *ny, int *nz,
462 int *ipc, double *rpc,
463 double *oC, double *cc, double *fc,
464 double *oE, double *oN, double *uC,
465 double *x, double *r) {
466
467 int i, j, k;
468
469 MAT3(oE, *nx, *ny, *nz);
470 MAT3(oN, *nx, *ny, *nz);
471 MAT3(uC, *nx, *ny, *nz);
472 MAT3(cc, *nx, *ny, *nz);
473 MAT3(fc, *nx, *ny, *nz);
474 MAT3(oC, *nx, *ny, *nz);
475 MAT3(x, *nx, *ny, *nz);
476 MAT3(r, *nx, *ny, *nz);
477
478 // Do it
479 #pragma omp parallel for private(i, j, k)
480 for (k=2; k<=*nz-1; k++) {
481 for (j=2; j<=*ny-1; j++) {
482 for(i=2; i<=*nx-1; i++) {
483 VAT3(r, i,j,k) = VAT3(fc, i, j, k)
484 + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
485 + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
486 + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
487 + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
488 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
489 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
490 - (VAT3(oC, i, j, k) + VAT3(cc, i, j, k)) * VAT3(x, i, j, k);
491 }
492 }
493 }
494}
495
496
497
498VPUBLIC void Vmresid27(int *nx, int *ny, int *nz,
499 int *ipc, double *rpc,
500 double *ac, double *cc, double *fc,
501 double *x, double *r) {
502
503 MAT2(ac, *nx * *ny * *nz, 1);
504
505 // Do in one step
506 Vmresid27_1s(nx,ny,nz,
507 ipc, rpc,
508 RAT2(ac, 1, 1), cc, fc,
509 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
510 RAT2(ac, 1, 5), RAT2(ac, 1, 6),
511 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
512 RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
513 x,r);
514}
515
516
517
518VPUBLIC void Vmresid27_1s(int *nx, int *ny, int *nz,
519 int *ipc, double *rpc,
520 double *oC, double *cc, double *fc,
521 double *oE, double *oN, double *uC,
522 double *oNE, double *oNW,
523 double *uE, double *uW, double *uN, double *uS,
524 double *uNE, double *uNW, double *uSE, double *uSW,
525 double *x, double *r) {
526
527 int i, j, k;
528
529 double tmpO, tmpU, tmpD;
530
531 MAT3(cc, *nx, *ny, *nz);
532 MAT3(fc, *nx, *ny, *nz);
533 MAT3(x, *nx, *ny, *nz);
534 MAT3(r, *nx, *ny, *nz);
535
536 MAT3(oC, *nx, *ny, *nz);
537 MAT3(oE, *nx, *ny, *nz);
538 MAT3(oN, *nx, *ny, *nz);
539 MAT3(oNE, *nx, *ny, *nz);
540 MAT3(oNW, *nx, *ny, *nz);
541
542 MAT3(uC, *nx, *ny, *nz);
543 MAT3(uE, *nx, *ny, *nz);
544 MAT3(uW, *nx, *ny, *nz);
545 MAT3(uN, *nx, *ny, *nz);
546 MAT3(uS, *nx, *ny, *nz);
547 MAT3(uNE, *nx, *ny, *nz);
548 MAT3(uNW, *nx, *ny, *nz);
549 MAT3(uSE, *nx, *ny, *nz);
550 MAT3(uSW, *nx, *ny, *nz);
551
552 #pragma omp parallel for private(i, j, k, tmpO, tmpU, tmpD)
553 for (k=2; k<=*nz-1; k++) {
554 for (j=2; j<=*ny-1; j++) {
555 for(i=2; i<=*nx-1; i++) {
556
557 tmpO =
558 + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
559 + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
560 + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
561 + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
562 + VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
563 + VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
564 + VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
565 + VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
566
567 tmpU =
568 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
569 + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
570 + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
571 + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
572 + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
573 + VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
574 + VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
575 + VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
576 + VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
577
578 tmpD =
579 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
580 + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
581 + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
582 + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
583 + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
584 + VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
585 + VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
586 + VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
587 + VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
588
589 VAT3(r, i, j, k) = VAT3(fc, i, j, k) + tmpO + tmpU + tmpD
590 - (VAT3(oC, i, j, k) + VAT3(cc, i, j, k)) * VAT3(x, i, j, k);
591 }
592 }
593 }
594}
595
596
597
598VPUBLIC void Vnmresid(int *nx, int *ny, int *nz,
599 int *ipc, double *rpc,
600 double *ac, double *cc, double *fc,
601 double *x, double *r, double *w1) {
602
603 int numdia;
604
605 // Do in oNe step ***
606 numdia = VAT(ipc, 11);
607 if (numdia == 7) {
608 Vnmresid7(nx, ny, nz, ipc, rpc, ac, cc, fc, x, r, w1);
609 } else if (numdia == 27) {
610 Vnmresid27(nx, ny, nz, ipc, rpc, ac, cc, fc, x, r, w1);
611 } else {
612 Vnm_print(2, "Vnmresid: invalid stencil type given...\n");
613 }
614}
615
616
617
618VPUBLIC void Vnmresid7(int *nx, int *ny, int *nz,
619 int *ipc, double *rpc,
620 double *ac, double *cc, double *fc,
621 double *x, double *r, double *w1) {
622
623 MAT2(ac, *nx * *ny * *nz, 1);
624
625 // Do in oNe step
626 Vnmresid7_1s(nx, ny, nz,
627 ipc, rpc,
628 RAT2(ac, 1, 1), cc, fc,
629 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
630 x, r, w1);
631}
632
633VPUBLIC void Vnmresid7_1s(int *nx, int *ny, int *nz,
634 int *ipc, double *rpc,
635 double *oC, double *cc, double *fc,
636 double *oE, double *oN, double *uC,
637 double *x, double *r, double *w1) {
638
639 int i, j, k;
640 int ipkey;
641
642 MAT3(oE, *nx, *ny, *nz);
643 MAT3(oN, *nx, *ny, *nz);
644 MAT3(uC, *nx, *ny, *nz);
645 MAT3(cc, *nx, *ny, *nz);
646 MAT3(fc, *nx, *ny, *nz);
647 MAT3(oC, *nx, *ny, *nz);
648 MAT3( x, *nx, *ny, *nz);
649 MAT3( r, *nx, *ny, *nz);
650 MAT3(w1, *nx, *ny, *nz);
651
652 // First get vector nonlinear term to avoid subroutine calls
653 ipkey = VAT(ipc, 10);
654 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
655
656 // The residual
657 for (k=2; k<=*nz-1; k++) {
658 for (j=2; j<=*ny-1; j++) {
659 for (i=2; i<=*nx-1; i++) {
660 VAT3(r, i, j, k) = VAT3(fc, i, j, k)
661 + VAT3(oN, i, j, k) * VAT3(x, i, j+1, k)
662 + VAT3(oN, i, j-1, k) * VAT3(x, i, j-1, k)
663 + VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
664 + VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
665 + VAT3(uC, i, j, k-1) * VAT3(x, i, j, k-1)
666 + VAT3(uC, i, j, k) * VAT3(x, i, j, k+1)
667 - VAT3(oC, i, j, k) * VAT3(x, i, j, k)
668 - VAT3(w1, i, j, k);
669 }
670 }
671 }
672}
673
674
675
676VPUBLIC void Vnmresid27(int *nx, int *ny, int *nz,
677 int *ipc, double *rpc,
678 double *ac, double *cc, double *fc,
679 double *x, double *r, double *w1) {
680
681 MAT2(ac, *nx * *ny * *nz, 1);
682
683 // Do in oNe step
684 Vnmresid27_1s(nx, ny, nz,
685 ipc, rpc,
686 RAT2(ac, 1, 1), cc, fc,
687 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
688 RAT2(ac, 1, 5), RAT2(ac, 1, 6),
689 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1, 10),
690 RAT2(ac, 1, 11), RAT2(ac, 1, 12), RAT2(ac, 1, 13), RAT2(ac, 1, 14),
691 x, r, w1);
692}
693
694
695
696VPUBLIC void Vnmresid27_1s(int *nx, int *ny, int *nz,
697 int *ipc, double *rpc,
698 double *oC, double *cc, double *fc,
699 double *oE, double *oN, double *uC,
700 double *oNE, double *oNW,
701 double *uE, double *uW, double *uN, double *uS,
702 double *uNE, double *uNW, double *uSE, double *uSW,
703 double *x, double *r, double *w1) {
704
705 int i, j, k;
706 int ipkey;
707 double tmpO, tmpU, tmpD;
708
709 MAT3( oC, *nx, *ny, *nz);
710 MAT3( cc, *nx, *ny, *nz);
711 MAT3( fc, *nx, *ny, *nz);
712 MAT3( oE, *nx, *ny, *nz);
713 MAT3( oN, *nx, *ny, *nz);
714 MAT3( uC, *nx, *ny, *nz);
715 MAT3(oNE, *nx, *ny, *nz);
716 MAT3(oNW, *nx, *ny, *nz);
717 MAT3( uE, *nx, *ny, *nz);
718 MAT3( uW, *nx, *ny, *nz);
719 MAT3( uN, *nx, *ny, *nz);
720 MAT3( uS, *nx, *ny, *nz);
721 MAT3(uNE, *nx, *ny, *nz);
722 MAT3(uNW, *nx, *ny, *nz);
723 MAT3(uSE, *nx, *ny, *nz);
724 MAT3(uSW, *nx, *ny, *nz);
725 MAT3( x, *nx, *ny, *nz);
726 MAT3( r, *nx, *ny, *nz);
727 MAT3( w1, *nx, *ny, *nz);
728
729 // First get vector noNlinear term to avoid subroutine calls
730 ipkey = VAT(ipc, 10);
731 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
732
733 // The residual
734 for (k=2; k<=*nz-1; k++) {
735 for (j=2; j<=*ny-1; j++) {
736 for (i=2; i<=*nx-1; i++) {
737
738 tmpO =
739 + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
740 + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
741 + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
742 + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
743 + VAT3(oNE, i, j, k) * VAT3(x, i+1, j+1, k)
744 + VAT3(oNW, i, j, k) * VAT3(x, i-1, j+1, k)
745 + VAT3(oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
746 + VAT3(oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
747
748 tmpU =
749 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
750 + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
751 + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
752 + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
753 + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
754 + VAT3(uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
755 + VAT3(uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
756 + VAT3(uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
757 + VAT3(uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
758
759 tmpD =
760 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
761 + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
762 + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
763 + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
764 + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
765 + VAT3(uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
766 + VAT3(uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
767 + VAT3(uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
768 + VAT3(uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
769
770 VAT3(r, i, j, k) =
771 + tmpO + tmpU + tmpD
772 + VAT3(fc, i, j, k)
773 - VAT3(oC, i, j, k) * VAT3(x, i, j, k)
774 - VAT3(w1, i, j, k);
775 }
776 }
777 }
778}
779
780
781
782VPUBLIC void Vrestrc(int *nxf, int *nyf, int *nzf,
783 int *nxc, int *nyc, int *nzc,
784 double *xin, double *xout, double *pc) {
785
786 MAT2(pc, *nxc * *nyc * *nzc, 1 );
787
788 Vrestrc2(nxf, nyf, nzf,
789 nxc, nyc, nzc,
790 xin, xout,
791 RAT2(pc, 1, 1), RAT2(pc, 1, 2), RAT2(pc, 1, 3), RAT2(pc, 1, 4), RAT2(pc, 1, 5),
792 RAT2(pc, 1, 6), RAT2(pc, 1, 7), RAT2(pc, 1, 8), RAT2(pc, 1, 9),
793 RAT2(pc, 1,10), RAT2(pc, 1,11), RAT2(pc, 1,12), RAT2(pc, 1,13), RAT2(pc, 1,14),
794 RAT2(pc, 1,15), RAT2(pc, 1,16), RAT2(pc, 1,17), RAT2(pc, 1,18),
795 RAT2(pc, 1,19), RAT2(pc, 1,20), RAT2(pc, 1,21), RAT2(pc, 1,22), RAT2(pc, 1,23),
796 RAT2(pc, 1,24), RAT2(pc, 1,25), RAT2(pc, 1,26), RAT2(pc, 1,27));
797}
798
799
800
801VEXTERNC void Vrestrc2(int *nxf, int *nyf, int *nzf,
802 int *nxc, int *nyc, int *nzc,
803 double *xin, double *xout,
804 double *oPC, double *oPN, double *oPS, double *oPE, double *oPW,
805 double *oPNE, double *oPNW, double *oPSE, double *oPSW,
806 double *uPC, double *uPN, double *uPS, double *uPE, double *uPW,
807 double *uPNE, double *uPNW, double *uPSE, double *uPSW,
808 double *dPC, double *dPN, double *dPS, double *dPE, double *dPW,
809 double *dPNE, double *dPNW, double *dPSE, double *dPSW) {
810
811 int i, j, k;
812 int ii, jj, kk;
813 int idimenshun = 3;
814
815 double tmpO, tmpU, tmpD;
816 double dimfac;
817
818 MAT3(xin, *nxf, *nyf, *nzf);
819 MAT3(xout, *nxc, *nyc, *nzc);
820
821 MAT3(oPC, *nxc, *nyc, *nzc);
822 MAT3(oPN, *nxc, *nyc, *nzc);
823 MAT3(oPS, *nxc, *nyc, *nzc);
824 MAT3(oPE, *nxc, *nyc, *nzc);
825 MAT3(oPW, *nxc, *nyc, *nzc);
826
827 MAT3(oPNE, *nxc, *nyc, *nzc);
828 MAT3(oPNW, *nxc, *nyc, *nzc);
829 MAT3(oPSE, *nxc, *nyc, *nzc);
830 MAT3(oPSW, *nxc, *nyc, *nzc);
831
832 MAT3(uPC, *nxc, *nyc, *nzc);
833 MAT3(uPN, *nxc, *nyc, *nzc);
834 MAT3(uPS, *nxc, *nyc, *nzc);
835 MAT3(uPE, *nxc, *nyc, *nzc);
836 MAT3(uPW, *nxc, *nyc, *nzc);
837
838 MAT3(uPNE, *nxc, *nyc, *nzc);
839 MAT3(uPNW, *nxc, *nyc, *nzc);
840 MAT3(uPSE, *nxc, *nyc, *nzc);
841 MAT3(uPSW, *nxc, *nyc, *nzc);
842
843 MAT3(dPC, *nxc, *nyc, *nzc);
844 MAT3(dPN, *nxc, *nyc, *nzc);
845 MAT3(dPS, *nxc, *nyc, *nzc);
846 MAT3(dPE, *nxc, *nyc, *nzc);
847 MAT3(dPW, *nxc, *nyc, *nzc);
848
849 MAT3(dPNE, *nxc, *nyc, *nzc);
850 MAT3(dPNW, *nxc, *nyc, *nzc);
851 MAT3(dPSE, *nxc, *nyc, *nzc);
852 MAT3(dPSW, *nxc, *nyc, *nzc);
853
854 // Verify correctness of the input boundary points
855 VfboundPMG00(nxf, nyf, nzf, xin);
856
857 dimfac = VPOW(2.0, idimenshun);
858
859 // Handle the interior points as average of 5 finer grid pts ***
860 #pragma omp parallel for private(k, kk, j, jj, i, ii, tmpO, tmpU, tmpD)
861 for (k=2; k<=*nzc-1; k++) {
862 kk = (k - 1) * 2 + 1;
863
864 for (j=2; j<=*nyc-1; j++) {
865 jj = (j - 1) * 2 + 1;
866
867 for (i=2; i<=*nxc-1; i++) {
868 ii = (i - 1) * 2 + 1;
869
870 // Compute the restriction
871 tmpO =
872 + VAT3( oPC, i, j, k) * VAT3(xin, ii, jj, kk)
873 + VAT3( oPN, i, j, k) * VAT3(xin, ii, jj+1, kk)
874 + VAT3( oPS, i, j, k) * VAT3(xin, ii, jj-1, kk)
875 + VAT3( oPE, i, j, k) * VAT3(xin, ii+1, jj, kk)
876 + VAT3( oPW, i, j, k) * VAT3(xin, ii-1, jj, kk)
877 + VAT3(oPNE, i, j, k) * VAT3(xin, ii+1, jj+1, kk)
878 + VAT3(oPNW, i, j, k) * VAT3(xin, ii-1, jj+1, kk)
879 + VAT3(oPSE, i, j, k) * VAT3(xin, ii+1, jj-1, kk)
880 + VAT3(oPSW, i, j, k) * VAT3(xin, ii-1, jj-1, kk);
881
882 tmpU =
883 + VAT3( uPC, i, j, k) * VAT3(xin, ii, jj, kk+1)
884 + VAT3( uPN, i, j, k) * VAT3(xin, ii, jj+1, kk+1)
885 + VAT3( uPS, i, j, k) * VAT3(xin, ii, jj-1, kk+1)
886 + VAT3( uPE, i, j, k) * VAT3(xin, ii+1, jj, kk+1)
887 + VAT3( uPW, i, j, k) * VAT3(xin, ii-1, jj, kk+1)
888 + VAT3(uPNE, i, j, k) * VAT3(xin, ii+1, jj+1, kk+1)
889 + VAT3(uPNW, i, j, k) * VAT3(xin, ii-1, jj+1, kk+1)
890 + VAT3(uPSE, i, j, k) * VAT3(xin, ii+1, jj-1, kk+1)
891 + VAT3(uPSW, i, j, k) * VAT3(xin, ii-1, jj-1, kk+1);
892
893 tmpD =
894 + VAT3( dPC, i, j, k) * VAT3(xin, ii, jj, kk-1)
895 + VAT3( dPN, i, j, k) * VAT3(xin, ii, jj+1, kk-1)
896 + VAT3( dPS, i, j, k) * VAT3(xin, ii, jj-1, kk-1)
897 + VAT3( dPE, i, j, k) * VAT3(xin, ii+1, jj, kk-1)
898 + VAT3( dPW, i, j, k) * VAT3(xin, ii-1, jj, kk-1)
899 + VAT3(dPNE, i, j, k) * VAT3(xin, ii+1, jj+1, kk-1)
900 + VAT3(dPNW, i, j, k) * VAT3(xin, ii-1, jj+1, kk-1)
901 + VAT3(dPSE, i, j, k) * VAT3(xin, ii+1, jj-1, kk-1)
902 + VAT3(dPSW, i, j, k) * VAT3(xin, ii-1, jj-1, kk-1);
903
904 VAT3(xout, i, j, k) = tmpO + tmpU + tmpD;
905 }
906 }
907 }
908
909 // Verify correctness of the output boundary points
910 VfboundPMG00(nxc, nyc, nzc, xout);
911}
912
913
914
915VPUBLIC void VinterpPMG(int *nxc, int *nyc, int *nzc,
916 int *nxf, int *nyf, int *nzf,
917 double *xin, double *xout,
918 double *pc) {
919
920 MAT2(pc, *nxc * *nyc * *nzc, 1);
921
922 VinterpPMG2(nxc, nyc, nzc,
923 nxf, nyf, nzf,
924 xin, xout,
925 RAT2(pc, 1, 1), RAT2(pc, 1, 2), RAT2(pc, 1, 3), RAT2(pc, 1, 4), RAT2(pc, 1, 5),
926 RAT2(pc, 1, 6), RAT2(pc, 1, 7), RAT2(pc, 1, 8), RAT2(pc, 1, 9),
927 RAT2(pc, 1,10), RAT2(pc, 1,11), RAT2(pc, 1,12), RAT2(pc, 1,13), RAT2(pc, 1,14),
928 RAT2(pc, 1,15), RAT2(pc, 1,16), RAT2(pc, 1,17), RAT2(pc, 1,18),
929 RAT2(pc, 1,19), RAT2(pc, 1,20), RAT2(pc, 1,21), RAT2(pc, 1,22), RAT2(pc, 1,23),
930 RAT2(pc, 1,24), RAT2(pc, 1,25), RAT2(pc, 1,26), RAT2(pc, 1,27));
931}
932
933
934
935VPUBLIC void VinterpPMG2(int *nxc, int *nyc, int *nzc,
936 int *nxf, int *nyf, int *nzf,
937 double *xin, double *xout,
938 double *oPC, double *oPN, double *oPS, double *oPE, double *oPW,
939 double *oPNE, double *oPNW, double *oPSE, double *oPSW,
940 double *uPC, double *uPN, double *uPS, double *uPE, double *uPW,
941 double *uPNE, double *uPNW, double *uPSE, double *uPSW,
942 double *dPC, double *dPN, double *dPS, double *dPE, double *dPW,
943 double *dPNE, double *dPNW, double *dPSE, double *dPSW) {
944
945 int i, j, k;
946 int ii, jj, kk;
947
948 MAT3( xin, *nxc, *nyc, *nzc);
949 MAT3(xout, *nxf, *nyf, *nzf);
950
951 MAT3( oPC, *nxc, *nyc, *nzc);
952 MAT3( oPN, *nxc, *nyc, *nzc);
953 MAT3( oPS, *nxc, *nyc, *nzc);
954 MAT3( oPE, *nxc, *nyc, *nzc);
955 MAT3( oPW, *nxc, *nyc, *nzc);
956
957 MAT3(oPNE, *nxc, *nyc, *nzc);
958 MAT3(oPNW, *nxc, *nyc, *nzc);
959 MAT3(oPSE, *nxc, *nyc, *nzc);
960 MAT3(oPSW, *nxc, *nyc, *nzc);
961
962 MAT3( uPC, *nxc, *nyc, *nzc);
963 MAT3( uPN, *nxc, *nyc, *nzc);
964 MAT3( uPS, *nxc, *nyc, *nzc);
965 MAT3( uPE, *nxc, *nyc, *nzc);
966 MAT3( uPW, *nxc, *nyc, *nzc);
967
968 MAT3(uPNE, *nxc, *nyc, *nzc);
969 MAT3(uPNW, *nxc, *nyc, *nzc);
970 MAT3(uPSE, *nxc, *nyc, *nzc);
971 MAT3(uPSW, *nxc, *nyc, *nzc);
972
973 MAT3( dPC, *nxc, *nyc, *nzc);
974 MAT3( dPN, *nxc, *nyc, *nzc);
975 MAT3( dPS, *nxc, *nyc, *nzc);
976 MAT3( dPE, *nxc, *nyc, *nzc);
977 MAT3( dPW, *nxc, *nyc, *nzc);
978
979 MAT3(dPNE, *nxc, *nyc, *nzc);
980 MAT3(dPNW, *nxc, *nyc, *nzc);
981 MAT3(dPSE, *nxc, *nyc, *nzc);
982 MAT3(dPSW, *nxc, *nyc, *nzc);
983
984 /* *********************************************************************
985 * Setup
986 * *********************************************************************/
987
988 // Verify correctness of the input boundary points ***
989 VfboundPMG00(nxc, nyc, nzc, xin);
990
991 // Do it
992 for (k=1; k<=*nzf-2; k+=2) {
993 kk = (k - 1) / 2 + 1;
994
995 for (j=1; j<=*nyf-2; j+=2) {
996 jj = (j - 1) / 2 + 1;
997
998 for (i=1; i<=*nxf-2; i+=2) {
999 ii = (i - 1) / 2 + 1;
1000
1001 /* ******************************************************** *
1002 * Type 1 -- Fine grid points common to a coarse grid point *
1003 * ******************************************************** */
1004
1005 // Copy coinciding points from coarse grid to fine grid
1006 VAT3(xout, i, j, k) = VAT3(xin, ii, jj, kk);
1007
1008 /* ******************************************************** *
1009 * type 2 -- fine grid points common to a coarse grid plane *
1010 * ******************************************************** */
1011
1012 // Fine grid pts common only to y-z planes on coarse grid
1013 // (intermediate pts between 2 grid points on x-row)
1014 VAT3(xout, i+1, j, k) = VAT3(oPE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1015 + VAT3(oPW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk);
1016
1017 // Fine grid pts common only to x-z planes on coarse grid
1018 // (intermediate pts between 2 grid points on a y-row)
1019 VAT3(xout, i, j+1, k) = VAT3(oPN, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1020 + VAT3(oPS, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk);
1021
1022 // Fine grid pts common only to x-y planes on coarse grid
1023 // (intermediate pts between 2 grid points on a z-row)
1024 VAT3(xout, i, j, k+1) = VAT3(uPC, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1025 + VAT3(dPC, ii, jj, kk+1) * VAT3(xin, ii, jj, kk+1);
1026
1027 /* ******************************************************* *
1028 * type 3 -- fine grid points common to a coarse grid line *
1029 * ******************************************************* */
1030
1031 // Fine grid pts common only to z planes on coarse grid
1032 // (intermediate pts between 4 grid pts on the xy-plane
1033
1034 VAT3(xout, i+1, j+1, k) = VAT3(oPNE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1035 + VAT3(oPNW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk)
1036 + VAT3(oPSE, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk)
1037 + VAT3(oPSW, ii+1, jj+1, kk) * VAT3(xin, ii+1, jj+1, kk);
1038
1039 // Fine grid pts common only to y planes on coarse grid
1040 // (intermediate pts between 4 grid pts on the xz-plane
1041 VAT3(xout, i+1, j, k+1) = VAT3(uPE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1042 + VAT3(uPW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk)
1043 + VAT3(dPE, ii, jj, kk+1) * VAT3(xin, ii, jj, kk+1)
1044 + VAT3(dPW, ii+1, jj, kk+1) * VAT3(xin, ii+1, jj, kk+1);
1045
1046 // Fine grid pts common only to x planes on coarse grid
1047 // (intermediate pts between 4 grid pts on the yz-plane***
1048 VAT3(xout, i, j+1, k+1) = VAT3(uPN, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1049 + VAT3(uPS, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk)
1050 + VAT3(dPN, ii, jj,kk+1) * VAT3(xin, ii, jj, kk+1)
1051 + VAT3(dPS, ii, jj+1,kk+1) * VAT3(xin, ii, jj+1, kk+1);
1052
1053 /* **************************************** *
1054 * type 4 -- fine grid points not common to *
1055 * coarse grid pts/lines/planes *
1056 * **************************************** */
1057
1058 // Completely interior points
1059 VAT3(xout, i+1,j+1,k+1) =
1060 + VAT3(uPNE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1061 + VAT3(uPNW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk)
1062 + VAT3(uPSE, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk)
1063 + VAT3(uPSW, ii+1, jj+1, kk) * VAT3(xin, ii+1, jj+1, kk)
1064 + VAT3(dPNE, ii, jj, kk+1) * VAT3(xin, ii, jj, kk+1)
1065 + VAT3(dPNW, ii+1, jj, kk+1) * VAT3(xin, ii+1, jj, kk+1)
1066 + VAT3(dPSE, ii, jj+1, kk+1) * VAT3(xin, ii, jj+1, kk+1)
1067 + VAT3(dPSW, ii+1, jj+1, kk+1) * VAT3(xin, ii+1, jj+1, kk+1);
1068 }
1069 }
1070 }
1071
1072 // Verify correctness of the output boundary points ***
1073 VfboundPMG00(nxf, nyf, nzf, xout);
1074}
1075
1076
1077
1078VPUBLIC void Vextrac(int *nxf, int *nyf, int *nzf,
1079 int *nxc, int *nyc, int *nzc,
1080 double *xin, double *xout) {
1081
1082 int i, j, k;
1083 int ii, jj, kk;
1084
1085 MAT3( xin, *nxf, *nyf, *nzf);
1086 MAT3(xout, *nxc, *nyc, *nzc);
1087
1088 // Verify correctness of the input boundary points
1089 VfboundPMG00(nxf, nyf, nzf, xin);
1090
1091 // Do it
1092 for (k=2; k<=*nzc-1; k++) {
1093 kk = (k - 1) * 2 + 1;
1094
1095 for (j=2; j<=*nyc-1; j++) {
1096 jj = (j - 1) * 2 + 1;
1097
1098 for (i=2; i<=*nxc-1; i++) {
1099 ii = (i - 1) * 2 + 1;
1100
1101 // Compute the restriction
1102 VAT3(xout, i, j, k) = VAT3(xin, ii, jj, kk);
1103 }
1104 }
1105 }
1106
1107 // Verify correctness of the output boundary points
1108 VfboundPMG00(nxc, nyc, nzc, xout);
1109}
VPUBLIC void Vmatvec(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *x, double *y)
Matrix-vector multiplication routines.
Definition matvecd.c:57
VPUBLIC void Vc_vec(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the nonlinearity (vector version)
Definition mypdec.c:121
VPUBLIC void VinterpPMG(int *nxc, int *nyc, int *nzc, int *nxf, int *nyf, int *nzf, double *xin, double *xout, double *pc)
Apply the prolongation operator.
Definition matvecd.c:915
VPUBLIC void VfboundPMG00(int *nx, int *ny, int *nz, double *x)
Initialize a grid function to have a zero boundary value.
Definition mikpckd.c:258
VPUBLIC void Vmresid(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *x, double *r)
Break the matrix data-structure into diagonals and then call the residual routine.
Definition matvecd.c:426
VPUBLIC void Vnmresid(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *x, double *r, double *w1)
Break the matrix data-structure into diagonals and then call the residual routine.
Definition matvecd.c:598
VPUBLIC void Vrestrc(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, double *xin, double *xout, double *pc)
Apply the restriction operator.
Definition matvecd.c:782
VPUBLIC void Vextrac(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, double *xin, double *xout)
Simple injection of a fine grid function into coarse grid.
Definition matvecd.c:1078
VEXTERNC void Vnmatvec(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *x, double *y, double *w1)
Break the matrix data-structure into diagonals and then call the matrix-vector routine.
Definition matvecd.c:232