APBS  1.5
buildAd.c
1 
50 #include "buildAd.h"
51 
52 VPUBLIC void VbuildA(int* nx, int* ny, int* nz,
53  int* ipkey, int* mgdisc, int* numdia,
54  int* ipc, double* rpc,
55  double* ac, double* cc, double* fc,
56  double* xf, double* yf, double* zf,
57  double* gxcf, double* gycf, double* gzcf,
58  double* a1cf, double* a2cf, double* a3cf,
59  double* ccf, double* fcf) {
60 
61  MAT2(ac, *nx * *ny * *nz, 14);
62 
63  if (*mgdisc == 0) {
64 
65  VbuildA_fv(nx, ny, nz,
66  ipkey, numdia,
67  ipc, rpc,
68  RAT2(ac, 1,1), cc, fc,
69  RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
70  xf, yf, zf,
71  gxcf, gycf, gzcf,
72  a1cf, a2cf, a3cf,
73  ccf, fcf);
74 
75  } else if (*mgdisc == 1) {
76 
77  VbuildA_fe(nx, ny, nz,
78  ipkey, numdia,
79  ipc, rpc,
80  RAT2(ac, 1, 1), cc, fc,
81  RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4), RAT2(ac, 1, 5), RAT2(ac, 1, 6),
82  RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
83  RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
84  xf, yf, zf,
85  gxcf, gycf, gzcf,
86  a1cf, a2cf, a3cf,
87  ccf, fcf);
88 
89  } else {
90 
91  Vnm_print(2, "VbuildA: Invalid discretization requested.\n");
93  exit(EXIT_FAILURE);
94 
95  }
96 }
97 
98 
99 
100 VPUBLIC void VbuildA_fv(int *nx, int *ny, int *nz,
101  int *ipkey, int *numdia,
102  int *ipc, double *rpc,
103  double *oC, double *cc, double *fc, double *oE, double *oN, double *uC,
104  double *xf, double *yf, double *zf,
105  double *gxcf, double *gycf, double *gzcf,
106  double *a1cf, double *a2cf, double *a3cf,
107  double *ccf, double *fcf) {
108 
109  int i, j, k; // @todo Document this function
110 
121  int ike, jke, kke;
122 
123 
124  int nxm1, nym1, nzm1;
125 
126 
127  double hx, hy, hz;
128 
129 
130  double hxm1, hym1, hzm1;
131 
132  double coef_fc;
133 
134  double bc_cond_e;
135  double bc_cond_w;
136  double bc_cond_n;
137  double bc_cond_s;
138  double bc_cond_u;
139  double bc_cond_d;
140  double coef_oE;
141  double coef_oN;
142  double coef_uC;
143  double coef_oEm1;
144  double coef_oNm1;
145  double coef_uCm1;
146 
147  double diag;
148 
149  MAT3( fc, *nx, *ny, *nz);
150  MAT3( fcf, *nx, *ny, *nz);
151  MAT3( cc, *nx, *ny, *nz);
152  MAT3( ccf, *nx, *ny, *nz);
153  MAT3( oC, *nx, *ny, *nz);
154  MAT3(a1cf, *nx, *ny, *nz);
155  MAT3(a2cf, *nx, *ny, *nz);
156  MAT3(a3cf, *nx, *ny, *nz);
157  MAT3( uC, *nx, *ny, *nz);
158  MAT3( oN, *nx, *ny, *nz);
159  MAT3( oE, *nx, *ny, *nz);
160  MAT3(gxcf, *ny, *nz, 2);
161  MAT3(gycf, *nx, *nz, 2);
162  MAT3(gzcf, *nx, *ny, 2);
163 
164  // Save the problem key with this operator. @todo: What?
165  VAT(ipc, 10) = *ipkey;
166 
167  // Note how many nonzeros in this discretization stencil
168  VAT(ipc, 11) = 7;
169  VAT(ipc, 12) = 1;
170  *numdia = 4;
171 
172  // Define n and determine number of mesh points
173  nxm1 = *nx - 1;
174  nym1 = *ny - 1;
175  nzm1 = *nz - 1;
176 
177  // Determine diag scale factor
178  // (would like something close to ones on the main diagonal)
179  // @todo: Make a more meaningful comment
180  diag = 1.0;
181 
182 
183 
184  /* *********************************************************************
185  * *** interior points ***
186  * ********************************************************************* */
187 
188  // build the operator
189  //fprintf(data, "%s\n", PRINT_FUNC);
190  for(k=2; k<=*nz-1; k++) {
191 
192  hzm1 = VAT(zf, k) - VAT(zf, k-1);
193  hz = VAT(zf, k+1) - VAT(zf, k);
194 
195  for(j=2; j<=*ny-1; j++) {
196 
197  hym1 = VAT(yf, j) - VAT(yf, j-1);
198  hy = VAT(yf, j+1) - VAT(yf, j);
199 
200  for(i=2; i<=*nx-1; i++) {
201 
202  hxm1 = VAT(xf, i) - VAT(xf, i-1);
203  hx = VAT(xf, i+1) - VAT(xf, i);
204 
205  // Calculate some coefficients
212  coef_oE = diag * (hym1 + hy) * (hzm1 + hz) / (4.0 * hx);
213  coef_oEm1 = diag * (hym1 + hy) * (hzm1 + hz) / (4.0 * hxm1);
214  coef_oN = diag * (hxm1 + hx) * (hzm1 + hz) / (4.0 * hy);
215  coef_oNm1 = diag * (hxm1 + hx) * (hzm1 + hz) / (4.0 * hym1);
216  coef_uC = diag * (hxm1 + hx) * (hym1 + hy) / (4.0 * hz);
217  coef_uCm1 = diag * (hxm1 + hx) * (hym1 + hy) / (4.0 * hzm1);
218  coef_fc = diag * (hxm1 + hx) * (hym1 + hy) * (hzm1 + hz) / 8.0;
219 
220  // Calculate the coefficient and source function
221  VAT3(fc, i, j, k) = coef_fc * VAT3(fcf, i, j, k);
222  VAT3(cc, i, j, k) = coef_fc * VAT3(ccf, i, j, k);
223  //fprintf(data, "%19.12E\n", VAT3(cc, i, j, k));
224 
225  // Calculate the diagonal for matvecs and smoothings
226 
227  VAT3(oC, i, j, k) = coef_oE * VAT3(a1cf, i, j, k) +
228  coef_oEm1 * VAT3(a1cf, i-1, j, k) +
229  coef_oN * VAT3(a2cf, i, j, k) +
230  coef_oNm1 * VAT3(a2cf, i, j-1, k) +
231  coef_uC * VAT3(a3cf, i, j, k) +
232  coef_uCm1 * VAT3(a3cf, i, j, k-1);
233 
234  //fprintf(data, "%19.12E\n", VAT3(oC, i, j, k));
235 
236  // Calculate the east neighbor
237  ike = VMIN2(1, VABS(i - nxm1));
238  VAT3(oE, i, j, k) = ike * coef_oE * VAT3(a1cf, i, j, k);
239  //fprintf(data, "%19.12E\n", VAT3(oE, i, j, k));
240  bc_cond_e = (1 - ike) * coef_oE * VAT3(a1cf, i, j, k) * VAT3(gxcf, j, k, 2);
241  VAT3(fc, i, j, k) += bc_cond_e;
242 
243  // Calculate the north neighbor
244  jke = VMIN2(1, VABS(j - nym1));
245  VAT3(oN, i, j, k) = jke * coef_oN * VAT3(a2cf, i, j, k);
246  //fprintf(data, "%19.12E\n", VAT3(oN, i, j, k));
247  bc_cond_n = (1 - jke) * coef_oN * VAT3(a2cf, i, j, k) * VAT3(gycf, i, k, 2);
248  VAT3(fc, i, j, k) += bc_cond_n;
249 
250  // Calculate the up neighbor
251  kke = VMIN2(1, VABS(k - nzm1));
252  VAT3(uC, i, j, k) = kke * coef_uC * VAT3(a3cf, i, j, k);
253  //fprintf(data, "%19.12E\n", VAT3(uC, i, j, k));
254  bc_cond_u = (1 - kke) * coef_uC * VAT3(a3cf, i, j, k) * VAT3(gzcf, i, j, 2);
255  VAT3(fc, i, j, k) += bc_cond_u;
256 
257  // Calculate the west neighbor (just handle b.c.)
258  ike = VMIN2(1, VABS(i - 2));
259  bc_cond_w = (1 - ike) * coef_oEm1 * VAT3(a1cf, i-1, j, k) * VAT3(gxcf, j, k, 1);
260  VAT3(fc, i, j, k) += bc_cond_w;
261 
262  // Calculate the south neighbor (just handle b.c.)
263  jke = VMIN2(1, VABS(j - 2));
264  bc_cond_s = (1 - jke) * coef_oNm1 * VAT3(a2cf, i, j-1, k) * VAT3(gycf, i, k, 1);
265  VAT3(fc, i, j, k) += bc_cond_s;
266 
267  // Calculate the down neighbor (just handle b.c.)
268  kke = VMIN2(1, VABS(k - 2));
269  bc_cond_d = (1 - kke) * coef_uCm1 * VAT3(a3cf, i, j, k-1) * VAT3(gzcf, i, j, 1);
270  VAT3(fc, i, j, k) += bc_cond_d;
271 
272  //fprintf(data, "%19.12E\n", VAT3(fc, i, j, k));
273  }
274  }
275  }
276 }
277 
278 
279 
280 VPUBLIC void VbuildA_fe(int *nx, int *ny, int *nz,
281  int *ipkey, int *numdia,
282  int *ipc, double *rpc,
283  double *oC, double *cc, double *fc,
284  double *oE, double *oN, double *uC,
285  double* oNE, double* oNW,
286  double* uE, double* uW,
287  double* uN, double* uS,
288  double* uNE, double* uNW, double* uSE, double* uSW,
289  double *xf, double *yf, double *zf,
290  double *gxcf, double *gycf, double *gzcf,
291  double *a1cf, double *a2cf, double *a3cf,
292  double *ccf, double *fcf) {
293  VABORT_MSG0("Untranslated Component: from buildAd.f");
294 }
double hx
Definition: vpmgp.h:87
int ipkey
Definition: vpmgp.h:109
double hy
Definition: vpmgp.h:88
int mgdisc
Definition: vpmgp.h:177
VPUBLIC void VbuildA(int *nx, int *ny, int *nz, int *ipkey, int *mgdisc, int *numdia, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *xf, double *yf, double *zf, double *gxcf, double *gycf, double *gzcf, double *a1cf, double *a2cf, double *a3cf, double *ccf, double *fcf)
Build the Laplacian.
Definition: buildAd.c:52
int ny
Definition: vpmgp.h:84
int nx
Definition: vpmgp.h:83
int nz
Definition: vpmgp.h:85