APBS  3.0.0
buildPd.c
1 
55 #include "buildPd.h"
56 
57 VPUBLIC 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 
96 VPUBLIC 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 
114 VEXTERNC 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 
165 for(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 
206 VPUBLIC 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 
231 VPUBLIC 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 
722 VPUBLIC 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 
749 VPUBLIC 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