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