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) {
66 MAT2(ac, *nx * *ny * *nz, 14);
70 VbuildA_fv(nx, ny, nz,
73 RAT2(ac, 1,1), cc, fc,
74 RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
80 }
else if (*mgdisc == 1) {
82 VbuildA_fe(nx, ny, nz,
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),
96 Vnm_print(2,
"VbuildA: Invalid discretization requested.\n");
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) {
129 int nxm1, nym1, nzm1;
135 double hxm1, hym1, hzm1;
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);
170 VAT(ipc, 10) = *ipkey;
195 for(k=2; k<=*nz-1; k++) {
197 hzm1 = VAT(zf, k) - VAT(zf, k-1);
198 hz = VAT(zf, k+1) - VAT(zf, k);
200 for(j=2; j<=*ny-1; j++) {
202 hym1 = VAT(yf, j) - VAT(yf, j-1);
203 hy = VAT(yf, j+1) - VAT(yf, j);
205 for(i=2; i<=*nx-1; i++) {
207 hxm1 = VAT(xf, i) - VAT(xf, i-1);
208 hx = VAT(xf, i+1) - VAT(xf, i);
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;
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);
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);
242 ike = VMIN2(1, VABS(i - nxm1));
243 VAT3(oE, i, j, k) = ike * coef_oE * VAT3(a1cf, 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;
249 jke = VMIN2(1, VABS(j - nym1));
250 VAT3(oN, i, j, k) = jke * coef_oN * VAT3(a2cf, 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;
256 kke = VMIN2(1, VABS(k - nzm1));
257 VAT3(uC, i, j, k) = kke * coef_uC * VAT3(a3cf, 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;
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;
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;
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;
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");
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.