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