APBS 3.0.0
Loading...
Searching...
No Matches
buildPd.c
1
55#include "buildPd.h"
56
57VPUBLIC void VbuildP(int *nxf, int *nyf, int *nzf,
58 int *nxc, int *nyc, int *nzc,
59 int *mgprol,
60 int *ipc, double *rpc,
61 double *pc, double *ac,
62 double *xf, double *yf, double *zf) {
63
64 int numdia;
65
66 MAT2(pc, *nxc * *nyc * *nzc, 1);
67 MAT2(ac, *nxf * *nyf * *nzf, 1);
68
69 if (*mgprol == 0) {
70
71 VbuildP_trilin(nxf, nyf, nzf,
72 nxc, nyc, nzc,
73 RAT2(pc, 1, 1),
74 xf, yf, zf);
75
76 } else if (*mgprol == 1) {
77
78 numdia = VAT(ipc, 11);
79
80 if (numdia == 7) {
81 VbuildP_op7(nxf, nyf, nzf,
82 nxc, nyc, nzc,
83 ipc, rpc,
84 RAT2(ac, 1, 1), RAT2(pc, 1, 1));
85 } else if (numdia == 27) {
86 VbuildP_op27(nxf, nyf, nzf,
87 nxc, nyc, nzc,
88 ipc, rpc,
89 RAT2(ac, 1, 1), RAT2(pc, 1, 1));
90 } else {
91 Vnm_print(2,"BUILDP: invalid stencil type given: %d\n", numdia);
92 }
93 }
94}
95
96VPUBLIC void VbuildP_trilin(int *nxf, int *nyf, int *nzf,
97 int *nxc, int *nyc, int *nzc,
98 double *pc,
99 double *xf, double *yf, double *zf) {
100
101 MAT2(pc, *nxc * *nyc * *nzc, 1);
102
103 VbuildPb_trilin(nxf, nyf, nzf,
104 nxc, nyc, nzc,
105 RAT2(pc, 1, 1),RAT2(pc, 1, 2),RAT2(pc, 1, 3),RAT2(pc, 1, 4),RAT2(pc, 1, 5),
106 RAT2(pc, 1, 6),RAT2(pc, 1, 7),RAT2(pc, 1, 8),RAT2(pc, 1, 9),
107 RAT2(pc, 1, 10),RAT2(pc, 1, 11),RAT2(pc, 1, 12),RAT2(pc, 1, 13),RAT2(pc, 1, 14),
108 RAT2(pc, 1, 15),RAT2(pc, 1, 16),RAT2(pc, 1, 17),RAT2(pc, 1, 18),
109 RAT2(pc, 1, 19),RAT2(pc, 1, 20),RAT2(pc, 1, 21),RAT2(pc, 1, 22),RAT2(pc, 1, 23),
110 RAT2(pc, 1, 24),RAT2(pc, 1, 25),RAT2(pc, 1, 26),RAT2(pc, 1, 27),
111 xf, yf, zf);
112}
113
114VEXTERNC void VbuildPb_trilin(int *nxf, int *nyf, int *nzf,
115 int *nxc, int *nyc, int *nzc,
116 double *oPC, double *oPN, double *oPS, double *oPE, double *oPW,
117 double *oPNE, double *oPNW, double *oPSE, double *oPSW,
118 double *uPC, double *uPN, double *uPS, double *uPE, double *uPW,
119 double *uPNE, double *uPNW, double *uPSE, double *uPSW,
120 double *dPC, double *dPN, double *dPS, double *dPE, double *dPW,
121 double *dPNE, double *dPNW, double *dPSE, double *dPSW,
122 double *xf, double *yf, double *zf) {
123
124 int i, j, k;
125
127 double won = 1.0;
128 double half = 1.0 / 2.0;
129 double quarter = 1.0 / 4.0;
130 double eighth = 1.0 / 8.0;
131
132 MAT3( oPC, *nxc, *nyc, *nzc);
133 MAT3( oPN, *nxc, *nyc, *nzc);
134 MAT3( oPS, *nxc, *nyc, *nzc);
135 MAT3( oPE, *nxc, *nyc, *nzc);
136 MAT3( oPW, *nxc, *nyc, *nzc);
137
138 MAT3(oPNE, *nxc, *nyc, *nzc);
139 MAT3(oPNW, *nxc, *nyc, *nzc);
140 MAT3(oPSE, *nxc, *nyc, *nzc);
141 MAT3(oPSW, *nxc, *nyc, *nzc);
142
143 MAT3( uPC, *nxc, *nyc, *nzc);
144 MAT3( uPN, *nxc, *nyc, *nzc);
145 MAT3( uPS, *nxc, *nyc, *nzc);
146 MAT3( uPE, *nxc, *nyc, *nzc);
147 MAT3( uPW, *nxc, *nyc, *nzc);
148
149 MAT3(uPNE, *nxc, *nyc, *nzc);
150 MAT3(uPNW, *nxc, *nyc, *nzc);
151 MAT3(uPSE, *nxc, *nyc, *nzc);
152 MAT3(uPSW, *nxc, *nyc, *nzc);
153
154 MAT3( dPC, *nxc, *nyc, *nzc);
155 MAT3( dPN, *nxc, *nyc, *nzc);
156 MAT3( dPS, *nxc, *nyc, *nzc);
157 MAT3( dPE, *nxc, *nyc, *nzc);
158 MAT3( dPW, *nxc, *nyc, *nzc);
159
160 MAT3(dPNE, *nxc, *nyc, *nzc);
161 MAT3(dPNW, *nxc, *nyc, *nzc);
162 MAT3(dPSE, *nxc, *nyc, *nzc);
163 MAT3(dPSW, *nxc, *nyc, *nzc);
164
165for(k=2; k<=*nzc-1; k++) {
166 for(j=2; j<=*nyc-1; j++) {
167 for(i=2; i<=*nxc-1; i++) {
168
169 VAT3(oPC, i,j,k) = won;
170
171 VAT3(oPN, i,j,k) = half;
172 VAT3(oPS, i,j,k) = half;
173 VAT3(oPE, i,j,k) = half;
174 VAT3(oPW, i,j,k) = half;
175 VAT3(uPC, i,j,k) = half;
176 VAT3(dPC, i,j,k) = half;
177
178 VAT3(oPNE, i,j,k) = quarter;
179 VAT3(oPNW, i,j,k) = quarter;
180 VAT3(oPSE, i,j,k) = quarter;
181 VAT3(oPSW, i,j,k) = quarter;
182 VAT3(dPN, i,j,k) = quarter;
183 VAT3(dPS, i,j,k) = quarter;
184 VAT3(dPE, i,j,k) = quarter;
185 VAT3(dPW, i,j,k) = quarter;
186 VAT3(uPN, i,j,k) = quarter;
187 VAT3(uPS, i,j,k) = quarter;
188 VAT3(uPE, i,j,k) = quarter;
189 VAT3(uPW, i,j,k) = quarter;
190
191 VAT3(dPNE, i,j,k) = eighth;
192 VAT3(dPNW, i,j,k) = eighth;
193 VAT3(dPSE, i,j,k) = eighth;
194 VAT3(dPSW, i,j,k) = eighth;
195 VAT3(uPNE, i,j,k) = eighth;
196 VAT3(uPNW, i,j,k) = eighth;
197 VAT3(uPSE, i,j,k) = eighth;
198 VAT3(uPSW, i,j,k) = eighth;
199 }
200 }
201 }
202}
203
204
205
206VPUBLIC void VbuildP_op7(int *nxf, int *nyf, int *nzf,
207 int *nxc, int *nyc, int *nzc,
208 int *ipc, double *rpc,
209 double *ac, double *pc) {
210
211 MAT2(ac, *nxf * *nyf * *nzf, 1);
212 MAT2(pc, *nxc * *nyc * *nzc, 1);
213
214 WARN_UNTESTED;
215
216 VbuildPb_op7(nxf, nyf, nzf,
217 nxc, nyc, nzc,
218 ipc, rpc,
219 RAT2(ac, 1, 1), RAT2(ac, 1, 2), RAT2(ac, 1, 3),
220 RAT2(ac, 1, 4),
221 RAT2(pc, 1, 1), RAT2(pc, 1, 2), RAT2(pc, 1, 3), RAT2(pc, 1, 4), RAT2(pc, 1, 5),
222 RAT2(pc, 1, 6), RAT2(pc, 1, 7), RAT2(pc, 1, 8), RAT2(pc, 1, 9),
223 RAT2(pc, 1, 10), RAT2(pc, 1, 11), RAT2(pc, 1, 12), RAT2(pc, 1, 13), RAT2(pc, 1, 14),
224 RAT2(pc, 1, 15), RAT2(pc, 1, 16), RAT2(pc, 1, 17), RAT2(pc, 1, 18),
225 RAT2(pc, 1, 19), RAT2(pc, 1, 20), RAT2(pc, 1, 21), RAT2(pc, 1, 22), RAT2(pc, 1, 23),
226 RAT2(pc, 1, 24), RAT2(pc, 1, 25), RAT2(pc, 1, 26), RAT2(pc, 1, 27));
227}
228
229
230
231VPUBLIC void VbuildPb_op7(int *nxf, int *nyf, int *nzf,
232 int *nxc, int *nyc, int *nzc,
233 int *ipc, double *rpc,
234 double *oC, double *oE, double *oN,
235 double *uC,
236 double *oPC, double *oPN, double *oPS, double *oPE, double *oPW,
237 double *oPNE, double *oPNW, double *oPSE, double *oPSW,
238 double *uPC, double *uPN, double *uPS, double *uPE, double *uPW,
239 double *uPNE, double *uPNW, double *uPSE, double *uPSW,
240 double *dPC, double *dPN, double *dPS, double *dPE, double *dPW,
241 double *dPNE, double *dPNW, double *dPSE, double *dPSW) {
242
243 int i, j, k;
244 int ii, jj, kk;
245 int im1, ip1;
246 int im2, ip2;
247 int jm1, jp1;
248 int jm2, jp2;
249 int km1, kp1;
250 int km2, kp2;
251 int iim1, iip1;
252 int jjm1, jjp1;
253 int kkm1, kkp1;
254
255 double won, half, quarter, eighth;
256
257 MAT3( oC, *nxf, *nyf, *nzf);
258 MAT3( oE, *nxf, *nyf, *nzf);
259 MAT3( oN, *nxf, *nyf, *nzf);
260 MAT3( uC, *nxf, *nyf, *nzf);
261 MAT3( oPC, *nxc, *nyc, *nzc);
262 MAT3( oPN, *nxc, *nyc, *nzc);
263 MAT3( oPS, *nxc, *nyc, *nzc);
264 MAT3( oPE, *nxc, *nyc, *nzc);
265 MAT3( oPW, *nxc, *nyc, *nzc);
266 MAT3(oPNE, *nxc, *nyc, *nzc);
267 MAT3(oPNW, *nxc, *nyc, *nzc);
268 MAT3(oPSE, *nxc, *nyc, *nzc);
269 MAT3(oPSW, *nxc, *nyc, *nzc);
270 MAT3( uPC, *nxc, *nyc, *nzc);
271 MAT3( uPN, *nxc, *nyc, *nzc);
272 MAT3( uPS, *nxc, *nyc, *nzc);
273 MAT3( uPE, *nxc, *nyc, *nzc);
274 MAT3( uPW, *nxc, *nyc, *nzc);
275 MAT3(uPNE, *nxc, *nyc, *nzc);
276 MAT3(uPNW, *nxc, *nyc, *nzc);
277 MAT3(uPSE, *nxc, *nyc, *nzc);
278 MAT3(uPSW, *nxc, *nyc, *nzc);
279 MAT3( dPC, *nxc, *nyc, *nzc);
280 MAT3( dPN, *nxc, *nyc, *nzc);
281 MAT3( dPS, *nxc, *nyc, *nzc);
282 MAT3( dPE, *nxc, *nyc, *nzc);
283 MAT3( dPW, *nxc, *nyc, *nzc);
284 MAT3(dPNE, *nxc, *nyc, *nzc);
285 MAT3(dPNW, *nxc, *nyc, *nzc);
286 MAT3(dPSE, *nxc, *nyc, *nzc);
287 MAT3(dPSW, *nxc, *nyc, *nzc);
288
289 WARN_UNTESTED;
290
291 // interpolation stencil ***
292 won = 1.0;
293 half = 1.0 / 2.0;
294 quarter = 1.0 / 4.0;
295 eighth = 1.0 / 8.0;
296
297 //fprintf(data, "%s\n", PRINT_FUNC);
298
299 for (kk = 2; kk < *nzc - 1; kk++) {
300 k = 2 * kk - 1;
301
302 for (jj = 2; jj < *nyc - 1; jj++) {
303 j = 2 * jj - 1;
304
305 for (ii = 2; ii < *nxc - 1; ii++) {
306 i = 2 * ii - 1;
307
308 // index computations ***
309 im1 = i - 1;
310 ip1 = i + 1;
311 im2 = i - 2;
312 ip2 = i + 2;
313 jm1 = j - 1;
314 jp1 = j + 1;
315 jm2 = j - 2;
316 jp2 = j + 2;
317 km1 = k - 1;
318 kp1 = k + 1;
319 km2 = k - 2;
320 kp2 = k + 2;
321 iim1 = ii - 1;
322 iip1 = ii + 1;
323 jjm1 = jj - 1;
324 jjp1 = jj + 1;
325 kkm1 = kk - 1;
326 kkp1 = kk + 1;
327
328 // *************************************************************
329 // *** > oPC;
330 // *************************************************************
331
332 VAT3( oPC, ii, jj, kk) = won;
333
334 //fprintf(data, "%19.12E\n", VAT3(oPC, ii, jj, kk));
335
336 // *************************************************************
337 // *** > oPN;
338 // *************************************************************
339
340 VAT3( oPN, ii, jj, kk) =
341 VAT3( oN, i, j, k) / ( VAT3( oC, i, jp1, k)
342 - VAT3( oE, im1, jp1, k)
343 - VAT3( oE, i, jp1, k)
344 - VAT3( uC, i, jp1, km1)
345 - VAT3( uC, i, jp1, k));
346
347 //fprintf(data, "%19.12E\n", VAT3(oPN, ii, jj, kk));
348
349 // *************************************************************
350 // *** > oPS;
351 // *************************************************************
352
353 VAT3( oPS, ii, jj, kk) =
354 VAT3( oN, i, jm1, k) / ( VAT3( oC, i, jm1, k)
355 - VAT3( oE, im1, jm1, k)
356 - VAT3( oE, i, jm1, k)
357 - VAT3( uC, i, jm1, km1)
358 - VAT3( uC, i, jm1, k));
359
360 //fprintf(data, "%19.12E\n", VAT3(oPS, ii, jj, kk));
361
362 // *************************************************************
363 // *** > oPE;
364 // *************************************************************
365
366 VAT3( oPE, ii, jj, kk) =
367 VAT3( oE, i, j, k) / ( VAT3( oC, ip1, j, k)
368 - VAT3( uC, ip1, j, km1)
369 - VAT3( uC, ip1, j, k)
370 - VAT3( oN, ip1, j, k)
371 - VAT3( oN, ip1, jm1, k));
372
373 //fprintf(data, "%19.12E\n", VAT3(oPE, ii, jj, kk));
374
375 // *************************************************************
376 // *** > oPW;
377 // *************************************************************
378
379 VAT3( oPW, ii, jj, kk) =
380 VAT3( oE, im1, j, k) / ( VAT3( oC, im1, j, k)
381 - VAT3( uC, im1, j, km1)
382 - VAT3( uC, im1, j, k)
383 - VAT3( oN, im1, j, k)
384 - VAT3( oN, im1, jm1, k));
385
386 //fprintf(data, "%19.12E\n", VAT3(oPW, ii, jj, kk));
387
388 // *************************************************************
389 // *** > oPNE;
390 // *************************************************************
391
392 VAT3(oPNE, ii, jj, kk) =
393 (
394 VAT3( oN, ip1, j, k) * VAT3( oPE, ii, jj, kk)
395 + VAT3( oE, i, jp1, k) * VAT3( oPN, ii, jj, kk)
396 ) / (
397 VAT3( oC, ip1, jp1, k)
398 - VAT3( uC, ip1, jp1, km1)
399 - VAT3( uC, ip1, jp1, k)
400 );
401
402 //fprintf(data, "%19.12E\n", VAT3(oPNE, ii, jj, kk));
403
404 // *************************************************************
405 // *** > oPNW;
406 // *************************************************************
407
408 VAT3(oPNW, ii, jj, kk) =
409 (
410 VAT3( oN, im1, j, k) * VAT3( oPW, ii, jj, kk)
411 + VAT3( oE, im1, jp1, k) * VAT3( oPN, ii, jj, kk)
412 ) / (
413 VAT3( oC, im1, jp1, k)
414 - VAT3( uC, im1, jp1, km1)
415 - VAT3( uC, im1, jp1, k)
416 );
417
418 //fprintf(data, "%19.12E\n", VAT3(oPNW, ii, jj, kk));
419
420 // *************************************************************
421 // *** > oPSE;
422 // *************************************************************
423
424 VAT3(oPSE, ii, jj, kk) =
425 (
426 VAT3( oN, ip1, jm1, k) * VAT3( oPE, ii, jj, kk)
427 + VAT3( oE, i, jm1, k) * VAT3( oPS, ii, jj, kk)
428 ) / (
429 VAT3( oC, ip1, jm1, k)
430 - VAT3( uC, ip1, jm1, km1)
431 - VAT3( uC, ip1, jm1, k)
432 );
433
434 //fprintf(data, "%19.12E\n", VAT3(oPSE, ii, jj, kk));
435
436 // *************************************************************
437 // *** > oPSW;
438 // *************************************************************
439
440 VAT3(oPSW, ii, jj, kk) =
441 (
442 VAT3( oN, im1, jm1, k) * VAT3( oPW, ii, jj, kk)
443 + VAT3( oE, im1, jm1, k) * VAT3( oPS, ii, jj, kk)
444 ) / (
445 VAT3( oC, im1, jm1, k)
446 - VAT3( uC, im1, jm1, km1)
447 - VAT3( uC, im1, jm1, k)
448 );
449
450 //fprintf(data, "%19.12E\n", VAT3(oPSW, ii, jj, kk));
451
452 // *************************************************************
453 // *** > dPC;
454 // *************************************************************
455
456 VAT3( dPC, ii, jj, kk) =
457 VAT3( uC, i, j, km1)
458 / (
459 VAT3( oC, i, j, km1)
460 - VAT3( oN, i, j, km1)
461 - VAT3( oN, i, jm1, km1)
462 - VAT3( oE, im1, j, km1)
463 - VAT3( oE, i, j, km1)
464 );
465
466 //fprintf(data, "%19.12E\n", VAT3(dPC, ii, jj, kk));
467
468 // *************************************************************
469 // *** > dPN;
470 // *************************************************************
471
472 VAT3( dPN, ii, jj, kk) =
473 (
474 VAT3( oN, i, j, km1) * VAT3( dPC, ii, jj, kk)
475 + VAT3( uC, i, jp1, km1) * VAT3( oPN, ii, jj, kk)
476 ) / (
477 VAT3( oC, i, jp1, km1)
478 - VAT3( oE, im1, jp1, km1)
479 - VAT3( oE, i, jp1, km1)
480 );
481
482 //fprintf(data, "%19.12E\n", VAT3(dPN, ii, jj, kk));
483
484 // *************************************************************
485 // *** > dPS;
486 // *************************************************************
487
488 VAT3( dPS, ii, jj, kk) =
489 (
490 VAT3( oN, i, jm1, km1) * VAT3( dPC, ii, jj, kk)
491 + VAT3( uC, i, jm1, km1) * VAT3( oPS, ii, jj, kk)
492 ) / (
493 VAT3( oC, i, jm1, km1)
494 - VAT3( oE, im1, jm1, km1)
495 - VAT3( oE, i, jm1, km1)
496 );
497
498 //fprintf(data, "%19.12E\n", VAT3(dPS, ii, jj, kk));
499
500 // *************************************************************
501 // *** > dPE;
502 // *************************************************************
503
504 VAT3( dPE, ii, jj, kk) =
505 (
506 VAT3( uC, ip1, j, km1) * VAT3( oPE, ii, jj, kk)
507 + VAT3( oE, i, j, km1) * VAT3( dPC, ii, jj, kk)
508 ) / (
509 VAT3( oC, ip1, j, km1)
510 - VAT3( oN, ip1, j, km1)
511 - VAT3( oN, ip1, jm1, km1)
512 );
513
514 //fprintf(data, "%19.12E\n", VAT3(dPE, ii, jj, kk));
515
516 // *************************************************************
517 // *** > dPW;
518 // *************************************************************
519
520 VAT3( dPW, ii, jj, kk) =
521 (
522 VAT3( uC, im1, j, km1) * VAT3( oPW, ii, jj, kk)
523 + VAT3( oE, im1, j, km1) * VAT3( dPC, ii, jj, kk)
524 ) / (
525 VAT3( oC, im1, j, km1)
526 - VAT3( oN, im1, j, km1)
527 - VAT3( oN, im1, jm1, km1)
528 );
529
530 //fprintf(data, "%19.12E\n", VAT3(dPW, ii, jj, kk));
531
532 // *************************************************************
533 // *** > dPNE;
534 // *************************************************************
535
536 VAT3(dPNE, ii, jj, kk) =
537 (
538 VAT3( uC, ip1, jp1, km1) * VAT3(oPNE, ii, jj, kk)
539 + VAT3( oE, i, jp1, km1) * VAT3( dPN, ii, jj, kk)
540 + VAT3( oN, ip1, j, km1) * VAT3( dPE, ii, jj, kk)
541 ) / VAT3( oC, ip1, jp1, km1);
542
543 //fprintf(data, "%19.12E\n", VAT3(dPNE, ii, jj, kk));
544
545 // *************************************************************
546 // *** > dPNW;
547 // *************************************************************
548
549 VAT3(dPNW, ii, jj, kk) =
550 (
551 VAT3( uC, im1, jp1, km1) * VAT3(oPNW, ii, jj, kk)
552 + VAT3( oE, im1, jp1, km1) * VAT3( dPN, ii, jj, kk)
553 + VAT3( oN, im1, j, km1) * VAT3( dPW, ii, jj, kk)
554 ) / VAT3( oC, im1, jp1, km1);
555
556 //fprintf(data, "%19.12E\n", VAT3(dPNW, ii, jj, kk));
557
558 // *************************************************************
559 // *** > dPSE;
560 // *************************************************************
561
562 VAT3(dPSE, ii, jj, kk) =
563 (
564 VAT3( uC, ip1, jm1, km1) * VAT3(oPSE, ii, jj, kk)
565 + VAT3( oE, i, jm1, km1) * VAT3( dPS, ii, jj, kk)
566 + VAT3( oN, ip1, jm1, km1) * VAT3( dPE, ii, jj, kk)
567 ) / VAT3( oC, ip1, jm1, km1);
568
569 //fprintf(data, "%19.12E\n", VAT3(dPSE, ii, jj, kk));
570
571 // *************************************************************
572 // *** > dPSW;
573 // *************************************************************
574
575 VAT3(dPSW, ii, jj, kk) =
576 (
577 VAT3( uC, im1, jm1, km1) * VAT3(oPSW, ii, jj, kk)
578 + VAT3( oE, im1, jm1, km1) * VAT3( dPS, ii, jj, kk)
579 + VAT3( oN, im1, jm1, km1) * VAT3( dPW, ii, jj, kk)
580 ) / VAT3( oC, im1, jm1, km1);
581
582 //fprintf(data, "%19.12E\n", VAT3(dPSW, ii, jj, kk));
583
584 // *************************************************************
585 // *** > uPC;
586 // *************************************************************
587
588 VAT3( uPC, ii, jj, kk) =
589 VAT3( uC, i, j, k)
590 / ( VAT3( oC, i, j, kp1)
591 - VAT3( oN, i, j, kp1)
592 - VAT3( oN, i, jm1, kp1)
593 - VAT3( oE, im1, j, kp1)
594 - VAT3( oE, i, j, kp1)
595 );
596
597 //fprintf(data, "%19.12E\n", VAT3(uPC, ii, jj, kk));
598
599 // *************************************************************
600 // *** > uPN;
601 // *************************************************************
602
603 VAT3( uPN, ii, jj, kk) =
604 (
605 VAT3( oN, i, j, kp1) * VAT3( uPC, ii, jj, kk)
606 + VAT3( uC, i, jp1, k) * VAT3( oPN, ii, jj, kk)
607 ) / (
608 VAT3( oC, i, jp1, kp1)
609 - VAT3( oE, im1, jp1, kp1)
610 - VAT3( oE, i, jp1, kp1)
611 );
612
613 //fprintf(data, "%19.12E\n", VAT3(uPN, ii, jj, kk));
614
615 // *************************************************************
616 // *** > uPS;
617 // *************************************************************
618
619 VAT3( uPS, ii, jj, kk) =
620 (
621 VAT3( oN, i, jm1, kp1) * VAT3( uPC, ii, jj, kk)
622 + VAT3( uC, i, jm1, k) * VAT3( oPS, ii, jj, kk)
623 ) / (
624 VAT3( oC, i, jm1, kp1)
625 - VAT3( oE, im1, jm1, kp1)
626 - VAT3( oE, i, jm1, kp1)
627 );
628
629 //fprintf(data, "%19.12E\n", VAT3(uPS, ii, jj, kk));
630
631 // *************************************************************
632 // *** > uPE;
633 // *************************************************************
634
635 VAT3( uPE, ii, jj, kk) =
636 (
637 VAT3( uC, ip1, j, k) * VAT3( oPE, ii, jj, kk)
638 + VAT3( oE, i, j, kp1) * VAT3( uPC, ii, jj, kk)
639 ) / (
640 VAT3( oC, ip1, j, kp1)
641 - VAT3( oN, ip1, j, kp1)
642 - VAT3( oN, ip1, jm1, kp1)
643 );
644
645 //fprintf(data, "%19.12E\n", VAT3(uPE, ii, jj, kk));
646
647 // *************************************************************
648 // *** > uPW;
649 // *************************************************************
650
651 VAT3( uPW, ii, jj, kk) =
652 (
653 VAT3( uC, im1, j, k) * VAT3( oPW, ii, jj, kk)
654 + VAT3( oE, im1, j, kp1) * VAT3( uPC, ii, jj, kk)
655 ) / (
656 VAT3( oC, im1, j, kp1)
657 - VAT3( oN, im1, j, kp1)
658 - VAT3( oN, im1, jm1, kp1)
659 );
660
661 //fprintf(data, "%19.12E\n", VAT3(uPW, ii, jj, kk));
662
663 // *************************************************************
664 // *** > uPNE;
665 // *************************************************************
666
667 VAT3(uPNE, ii, jj, kk) =
668 (
669 VAT3( uC, ip1, jp1, k) * VAT3(oPNE, ii, jj, kk)
670 + VAT3( oE, i, jp1, kp1) * VAT3( uPN, ii, jj, kk)
671 + VAT3( oN, ip1, j, kp1) * VAT3( uPE, ii, jj, kk)
672 ) / VAT3( oC, ip1, jp1, kp1);
673
674 //fprintf(data, "%19.12E\n", VAT3(uPNE, ii, jj, kk));
675
676 // *************************************************************
677 // *** > uPNW;
678 // *************************************************************
679
680 VAT3(uPNW, ii, jj, kk) =
681 (
682 VAT3( uC, im1, jp1, k) * VAT3(oPNW, ii, jj, kk)
683 + VAT3( oE, im1, jp1, kp1) * VAT3( uPN, ii, jj, kk)
684 + VAT3( oN, im1, j, kp1) * VAT3( uPW, ii, jj, kk)
685 ) / VAT3( oC, im1, jp1, kp1);
686
687 //fprintf(data, "%19.12E\n", VAT3(uPNW, ii, jj, kk));
688
689 // *************************************************************
690 // *** > uPSE;
691 // *************************************************************
692
693 VAT3(uPSE, ii, jj, kk) =
694 (
695 VAT3( uC, ip1, jm1, k) * VAT3(oPSE, ii, jj, kk)
696 + VAT3( oE, i, jm1, kp1) * VAT3( uPS, ii, jj, kk)
697 + VAT3( oN, ip1, jm1, kp1) * VAT3( uPE, ii, jj, kk)
698 ) / VAT3( oC, ip1, jm1, kp1);
699
700 //fprintf(data, "%19.12E\n", VAT3(uPSE, ii, jj, kk));
701
702 // *************************************************************
703 // *** > uPSW;
704 // *************************************************************
705
706 VAT3(uPSW, ii, jj, kk) =
707 (
708 VAT3( uC, im1, jm1, k) * VAT3(oPSW, ii, jj, kk)
709 + VAT3( oE, im1, jm1, kp1) * VAT3( uPS, ii, jj, kk)
710 + VAT3( oN, im1, jm1, kp1) * VAT3( uPW, ii, jj, kk)
711 ) / VAT3( oC, im1, jm1, kp1);
712
713 //fprintf(data, "%19.12E\n", VAT3(uPSW, ii, jj, kk));
714
715 }
716 }
717 }
718}
719
720
721
722VPUBLIC void VbuildP_op27(int *nxf, int *nyf, int *nzf,
723 int *nxc, int *nyc, int *nzc,
724 int *ipc, double *rpc,
725 double *ac, double *pc) {
726
727 MAT2(ac, *nxf * *nyf * *nzf, 1);
728 MAT2(pc, *nxc * *nyc * *nzc, 1);
729
730 WARN_UNTESTED;
731
732 VbuildPb_op27(nxf, nyf, nzf,
733 nxc, nyc, nzc,
734 ipc, rpc,
735
736 RAT2(ac, 1, 1), RAT2(ac, 1, 2), RAT2(ac, 1, 3),
737 RAT2(ac, 1, 4),
738 RAT2(ac, 1, 5), RAT2(ac, 1, 6),
739 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1, 10),
740 RAT2(ac, 1, 11), RAT2(ac, 1, 12), RAT2(ac, 1, 13), RAT2(ac, 1, 14),
741 RAT2(pc, 1, 1), RAT2(pc, 1, 2), RAT2(pc, 1, 3), RAT2(pc, 1, 4), RAT2(pc, 1, 5),
742 RAT2(pc, 1, 6), RAT2(pc, 1, 7), RAT2(pc, 1, 8), RAT2(pc, 1, 9),
743 RAT2(pc, 1, 10), RAT2(pc, 1, 11), RAT2(pc, 1, 12), RAT2(pc, 1, 13), RAT2(pc, 1, 14),
744 RAT2(pc, 1, 15), RAT2(pc, 1, 16), RAT2(pc, 1, 17), RAT2(pc, 1, 18),
745 RAT2(pc, 1, 19), RAT2(pc, 1, 20), RAT2(pc, 1, 21), RAT2(pc, 1, 22), RAT2(pc, 1, 23),
746 RAT2(pc, 1, 24), RAT2(pc, 1, 25), RAT2(pc, 1, 26), RAT2(pc, 1, 27));
747}
748
749VPUBLIC void VbuildPb_op27(int *nxf, int *nyf, int *nzf,
750 int *nxc, int *nyc, int *nzc,
751 int *ipc, double *rpc,
752 double *oC, double *oE, double *oN,
753 double *uC,
754 double *oNE, double *oNW,
755 double *uE, double *uW, double *uN, double *uS,
756 double *uNE, double *uNW, double *uSE, double *uSW,
757 double *oPC, double *oPN, double *oPS, double *oPE, double *oPW,
758 double *oPNE, double *oPNW, double *oPSE, double *oPSW,
759 double *uPC, double *uPN, double *uPS, double *uPE, double *uPW,
760 double *uPNE, double *uPNW, double *uPSE, double *uPSW,
761 double *dPC, double *dPN, double *dPS, double *dPE, double *dPW,
762 double *dPNE, double *dPNW, double *dPSE, double *dPSW) {
763
764 int i, j, k;
765 int ii, jj, kk;
766 int im1, ip1;
767 int im2, ip2;
768 int jm1, jp1;
769 int jm2, jp2;
770 int km1, kp1;
771 int km2, kp2;
772 int iim1, iip1;
773 int jjm1, jjp1;
774 int kkm1, kkp1;
775
776 double won, half, quarter, eighth;
777
778 MAT3( oC, *nxf, *nyf, *nzf);
779 MAT3( oE, *nxf, *nyf, *nzf);
780 MAT3( oN, *nxf, *nyf, *nzf);
781 MAT3( uC, *nxf, *nyf, *nzf);
782 MAT3( oNE, *nxf, *nyf, *nzf);
783 MAT3( oNW, *nxf, *nyf, *nzf);
784 MAT3( uE, *nxf, *nyf, *nzf);
785 MAT3( uW, *nxf, *nyf, *nzf);
786 MAT3( uN, *nxf, *nyf, *nzf);
787 MAT3( uS, *nxf, *nyf, *nzf);
788 MAT3( uNE, *nxf, *nyf, *nzf);
789 MAT3( uNW, *nxf, *nyf, *nzf);
790 MAT3( uSE, *nxf, *nyf, *nzf);
791 MAT3( uSW, *nxf, *nyf, *nzf);
792 MAT3( oPC, *nxc, *nyc, *nzc);
793 MAT3( oPN, *nxc, *nyc, *nzc);
794 MAT3( oPS, *nxc, *nyc, *nzc);
795 MAT3( oPE, *nxc, *nyc, *nzc);
796 MAT3( oPW, *nxc, *nyc, *nzc);
797 MAT3(oPNE, *nxc, *nyc, *nzc);
798 MAT3(oPNW, *nxc, *nyc, *nzc);
799 MAT3(oPSE, *nxc, *nyc, *nzc);
800 MAT3(oPSW, *nxc, *nyc, *nzc);
801 MAT3( uPC, *nxc, *nyc, *nzc);
802 MAT3( uPN, *nxc, *nyc, *nzc);
803 MAT3( uPS, *nxc, *nyc, *nzc);
804 MAT3( uPE, *nxc, *nyc, *nzc);
805 MAT3( uPW, *nxc, *nyc, *nzc);
806 MAT3(uPNE, *nxc, *nyc, *nzc);
807 MAT3(uPNW, *nxc, *nyc, *nzc);
808 MAT3(uPSE, *nxc, *nyc, *nzc);
809 MAT3(uPSW, *nxc, *nyc, *nzc);
810 MAT3( dPC, *nxc, *nyc, *nzc);
811 MAT3( dPN, *nxc, *nyc, *nzc);
812 MAT3( dPS, *nxc, *nyc, *nzc);
813 MAT3( dPE, *nxc, *nyc, *nzc);
814 MAT3( dPW, *nxc, *nyc, *nzc);
815 MAT3(dPNE, *nxc, *nyc, *nzc);
816 MAT3(dPNW, *nxc, *nyc, *nzc);
817 MAT3(dPSE, *nxc, *nyc, *nzc);
818 MAT3(dPSW, *nxc, *nyc, *nzc);
819
820 WARN_UNTESTED;
821
822 // Interpolation Stencil
823 won = 1.0;
824 half = 1.0 / 2.0;
825 quarter = 1.0 / 4.0;
826 eighth = 1.0 / 8.0;
827
828 //fprintf(data, "%s\n", PRINT_FUNC);
829
830 for (kk = 2; kk <= *nzc - 1; kk++) {
831 k = 2 * kk - 1;
832
833 for (jj = 2; jj <= *nyc - 1; jj++) {
834 j = 2 * jj - 1;
835
836 for (ii = 2; ii <= *nxc - 1; ii++) {
837 i = 2 * ii - 1;
838
839 // Index Computations
840 im1 = i - 1;
841 ip1 = i + 1;
842 im2 = i - 2;
843 ip2 = i + 2;
844 jm1 = j - 1;
845 jp1 = j + 1;
846 jm2 = j - 2;
847 jp2 = j + 2;
848 km1 = k - 1;
849 kp1 = k + 1;
850 km2 = k - 2;
851 kp2 = k + 2;
852 iim1 = ii - 1;
853 iip1 = ii + 1;
854 jjm1 = jj - 1;
855 jjp1 = jj + 1;
856 kkm1 = kk - 1;
857 kkp1 = kk + 1;
858
859 //* **********************************************************
860 //* *** > oPC;
861 //* **********************************************************
862
863 VAT3( oPC, ii, jj, kk) = won;
864
865 //fprintf(data, "%19.12E\n", VAT3(oPC, ii, jj, kk));
866
867 //* **********************************************************
868 //* *** > oPN;
869 //* **********************************************************
870
871 VAT3( oPN, ii, jj, kk) =
872 (
873 VAT3( uNE, im1, j, km1)
874 + VAT3( uN, i, j, km1)
875 + VAT3( uNW, ip1, j, km1)
876 + VAT3( oNE, im1, j, k)
877 + VAT3( oN, i, j, k)
878 + VAT3( oNW, ip1, j, k)
879 + VAT3( uSW, i, jp1, k)
880 + VAT3( uS, i, jp1, k)
881 + VAT3( uSE, i, jp1, k)
882 ) / (
883 VAT3( oC, i, jp1, k)
884 - VAT3( oE, im1, jp1, k)
885 - VAT3( oE, i, jp1, k)
886 - VAT3( uC, i, jp1, km1)
887 - VAT3( uE, im1, jp1, km1)
888 - VAT3( uW, ip1, jp1, km1)
889 - VAT3( uC, i, jp1, k)
890 - VAT3( uW, i, jp1, k)
891 - VAT3( uE, i, jp1, k)
892 );
893
894 //fprintf(data, "%19.12E\n", VAT3(oPN, ii, jj, kk));
895
896 //* **********************************************************
897 //* *** > oPS;
898 //* **********************************************************
899
900 VAT3( oPS, ii, jj, kk) =
901 (
902 VAT3( uSE, im1, j, km1)
903 + VAT3( uS, i, j, km1)
904 + VAT3( uSW, ip1, j, km1)
905 + VAT3( oNW, i, jm1, k)
906 + VAT3( oN, i, jm1, k)
907 + VAT3( oNE, i, jm1, k)
908 + VAT3( uNW, i, jm1, k)
909 + VAT3( uN, i, jm1, k)
910 + VAT3( uNE, i, jm1, k)
911 ) / (
912 VAT3( oC, i, jm1, k)
913 - VAT3( oE, im1, jm1, k)
914 - VAT3( oE, i, jm1, k)
915 - VAT3( uC, i, jm1, km1)
916 - VAT3( uE, im1, jm1, km1)
917 - VAT3( uW, ip1, jm1, km1)
918 - VAT3( uC, i, jm1, k)
919 - VAT3( uW, i, jm1, k)
920 - VAT3( uE, i, jm1, k)
921 );
922
923 //fprintf(data, "%19.12E\n", VAT3(oPS, ii, jj, kk));
924
925 //* **********************************************************
926 //* *** > oPE;
927 //* **********************************************************
928
929 VAT3( oPE, ii, jj, kk) =
930 (
931 VAT3( uSE, i, jp1, km1)
932 + VAT3( oNW, ip1, j, k)
933 + VAT3( uNW, ip1, j, k)
934 + VAT3( uE, i, j, km1)
935 + VAT3( oE, i, j, k)
936 + VAT3( uW, ip1, j, k)
937 + VAT3( uNE, i, jm1, km1)
938 + VAT3( oNE, i, jm1, k)
939 + VAT3( uSW, ip1, j, k)
940 ) / (
941 VAT3( oC, ip1, j, k)
942 - VAT3( uC, ip1, j, km1)
943 - VAT3( uC, ip1, j, k)
944 - VAT3( oN, ip1, j, k)
945 - VAT3( uS, ip1, jp1, km1)
946 - VAT3( uN, ip1, j, k)
947 - VAT3( oN, ip1, jm1, k)
948 - VAT3( uN, ip1, jm1, km1)
949 - VAT3( uS, ip1, j, k)
950 );
951
952 //fprintf(data, "%19.12E\n", VAT3(oPE, ii, jj, kk));
953
954 //* **********************************************************
955 //* *** > oPW;
956 //* **********************************************************
957
958 VAT3( oPW, ii, jj, kk) =
959 (
960 VAT3( uSW, i, jp1, km1)
961 + VAT3( oNE, im1, j, k)
962 + VAT3( uNE, im1, j, k)
963 + VAT3( uW, i, j, km1)
964 + VAT3( oE, im1, j, k)
965 + VAT3( uE, im1, j, k)
966 + VAT3( uNW, i, jm1, km1)
967 + VAT3( oNW, i, jm1, k)
968 + VAT3( uSE, im1, j, k)
969 ) / (
970 VAT3( oC, im1, j, k)
971 - VAT3( uC, im1, j, km1)
972 - VAT3( uC, im1, j, k)
973 - VAT3( oN, im1, j, k)
974 - VAT3( uS, im1, jp1, km1)
975 - VAT3( uN, im1, j, k)
976 - VAT3( oN, im1, jm1, k)
977 - VAT3( uN, im1, jm1, km1)
978 - VAT3( uS, im1, j, k)
979 );
980
981 //fprintf(data, "%19.12E\n", VAT3(oPW, ii, jj, kk));
982
983 //* **********************************************************
984 //* *** > oPNE;
985 //* **********************************************************
986
987 VAT3(oPNE, ii, jj, kk) =
988 (
989 VAT3( uNE, i, j, km1)
990 + VAT3( oNE, i, j, k)
991 + VAT3( uSW, ip1, jp1, k)
992 + (
993 VAT3( uN, ip1, j, km1)
994 + VAT3( oN, ip1, j, k)
995 + VAT3( uS, ip1, jp1, k)
996 )
997 * VAT3( oPE, ii, jj, kk)
998 + (
999 VAT3( uE, i, jp1, km1)
1000 + VAT3( oE, i, jp1, k)
1001 + VAT3( uW, ip1, jp1, k)
1002 )
1003 * VAT3( oPN, ii, jj, kk)
1004 ) / (
1005 VAT3( oC, ip1, jp1, k)
1006 - VAT3( uC, ip1, jp1, km1)
1007 - VAT3( uC, ip1, jp1, k)
1008 );
1009
1010 //fprintf(data, "%19.12E\n", VAT3(oPNE, ii, jj, kk));
1011
1012 //* **********************************************************
1013 //* *** > oPNW;
1014 //* **********************************************************
1015
1016 VAT3(oPNW, ii, jj, kk) =
1017 (
1018 VAT3( uNW, i, j, km1)
1019 + VAT3( oNW, i, j, k)
1020 + VAT3( uSE, im1, jp1, k)
1021 + (
1022 VAT3( uN, im1, j, km1)
1023 + VAT3( oN, im1, j, k)
1024 + VAT3( uS, im1, jp1, k)
1025 )
1026 * VAT3( oPW, ii, jj, kk)
1027 + (
1028 VAT3( uW, i, jp1, km1)
1029 + VAT3( oE, im1, jp1, k)
1030 + VAT3( uE, im1, jp1, k)
1031 )
1032 * VAT3( oPN, ii, jj, kk)
1033 ) / (
1034 VAT3( oC, im1, jp1, k)
1035 - VAT3( uC, im1, jp1, km1)
1036 - VAT3( uC, im1, jp1, k)
1037 );
1038
1039 //fprintf(data, "%19.12E\n", VAT3(oPNW, ii, jj, kk));
1040
1041 //* **********************************************************
1042 //* *** > oPSE;
1043 //* **********************************************************
1044
1045 VAT3(oPSE, ii, jj, kk) =
1046 (
1047 VAT3( uSE, i, j, km1)
1048 + VAT3( oNW, ip1, jm1, k)
1049 + VAT3( uNW, ip1, jm1, k)
1050 + (
1051 VAT3( uS, ip1, j, km1)
1052 + VAT3( oN, ip1, jm1, k)
1053 + VAT3( uN, ip1, jm1, k)
1054 )
1055 * VAT3( oPE, ii, jj, kk)
1056 + (
1057 VAT3( uE, i, jm1, km1)
1058 + VAT3( oE, i, jm1, k)
1059 + VAT3( uW, ip1, jm1, k)
1060 )
1061 * VAT3( oPS, ii, jj, kk)
1062 ) / (
1063 VAT3( oC, ip1, jm1, k)
1064 - VAT3( uC, ip1, jm1, km1)
1065 - VAT3( uC, ip1, jm1, k)
1066 );
1067
1068 //fprintf(data, "%19.12E\n", VAT3(oPSE, ii, jj, kk));
1069
1070 //* **********************************************************
1071 //* *** > oPSW;
1072 //* **********************************************************
1073
1074 VAT3(oPSW, ii, jj, kk) =
1075 (
1076 VAT3( uSW, i, j, km1)
1077 + VAT3( oNE, im1, jm1, k)
1078 + VAT3( uNE, im1, jm1, k)
1079 + (
1080 VAT3( uS, im1, j, km1)
1081 + VAT3( oN, im1, jm1, k)
1082 + VAT3( uN, im1, jm1, k)
1083 )
1084 * VAT3( oPW, ii, jj, kk)
1085 + (
1086 VAT3( uW, i, jm1, km1)
1087 + VAT3( oE, im1, jm1, k)
1088 + VAT3( uE, im1, jm1, k)
1089 )
1090 * VAT3( oPS, ii, jj, kk)
1091 ) / (
1092 VAT3( oC, im1, jm1, k)
1093 - VAT3( uC, im1, jm1, km1)
1094 - VAT3( uC, im1, jm1, k)
1095 );
1096
1097 //fprintf(data, "%19.12E\n", VAT3(oPSW, ii, jj, kk));
1098
1099 //* **********************************************************
1100 //* *** > dPC;
1101 //* **********************************************************
1102
1103 VAT3( dPC, ii, jj, kk) =
1104 (
1105 VAT3( uNW, i, j, km1)
1106 + VAT3( uW, i, j, km1)
1107 + VAT3( uSW, i, j, km1)
1108 + VAT3( uN, i, j, km1)
1109 + VAT3( uC, i, j, km1)
1110 + VAT3( uS, i, j, km1)
1111 + VAT3( uNE, i, j, km1)
1112 + VAT3( uE, i, j, km1)
1113 + VAT3( uSE, i, j, km1)
1114 ) / (
1115 VAT3( oC, i, j, km1)
1116 - VAT3( oN, i, j, km1)
1117 - VAT3( oN, i, jm1, km1)
1118 - VAT3( oNW, i, j, km1)
1119 - VAT3( oE, im1, j, km1)
1120 - VAT3( oNE, im1, jm1, km1)
1121 - VAT3( oNE, i, j, km1)
1122 - VAT3( oE, i, j, km1)
1123 - VAT3( oNW, ip1, jm1, km1)
1124 );
1125
1126 //fprintf(data, "%19.12E\n", VAT3(dPC, ii, jj, kk));
1127
1128 //* **********************************************************
1129 //* *** > dPN;
1130 //* **********************************************************
1131
1132 VAT3( dPN, ii, jj, kk) =
1133 (
1134 VAT3( uSW, i, jp1, km1)
1135 + VAT3( uS, i, jp1, km1)
1136 + VAT3( uSE, i, jp1, km1)
1137 + (
1138 VAT3( oNE, im1, j, km1)
1139 + VAT3( oN, i, j, km1)
1140 + VAT3( oNW, ip1, j, km1)
1141 )
1142 * VAT3( dPC, ii, jj, kk)
1143 + (
1144 VAT3( uW, i, jp1, km1)
1145 + VAT3( uC, i, jp1, km1)
1146 + VAT3( uE, i, jp1, km1)
1147 )
1148 * VAT3( oPN, ii, jj, kk)
1149 ) / (
1150 VAT3( oC, i, jp1, km1)
1151 - VAT3( oE, im1, jp1, km1)
1152 - VAT3( oE, i, jp1, km1)
1153 );
1154
1155 //fprintf(data, "%19.12E\n", VAT3(dPN, ii, jj, kk));
1156
1157 //* **********************************************************
1158 //* *** > dPS;
1159 //* **********************************************************
1160
1161 VAT3( dPS, ii, jj, kk) =
1162 (
1163 VAT3( uNW, i, jm1, km1)
1164 + VAT3( uN, i, jm1, km1)
1165 + VAT3( uNE, i, jm1, km1)
1166 + (
1167 VAT3( oNW, i, jm1, km1)
1168 + VAT3( oN, i, jm1, km1)
1169 + VAT3( oNE, i, jm1, km1)
1170 )
1171 * VAT3( dPC, ii, jj, kk)
1172 + (
1173 VAT3( uW, i, jm1, km1)
1174 + VAT3( uC, i, jm1, km1)
1175 + VAT3( uE, i, jm1, km1)
1176 )
1177 * VAT3( oPS, ii, jj, kk)
1178 ) / (
1179 VAT3( oC, i, jm1, km1)
1180 - VAT3( oE, im1, jm1, km1)
1181 - VAT3( oE, i, jm1, km1)
1182 );
1183
1184 //fprintf(data, "%19.12E\n", VAT3(dPS, ii, jj, kk));
1185
1186 //* **********************************************************
1187 //* *** > dPE;
1188 //* **********************************************************
1189
1190 VAT3( dPE, ii, jj, kk) =
1191 (
1192 VAT3( uNW, ip1, j, km1)
1193 + VAT3( uW, ip1, j, km1)
1194 + VAT3( uSW, ip1, j, km1)
1195 + (
1196 VAT3( uN, ip1, j, km1)
1197 + VAT3( uC, ip1, j, km1)
1198 + VAT3( uS, ip1, j, km1)
1199 )
1200 * VAT3( oPE, ii, jj, kk)
1201 + (
1202 VAT3( oNW, ip1, j, km1)
1203 + VAT3( oE, i, j, km1)
1204 + VAT3( oNE, i, jm1, km1)
1205 )
1206 * VAT3( dPC, ii, jj, kk)
1207 ) / (
1208 VAT3( oC, ip1, j, km1)
1209 - VAT3( oN, ip1, j, km1)
1210 - VAT3( oN, ip1, jm1, km1)
1211 );
1212
1213 //fprintf(data, "%19.12E\n", VAT3(dPE, ii, jj, kk));
1214
1215 //* **********************************************************
1216 //* *** > dPW;
1217 //* **********************************************************
1218
1219 VAT3( dPW, ii, jj, kk) =
1220 (
1221 VAT3( uNE, im1, j, km1)
1222 + VAT3( uE, im1, j, km1)
1223 + VAT3( uSE, im1, j, km1)
1224 + (
1225 VAT3( uN, im1, j, km1)
1226 + VAT3( uC, im1, j, km1)
1227 + VAT3( uS, im1, j, km1)
1228 )
1229 * VAT3( oPW, ii, jj, kk)
1230 + (
1231 VAT3( oNE, im1, j, km1)
1232 + VAT3( oE, im1, j, km1)
1233 + VAT3( oNW, i, jm1, km1)
1234 )
1235 * VAT3( dPC, ii, jj, kk)
1236 ) / (
1237 VAT3( oC, im1, j, km1)
1238 - VAT3( oN, im1, j, km1)
1239 - VAT3( oN, im1, jm1, km1)
1240 );
1241
1242 //fprintf(data, "%19.12E\n", VAT3(dPW, ii, jj, kk));
1243
1244 //* **********************************************************
1245 //* *** > dPNE;
1246 //* **********************************************************
1247
1248 VAT3(dPNE, ii, jj, kk) =
1249 (
1250 VAT3( uSW, ip1, jp1, km1)
1251 + VAT3( uW, ip1, jp1, km1)
1252 * VAT3( oPN, ii, jj, kk)
1253 + VAT3( uS, ip1, jp1, km1)
1254 * VAT3( oPE, ii, jj, kk)
1255 + VAT3( uC, ip1, jp1, km1)
1256 * VAT3(oPNE, ii, jj, kk)
1257 + VAT3( oNE, i, j, km1)
1258 * VAT3( dPC, ii, jj, kk)
1259 + VAT3( oE, i, jp1, km1)
1260 * VAT3( dPN, ii, jj, kk)
1261 + VAT3( oN, ip1, j, km1)
1262 * VAT3( dPE, ii, jj, kk)
1263 )
1264 / VAT3( oC, ip1, jp1, km1);
1265
1266 //fprintf(data, "%19.12E\n", VAT3(dPNE, ii, jj, kk));
1267
1268 //* **********************************************************
1269 //* *** > dPNW;
1270 //* **********************************************************
1271
1272 VAT3(dPNW, ii, jj, kk) =
1273 (
1274 VAT3( uSE, im1, jp1, km1)
1275 + VAT3( uE, im1, jp1, km1)
1276 * VAT3( oPN, ii, jj, kk)
1277 + VAT3( uS, im1, jp1, km1)
1278 * VAT3( oPW, ii, jj, kk)
1279 + VAT3( uC, im1, jp1, km1)
1280 * VAT3(oPNW, ii, jj, kk)
1281 + VAT3( oNW, i, j, km1)
1282 * VAT3( dPC, ii, jj, kk)
1283 + VAT3( oE, im1, jp1, km1)
1284 * VAT3( dPN, ii, jj, kk)
1285 + VAT3( oN, im1, j, km1)
1286 * VAT3( dPW, ii, jj, kk)
1287 )
1288 / VAT3( oC, im1, jp1, km1);
1289
1290 //fprintf(data, "%19.12E\n", VAT3(dPNW, ii, jj, kk));
1291
1292 //* **********************************************************
1293 //* *** > dPSE;
1294 //* **********************************************************
1295
1296 VAT3(dPSE, ii, jj, kk) =
1297 (
1298 VAT3( uNW, ip1, jm1, km1)
1299 + VAT3( uW, ip1, jm1, km1)
1300 * VAT3( oPS, ii, jj, kk)
1301 + VAT3( uN, ip1, jm1, km1)
1302 * VAT3( oPE, ii, jj, kk)
1303 + VAT3( uC, ip1, jm1, km1)
1304 * VAT3(oPSE, ii, jj, kk)
1305 + VAT3( oNW, ip1, jm1, km1)
1306 * VAT3( dPC, ii, jj, kk)
1307 + VAT3( oE, i, jm1, km1)
1308 * VAT3( dPS, ii, jj, kk)
1309 + VAT3( oN, ip1, jm1, km1)
1310 * VAT3( dPE, ii, jj, kk)
1311 )
1312 / VAT3( oC, ip1, jm1, km1);
1313
1314 //fprintf(data, "%19.12E\n", VAT3(dPSE, ii, jj, kk));
1315
1316 //* **********************************************************
1317 //* *** > dPSW;
1318 //* **********************************************************
1319
1320 VAT3(dPSW, ii, jj, kk) =
1321 (
1322 VAT3( uNE, im1, jm1, km1)
1323 + VAT3( uE, im1, jm1, km1)
1324 * VAT3( oPS, ii, jj, kk)
1325 + VAT3( uN, im1, jm1, km1)
1326 * VAT3( oPW, ii, jj, kk)
1327 + VAT3( uC, im1, jm1, km1)
1328 * VAT3(oPSW, ii, jj, kk)
1329 + VAT3( oNE, im1, jm1, km1)
1330 * VAT3( dPC, ii, jj, kk)
1331 + VAT3( oE, im1, jm1, km1)
1332 * VAT3( dPS, ii, jj, kk)
1333 + VAT3( oN, im1, jm1, km1)
1334 * VAT3( dPW, ii, jj, kk)
1335 )
1336 / VAT3( oC, im1, jm1, km1);
1337
1338 //fprintf(data, "%19.12E\n", VAT3(dPSW, ii, jj, kk));
1339
1340 //* **********************************************************
1341 //* *** > uPC;
1342 //* **********************************************************
1343
1344 VAT3( uPC, ii, jj, kk) =
1345 (
1346 VAT3( uSE, im1, jp1, k)
1347 + VAT3( uE, im1, j, k)
1348 + VAT3( uNE, im1, jm1, k)
1349 + VAT3( uS, i, jp1, k)
1350 + VAT3( uC, i, j, k)
1351 + VAT3( uN, i, jm1, k)
1352 + VAT3( uSW, ip1, jp1, k)
1353 + VAT3( uW, ip1, j, k)
1354 + VAT3( uNW, ip1, jm1, k)
1355 ) / (
1356 VAT3( oC, i, j, kp1)
1357 - VAT3( oN, i, j, kp1)
1358 - VAT3( oN, i, jm1, kp1)
1359 - VAT3( oNW, i, j, kp1)
1360 - VAT3( oE, im1, j, kp1)
1361 - VAT3( oNE, im1, jm1, kp1)
1362 - VAT3( oNE, i, j, kp1)
1363 - VAT3( oE, i, j, kp1)
1364 - VAT3( oNW, ip1, jm1, kp1)
1365 );
1366
1367 //fprintf(data, "%19.12E\n", VAT3(uPC, ii, jj, kk));
1368
1369 //* **********************************************************
1370 //* *** > uPN;
1371 //* **********************************************************
1372
1373 VAT3( uPN, ii, jj, kk) =
1374 (
1375 VAT3( uNE, im1, j, k)
1376 + VAT3( uN, i, j, k)
1377 + VAT3( uNW, ip1, j, k)
1378 + (
1379 VAT3( oNE, im1, j, kp1)
1380 + VAT3( oN, i, j, kp1)
1381 + VAT3( oNW, ip1, j, kp1)
1382 )
1383 * VAT3( uPC, ii, jj, kk)
1384 + (
1385 VAT3( uE, im1, jp1, k)
1386 + VAT3( uC, i, jp1, k)
1387 + VAT3( uW, ip1, jp1, k)
1388 )
1389 * VAT3( oPN, ii, jj, kk)
1390 ) / (
1391 VAT3( oC, i, jp1, kp1)
1392 - VAT3( oE, im1, jp1, kp1)
1393 - VAT3( oE, i, jp1, kp1)
1394 );
1395
1396 //fprintf(data, "%19.12E\n", VAT3(uPN, ii, jj, kk));
1397
1398 //* **********************************************************
1399 //* *** > uPS;
1400 //* **********************************************************
1401
1402 VAT3( uPS, ii, jj, kk) =
1403 (
1404 VAT3( uSE, im1, j, k)
1405 + VAT3( uS, i, j, k)
1406 + VAT3( uSW, ip1, j, k)
1407 + (
1408 VAT3( oNW, i, jm1, kp1)
1409 + VAT3( oN, i, jm1, kp1)
1410 + VAT3( oNE, i, jm1, kp1)
1411 )
1412 * VAT3( uPC, ii, jj, kk)
1413 + (
1414 VAT3( uE, im1, jm1, k)
1415 + VAT3( uC, i, jm1, k)
1416 + VAT3( uW, ip1, jm1, k)
1417 )
1418 * VAT3( oPS, ii, jj, kk)
1419 ) / (
1420 VAT3( oC, i, jm1, kp1)
1421 - VAT3( oE, im1, jm1, kp1)
1422 - VAT3( oE, i, jm1, kp1)
1423 );
1424
1425 //fprintf(data, "%19.12E\n", VAT3(uPS, ii, jj, kk));
1426
1427 //* **********************************************************
1428 //* *** > uPE;
1429 //* **********************************************************
1430
1431 VAT3( uPE, ii, jj, kk) =
1432 (
1433 VAT3( uSE, i, jp1, k)
1434 + VAT3( uS, ip1, jp1, k)
1435 + VAT3( uNE, i, jm1, k)
1436 + (
1437 VAT3( uS, ip1, jp1, k)
1438 + VAT3( uC, ip1, j, k)
1439 + VAT3( uN, ip1, jm1, k)
1440 )
1441 * VAT3( oPE, ii, jj, kk)
1442 + (
1443 VAT3( oNW, ip1, j, kp1)
1444 + VAT3( oE, i, j, kp1)
1445 + VAT3( oNE, i, jm1, kp1)
1446 )
1447 * VAT3( uPC, ii, jj, kk)
1448 ) / (
1449 VAT3( oC, ip1, j, kp1)
1450 - VAT3( oN, ip1, j, kp1)
1451 - VAT3( oN, ip1, jm1, kp1)
1452 );
1453
1454 //fprintf(data, "%19.12E\n", VAT3(uPE, ii, jj, kk));
1455
1456 //* **********************************************************
1457 //* *** > uPW;
1458 //* **********************************************************
1459
1460 VAT3( uPW, ii, jj, kk) =
1461 (
1462 VAT3( uSW, i, jp1, k)
1463 + VAT3( uW, i, j, k)
1464 + VAT3( uNW, i, jm1, k)
1465 + (
1466 VAT3( uS, im1, jp1, k)
1467 + VAT3( uC, im1, j, k)
1468 + VAT3( uN, im1, jm1, k)
1469 )
1470 * VAT3( oPW, ii, jj, kk)
1471 + (
1472 VAT3( oNE, im1, j, kp1)
1473 + VAT3( oE, im1, j, kp1)
1474 + VAT3( oNW, i, jm1, kp1)
1475 )
1476 * VAT3( uPC, ii, jj, kk)
1477 ) / (
1478 VAT3( oC, im1, j, kp1)
1479 - VAT3( oN, im1, j, kp1)
1480 - VAT3( oN, im1, jm1, kp1)
1481 );
1482
1483 //fprintf(data, "%19.12E\n", VAT3(uPW, ii, jj, kk));
1484
1485 //* **********************************************************
1486 //* *** > uPNE;
1487 //* **********************************************************
1488
1489 VAT3(uPNE, ii, jj, kk) =
1490 (
1491 VAT3( uNE, i, j, k)
1492 + VAT3( uE, i, jp1, k)
1493 * VAT3( oPN, ii, jj, kk)
1494 + VAT3( uN, ip1, j, k)
1495 * VAT3( oPE, ii, jj, kk)
1496 + VAT3( uC, ip1, jp1, k)
1497 * VAT3(oPNE, ii, jj, kk)
1498 + VAT3( oNE, i, j, kp1)
1499 * VAT3( uPC, ii, jj, kk)
1500 + VAT3( oE, i, jp1, kp1)
1501 * VAT3( uPN, ii, jj, kk)
1502 + VAT3( oN, ip1, j, kp1)
1503 * VAT3( uPE, ii, jj, kk)
1504 )
1505 / VAT3( oC, ip1, jp1, kp1);
1506
1507 //fprintf(data, "%19.12E\n", VAT3(uPNE, ii, jj, kk));
1508
1509 //* **********************************************************
1510 //* *** > uPNW;
1511 //* **********************************************************
1512
1513 VAT3(uPNW, ii, jj, kk) =
1514 (
1515 VAT3( uNW, i, j, k)
1516 + VAT3( uW, i, jp1, k)
1517 * VAT3( oPN, ii, jj, kk)
1518 + VAT3( uN, im1, j, k)
1519 * VAT3( oPW, ii, jj, kk)
1520 + VAT3( uC, im1, jp1, k)
1521 * VAT3(oPNW, ii, jj, kk)
1522 + VAT3( oNW, i, j, kp1)
1523 * VAT3( uPC, ii, jj, kk)
1524 + VAT3( oE, im1, jp1, kp1)
1525 * VAT3( uPN, ii, jj, kk)
1526 + VAT3( oN, im1, j, kp1)
1527 * VAT3( uPW, ii, jj, kk)
1528 )
1529 / VAT3( oC, im1, jp1, kp1);
1530
1531 //fprintf(data, "%19.12E\n", VAT3(uPNW, ii, jj, kk));
1532
1533 //* **********************************************************
1534 //* *** > uPSE;
1535 //* **********************************************************
1536
1537 VAT3(uPSE, ii, jj, kk) =
1538 (
1539 VAT3( uSE, i, j, k)
1540 + VAT3( uE, i, jm1, k)
1541 * VAT3( oPS, ii, jj, kk)
1542 + VAT3( uS, ip1, j, k)
1543 * VAT3( oPE, ii, jj, kk)
1544 + VAT3( uC, ip1, jm1, k)
1545 * VAT3(oPSE, ii, jj, kk)
1546 + VAT3( oNW, ip1, jm1, kp1)
1547 * VAT3( uPC, ii, jj, kk)
1548 + VAT3( oE, i, jm1, kp1)
1549 * VAT3( uPS, ii, jj, kk)
1550 + VAT3( oN, ip1, jm1, kp1)
1551 * VAT3( uPE, ii, jj, kk)
1552 )
1553 / VAT3( oC, ip1, jm1, kp1);
1554
1555 //fprintf(data, "%19.12E\n", VAT3(uPSE, ii, jj, kk));
1556
1557 //* **********************************************************
1558 //* *** > uPSW;
1559 //* **********************************************************
1560
1561 VAT3(uPSW, ii, jj, kk) =
1562 (
1563 VAT3( uSW, i, j, k)
1564 + VAT3( uW, i, jm1, k)
1565 * VAT3( oPS, ii, jj, kk)
1566 + VAT3( uS, im1, j, k)
1567 * VAT3( oPW, ii, jj, kk)
1568 + VAT3( uC, im1, jm1, k)
1569 * VAT3(oPSW, ii, jj, kk)
1570 + VAT3( oNE, im1, jm1, kp1)
1571 * VAT3( uPC, ii, jj, kk)
1572 + VAT3( oE, im1, jm1, kp1)
1573 * VAT3( uPS, ii, jj, kk)
1574 + VAT3( oN, im1, jm1, kp1)
1575 * VAT3( uPW, ii, jj, kk)
1576 )
1577 / VAT3( oC, im1, jm1, kp1);
1578
1579 //fprintf(data, "%19.12E\n", VAT3(uPSW, ii, jj, kk));
1580
1581 }
1582 }
1583 }
1584}
VPUBLIC void VbuildP(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, int *mgprol, int *ipc, double *rpc, double *pc, double *ac, double *xf, double *yf, double *zf)
Builds prolongation matrix.
Definition buildPd.c:57