APBS 3.0.0
Loading...
Searching...
No Matches
buildAd.c
1
55#include "buildAd.h"
56
57VPUBLIC 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
105VPUBLIC 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
285VPUBLIC 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