Actual source code: ex30.c
1: static char help[] = "Test DMStagMatGetValuesStencil() in 3D\n\n";
3: #include <petscdm.h>
4: #include <petscksp.h>
5: #include <petscdmstag.h>
7: /* Shorter, more convenient names for DMStagStencilLocation entries */
8: #define BACK_DOWN_LEFT DMSTAG_BACK_DOWN_LEFT
9: #define BACK_DOWN DMSTAG_BACK_DOWN
10: #define BACK_DOWN_RIGHT DMSTAG_BACK_DOWN_RIGHT
11: #define BACK_LEFT DMSTAG_BACK_LEFT
12: #define BACK DMSTAG_BACK
13: #define BACK_RIGHT DMSTAG_BACK_RIGHT
14: #define BACK_UP_LEFT DMSTAG_BACK_UP_LEFT
15: #define BACK_UP DMSTAG_BACK_UP
16: #define BACK_UP_RIGHT DMSTAG_BACK_UP_RIGHT
17: #define DOWN_LEFT DMSTAG_DOWN_LEFT
18: #define DOWN DMSTAG_DOWN
19: #define DOWN_RIGHT DMSTAG_DOWN_RIGHT
20: #define LEFT DMSTAG_LEFT
21: #define ELEMENT DMSTAG_ELEMENT
22: #define RIGHT DMSTAG_RIGHT
23: #define UP_LEFT DMSTAG_UP_LEFT
24: #define UP DMSTAG_UP
25: #define UP_RIGHT DMSTAG_UP_RIGHT
26: #define FRONT_DOWN_LEFT DMSTAG_FRONT_DOWN_LEFT
27: #define FRONT_DOWN DMSTAG_FRONT_DOWN
28: #define FRONT_DOWN_RIGHT DMSTAG_FRONT_DOWN_RIGHT
29: #define FRONT_LEFT DMSTAG_FRONT_LEFT
30: #define FRONT DMSTAG_FRONT
31: #define FRONT_RIGHT DMSTAG_FRONT_RIGHT
32: #define FRONT_UP_LEFT DMSTAG_FRONT_UP_LEFT
33: #define FRONT_UP DMSTAG_FRONT_UP
34: #define FRONT_UP_RIGHT DMSTAG_FRONT_UP_RIGHT
36: static PetscErrorCode CreateMat(DM, Mat *);
37: static PetscErrorCode CheckMat(DM, Mat);
39: int main(int argc, char **argv)
40: {
41: DM dmSol;
42: Mat A;
45: PetscInitialize(&argc, &argv, (char *)0, help);
46: {
47: const PetscInt dof0 = 0, dof1 = 0, dof2 = 1, dof3 = 1; /* 1 dof on each face and element center */
48: const PetscInt stencilWidth = 1;
49: DMStagCreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 4, 5, 6, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, dof3, DMSTAG_STENCIL_BOX, stencilWidth, NULL, NULL, NULL, &dmSol);
50: DMSetFromOptions(dmSol);
51: DMSetUp(dmSol);
52: DMStagSetUniformCoordinatesExplicit(dmSol, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
53: }
54: CreateMat(dmSol, &A);
55: CheckMat(dmSol, A);
56: MatDestroy(&A);
57: DMDestroy(&dmSol);
58: PetscFinalize();
59: return 0;
60: }
62: static PetscErrorCode CreateMat(DM dmSol, Mat *pA)
63: {
64: Vec coordLocal;
65: Mat A;
66: PetscInt startx, starty, startz, N[3], nx, ny, nz, ex, ey, ez, d;
67: PetscInt icp[3], icux[3], icuy[3], icuz[3], icux_right[3], icuy_up[3], icuz_front[3];
68: PetscBool isLastRankx, isLastRanky, isLastRankz, isFirstRankx, isFirstRanky, isFirstRankz;
69: PetscReal hx, hy, hz;
70: DM dmCoord;
71: PetscScalar ****arrCoord;
74: DMCreateMatrix(dmSol, pA);
75: A = *pA;
76: DMStagGetCorners(dmSol, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL);
77: DMStagGetGlobalSizes(dmSol, &N[0], &N[1], &N[2]);
79: DMStagGetIsLastRank(dmSol, &isLastRankx, &isLastRanky, &isLastRankz);
80: DMStagGetIsFirstRank(dmSol, &isFirstRankx, &isFirstRanky, &isFirstRankz);
81: hx = 1.0 / N[0];
82: hy = 1.0 / N[1];
83: hz = 1.0 / N[2];
84: DMGetCoordinateDM(dmSol, &dmCoord);
85: DMGetCoordinatesLocal(dmSol, &coordLocal);
86: DMStagVecGetArrayRead(dmCoord, coordLocal, &arrCoord);
87: for (d = 0; d < 3; ++d) {
88: DMStagGetLocationSlot(dmCoord, ELEMENT, d, &icp[d]);
89: DMStagGetLocationSlot(dmCoord, LEFT, d, &icux[d]);
90: DMStagGetLocationSlot(dmCoord, DOWN, d, &icuy[d]);
91: DMStagGetLocationSlot(dmCoord, BACK, d, &icuz[d]);
92: DMStagGetLocationSlot(dmCoord, RIGHT, d, &icux_right[d]);
93: DMStagGetLocationSlot(dmCoord, UP, d, &icuy_up[d]);
94: DMStagGetLocationSlot(dmCoord, FRONT, d, &icuz_front[d]);
95: }
97: for (ez = startz; ez < startz + nz; ++ez) {
98: for (ey = starty; ey < starty + ny; ++ey) {
99: for (ex = startx; ex < startx + nx; ++ex) {
100: if (ex == N[0] - 1) {
101: /* Right Boundary velocity Dirichlet */
102: DMStagStencil row;
103: const PetscScalar valA = 1.0;
104: row.i = ex;
105: row.j = ey;
106: row.k = ez;
107: row.loc = RIGHT;
108: row.c = 0;
109: DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
110: }
111: if (ey == N[1] - 1) {
112: /* Top boundary velocity Dirichlet */
113: DMStagStencil row;
114: const PetscScalar valA = 1.0;
115: row.i = ex;
116: row.j = ey;
117: row.k = ez;
118: row.loc = UP;
119: row.c = 0;
120: DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
121: }
122: if (ez == N[2] - 1) {
123: /* Top boundary velocity Dirichlet */
124: DMStagStencil row;
125: const PetscScalar valA = 1.0;
126: row.i = ex;
127: row.j = ey;
128: row.k = ez;
129: row.loc = FRONT;
130: row.c = 0;
131: DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
132: }
134: /* Equation on left face of this element */
135: if (ex == 0) {
136: /* Left velocity Dirichlet */
137: DMStagStencil row;
138: const PetscScalar valA = 1.0;
139: row.i = ex;
140: row.j = ey;
141: row.k = ez;
142: row.loc = LEFT;
143: row.c = 0;
144: DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
145: } else {
146: /* X-momentum interior equation : (u_xx + u_yy + u_zz) - p_x = f^x */
147: DMStagStencil row, col[9];
148: PetscScalar valA[9];
149: PetscInt nEntries;
151: row.i = ex;
152: row.j = ey;
153: row.k = ez;
154: row.loc = LEFT;
155: row.c = 0;
156: if (ey == 0) {
157: if (ez == 0) {
158: nEntries = 7;
159: col[0].i = ex;
160: col[0].j = ey;
161: col[0].k = ez;
162: col[0].loc = LEFT;
163: col[0].c = 0;
164: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
165: /* Missing down term */
166: col[1].i = ex;
167: col[1].j = ey + 1;
168: col[1].k = ez;
169: col[1].loc = LEFT;
170: col[1].c = 0;
171: valA[1] = 1.0 / (hy * hy);
172: col[2].i = ex - 1;
173: col[2].j = ey;
174: col[2].k = ez;
175: col[2].loc = LEFT;
176: col[2].c = 0;
177: valA[2] = 1.0 / (hx * hx);
178: col[3].i = ex + 1;
179: col[3].j = ey;
180: col[3].k = ez;
181: col[3].loc = LEFT;
182: col[3].c = 0;
183: valA[3] = 1.0 / (hx * hx);
184: /* Missing back term */
185: col[4].i = ex;
186: col[4].j = ey;
187: col[4].k = ez + 1;
188: col[4].loc = LEFT;
189: col[4].c = 0;
190: valA[4] = 1.0 / (hz * hz);
191: col[5].i = ex - 1;
192: col[5].j = ey;
193: col[5].k = ez;
194: col[5].loc = ELEMENT;
195: col[5].c = 0;
196: valA[5] = 1.0 / hx;
197: col[6].i = ex;
198: col[6].j = ey;
199: col[6].k = ez;
200: col[6].loc = ELEMENT;
201: col[6].c = 0;
202: valA[6] = -1.0 / hx;
203: } else if (ez == N[2] - 1) {
204: nEntries = 7;
205: col[0].i = ex;
206: col[0].j = ey;
207: col[0].k = ez;
208: col[0].loc = LEFT;
209: col[0].c = 0;
210: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
211: /* Missing down term */
212: col[1].i = ex;
213: col[1].j = ey + 1;
214: col[1].k = ez;
215: col[1].loc = LEFT;
216: col[1].c = 0;
217: valA[1] = 1.0 / (hy * hy);
218: col[2].i = ex - 1;
219: col[2].j = ey;
220: col[2].k = ez;
221: col[2].loc = LEFT;
222: col[2].c = 0;
223: valA[2] = 1.0 / (hx * hx);
224: col[3].i = ex + 1;
225: col[3].j = ey;
226: col[3].k = ez;
227: col[3].loc = LEFT;
228: col[3].c = 0;
229: valA[3] = 1.0 / (hx * hx);
230: col[4].i = ex;
231: col[4].j = ey;
232: col[4].k = ez - 1;
233: col[4].loc = LEFT;
234: col[4].c = 0;
235: valA[4] = 1.0 / (hz * hz);
236: /* Missing front term */
237: col[5].i = ex - 1;
238: col[5].j = ey;
239: col[5].k = ez;
240: col[5].loc = ELEMENT;
241: col[5].c = 0;
242: valA[5] = 1.0 / hx;
243: col[6].i = ex;
244: col[6].j = ey;
245: col[6].k = ez;
246: col[6].loc = ELEMENT;
247: col[6].c = 0;
248: valA[6] = -1.0 / hx;
249: } else {
250: nEntries = 8;
251: col[0].i = ex;
252: col[0].j = ey;
253: col[0].k = ez;
254: col[0].loc = LEFT;
255: col[0].c = 0;
256: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
257: /* Missing down term */
258: col[1].i = ex;
259: col[1].j = ey + 1;
260: col[1].k = ez;
261: col[1].loc = LEFT;
262: col[1].c = 0;
263: valA[1] = 1.0 / (hy * hy);
264: col[2].i = ex - 1;
265: col[2].j = ey;
266: col[2].k = ez;
267: col[2].loc = LEFT;
268: col[2].c = 0;
269: valA[2] = 1.0 / (hx * hx);
270: col[3].i = ex + 1;
271: col[3].j = ey;
272: col[3].k = ez;
273: col[3].loc = LEFT;
274: col[3].c = 0;
275: valA[3] = 1.0 / (hx * hx);
276: col[4].i = ex;
277: col[4].j = ey;
278: col[4].k = ez - 1;
279: col[4].loc = LEFT;
280: col[4].c = 0;
281: valA[4] = 1.0 / (hz * hz);
282: col[5].i = ex;
283: col[5].j = ey;
284: col[5].k = ez + 1;
285: col[5].loc = LEFT;
286: col[5].c = 0;
287: valA[5] = 1.0 / (hz * hz);
288: col[6].i = ex - 1;
289: col[6].j = ey;
290: col[6].k = ez;
291: col[6].loc = ELEMENT;
292: col[6].c = 0;
293: valA[6] = 1.0 / hx;
294: col[7].i = ex;
295: col[7].j = ey;
296: col[7].k = ez;
297: col[7].loc = ELEMENT;
298: col[7].c = 0;
299: valA[7] = -1.0 / hx;
300: }
301: } else if (ey == N[1] - 1) {
302: if (ez == 0) {
303: nEntries = 7;
304: col[0].i = ex;
305: col[0].j = ey;
306: col[0].k = ez;
307: col[0].loc = LEFT;
308: col[0].c = 0;
309: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
310: col[1].i = ex;
311: col[1].j = ey - 1;
312: col[1].k = ez;
313: col[1].loc = LEFT;
314: col[1].c = 0;
315: valA[1] = 1.0 / (hy * hy);
316: /* Missing up term */
317: col[2].i = ex - 1;
318: col[2].j = ey;
319: col[2].k = ez;
320: col[2].loc = LEFT;
321: col[2].c = 0;
322: valA[2] = 1.0 / (hx * hx);
323: col[3].i = ex + 1;
324: col[3].j = ey;
325: col[3].k = ez;
326: col[3].loc = LEFT;
327: col[3].c = 0;
328: valA[3] = 1.0 / (hx * hx);
329: /* Missing back entry */
330: col[4].i = ex;
331: col[4].j = ey;
332: col[4].k = ez + 1;
333: col[4].loc = LEFT;
334: col[4].c = 0;
335: valA[4] = 1.0 / (hz * hz);
336: col[5].i = ex - 1;
337: col[5].j = ey;
338: col[5].k = ez;
339: col[5].loc = ELEMENT;
340: col[5].c = 0;
341: valA[5] = 1.0 / hx;
342: col[6].i = ex;
343: col[6].j = ey;
344: col[6].k = ez;
345: col[6].loc = ELEMENT;
346: col[6].c = 0;
347: valA[6] = -1.0 / hx;
348: } else if (ez == N[2] - 1) {
349: nEntries = 7;
350: col[0].i = ex;
351: col[0].j = ey;
352: col[0].k = ez;
353: col[0].loc = LEFT;
354: col[0].c = 0;
355: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
356: col[1].i = ex;
357: col[1].j = ey - 1;
358: col[1].k = ez;
359: col[1].loc = LEFT;
360: col[1].c = 0;
361: valA[1] = 1.0 / (hy * hy);
362: /* Missing up term */
363: col[2].i = ex - 1;
364: col[2].j = ey;
365: col[2].k = ez;
366: col[2].loc = LEFT;
367: col[2].c = 0;
368: valA[2] = 1.0 / (hx * hx);
369: col[3].i = ex + 1;
370: col[3].j = ey;
371: col[3].k = ez;
372: col[3].loc = LEFT;
373: col[3].c = 0;
374: valA[3] = 1.0 / (hx * hx);
375: col[4].i = ex;
376: col[4].j = ey;
377: col[4].k = ez - 1;
378: col[4].loc = LEFT;
379: col[4].c = 0;
380: valA[4] = 1.0 / (hz * hz);
381: /* Missing front term */
382: col[5].i = ex - 1;
383: col[5].j = ey;
384: col[5].k = ez;
385: col[5].loc = ELEMENT;
386: col[5].c = 0;
387: valA[5] = 1.0 / hx;
388: col[6].i = ex;
389: col[6].j = ey;
390: col[6].k = ez;
391: col[6].loc = ELEMENT;
392: col[6].c = 0;
393: valA[6] = -1.0 / hx;
394: } else {
395: nEntries = 8;
396: col[0].i = ex;
397: col[0].j = ey;
398: col[0].k = ez;
399: col[0].loc = LEFT;
400: col[0].c = 0;
401: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
402: col[1].i = ex;
403: col[1].j = ey - 1;
404: col[1].k = ez;
405: col[1].loc = LEFT;
406: col[1].c = 0;
407: valA[1] = 1.0 / (hy * hy);
408: /* Missing up term */
409: col[2].i = ex - 1;
410: col[2].j = ey;
411: col[2].k = ez;
412: col[2].loc = LEFT;
413: col[2].c = 0;
414: valA[2] = 1.0 / (hx * hx);
415: col[3].i = ex + 1;
416: col[3].j = ey;
417: col[3].k = ez;
418: col[3].loc = LEFT;
419: col[3].c = 0;
420: valA[3] = 1.0 / (hx * hx);
421: col[4].i = ex;
422: col[4].j = ey;
423: col[4].k = ez - 1;
424: col[4].loc = LEFT;
425: col[4].c = 0;
426: valA[4] = 1.0 / (hz * hz);
427: col[5].i = ex;
428: col[5].j = ey;
429: col[5].k = ez + 1;
430: col[5].loc = LEFT;
431: col[5].c = 0;
432: valA[5] = 1.0 / (hz * hz);
433: col[6].i = ex - 1;
434: col[6].j = ey;
435: col[6].k = ez;
436: col[6].loc = ELEMENT;
437: col[6].c = 0;
438: valA[6] = 1.0 / hx;
439: col[7].i = ex;
440: col[7].j = ey;
441: col[7].k = ez;
442: col[7].loc = ELEMENT;
443: col[7].c = 0;
444: valA[7] = -1.0 / hx;
445: }
446: } else if (ez == 0) {
447: nEntries = 8;
448: col[0].i = ex;
449: col[0].j = ey;
450: col[0].k = ez;
451: col[0].loc = LEFT;
452: col[0].c = 0;
453: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
454: col[1].i = ex;
455: col[1].j = ey - 1;
456: col[1].k = ez;
457: col[1].loc = LEFT;
458: col[1].c = 0;
459: valA[1] = 1.0 / (hy * hy);
460: col[2].i = ex;
461: col[2].j = ey + 1;
462: col[2].k = ez;
463: col[2].loc = LEFT;
464: col[2].c = 0;
465: valA[2] = 1.0 / (hy * hy);
466: col[3].i = ex - 1;
467: col[3].j = ey;
468: col[3].k = ez;
469: col[3].loc = LEFT;
470: col[3].c = 0;
471: valA[3] = 1.0 / (hx * hx);
472: col[4].i = ex + 1;
473: col[4].j = ey;
474: col[4].k = ez;
475: col[4].loc = LEFT;
476: col[4].c = 0;
477: valA[4] = 1.0 / (hx * hx);
478: /* Missing back term */
479: col[5].i = ex;
480: col[5].j = ey;
481: col[5].k = ez + 1;
482: col[5].loc = LEFT;
483: col[5].c = 0;
484: valA[5] = 1.0 / (hz * hz);
485: col[6].i = ex - 1;
486: col[6].j = ey;
487: col[6].k = ez;
488: col[6].loc = ELEMENT;
489: col[6].c = 0;
490: valA[6] = 1.0 / hx;
491: col[7].i = ex;
492: col[7].j = ey;
493: col[7].k = ez;
494: col[7].loc = ELEMENT;
495: col[7].c = 0;
496: valA[7] = -1.0 / hx;
497: } else if (ez == N[2] - 1) {
498: nEntries = 8;
499: col[0].i = ex;
500: col[0].j = ey;
501: col[0].k = ez;
502: col[0].loc = LEFT;
503: col[0].c = 0;
504: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
505: col[1].i = ex;
506: col[1].j = ey - 1;
507: col[1].k = ez;
508: col[1].loc = LEFT;
509: col[1].c = 0;
510: valA[1] = 1.0 / (hy * hy);
511: col[2].i = ex;
512: col[2].j = ey + 1;
513: col[2].k = ez;
514: col[2].loc = LEFT;
515: col[2].c = 0;
516: valA[2] = 1.0 / (hy * hy);
517: col[3].i = ex - 1;
518: col[3].j = ey;
519: col[3].k = ez;
520: col[3].loc = LEFT;
521: col[3].c = 0;
522: valA[3] = 1.0 / (hx * hx);
523: col[4].i = ex + 1;
524: col[4].j = ey;
525: col[4].k = ez;
526: col[4].loc = LEFT;
527: col[4].c = 0;
528: valA[4] = 1.0 / (hx * hx);
529: col[5].i = ex;
530: col[5].j = ey;
531: col[5].k = ez - 1;
532: col[5].loc = LEFT;
533: col[5].c = 0;
534: valA[5] = 1.0 / (hz * hz);
535: /* Missing front term */
536: col[6].i = ex - 1;
537: col[6].j = ey;
538: col[6].k = ez;
539: col[6].loc = ELEMENT;
540: col[6].c = 0;
541: valA[6] = 1.0 / hx;
542: col[7].i = ex;
543: col[7].j = ey;
544: col[7].k = ez;
545: col[7].loc = ELEMENT;
546: col[7].c = 0;
547: valA[7] = -1.0 / hx;
548: } else {
549: nEntries = 9;
550: col[0].i = ex;
551: col[0].j = ey;
552: col[0].k = ez;
553: col[0].loc = LEFT;
554: col[0].c = 0;
555: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
556: col[1].i = ex;
557: col[1].j = ey - 1;
558: col[1].k = ez;
559: col[1].loc = LEFT;
560: col[1].c = 0;
561: valA[1] = 1.0 / (hy * hy);
562: col[2].i = ex;
563: col[2].j = ey + 1;
564: col[2].k = ez;
565: col[2].loc = LEFT;
566: col[2].c = 0;
567: valA[2] = 1.0 / (hy * hy);
568: col[3].i = ex - 1;
569: col[3].j = ey;
570: col[3].k = ez;
571: col[3].loc = LEFT;
572: col[3].c = 0;
573: valA[3] = 1.0 / (hx * hx);
574: col[4].i = ex + 1;
575: col[4].j = ey;
576: col[4].k = ez;
577: col[4].loc = LEFT;
578: col[4].c = 0;
579: valA[4] = 1.0 / (hx * hx);
580: col[5].i = ex;
581: col[5].j = ey;
582: col[5].k = ez - 1;
583: col[5].loc = LEFT;
584: col[5].c = 0;
585: valA[5] = 1.0 / (hz * hz);
586: col[6].i = ex;
587: col[6].j = ey;
588: col[6].k = ez + 1;
589: col[6].loc = LEFT;
590: col[6].c = 0;
591: valA[6] = 1.0 / (hz * hz);
592: col[7].i = ex - 1;
593: col[7].j = ey;
594: col[7].k = ez;
595: col[7].loc = ELEMENT;
596: col[7].c = 0;
597: valA[7] = 1.0 / hx;
598: col[8].i = ex;
599: col[8].j = ey;
600: col[8].k = ez;
601: col[8].loc = ELEMENT;
602: col[8].c = 0;
603: valA[8] = -1.0 / hx;
604: }
605: DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES);
606: }
608: /* Equation on bottom face of this element */
609: if (ey == 0) {
610: /* Bottom boundary velocity Dirichlet */
611: DMStagStencil row;
612: const PetscScalar valA = 1.0;
613: row.i = ex;
614: row.j = ey;
615: row.k = ez;
616: row.loc = DOWN;
617: row.c = 0;
618: DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
619: } else {
620: /* Y-momentum equation, (v_xx + v_yy + v_zz) - p_y = f^y */
621: DMStagStencil row, col[9];
622: PetscScalar valA[9];
623: PetscInt nEntries;
625: row.i = ex;
626: row.j = ey;
627: row.k = ez;
628: row.loc = DOWN;
629: row.c = 0;
630: if (ex == 0) {
631: if (ez == 0) {
632: nEntries = 7;
633: col[0].i = ex;
634: col[0].j = ey;
635: col[0].k = ez;
636: col[0].loc = DOWN;
637: col[0].c = 0;
638: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
639: col[1].i = ex;
640: col[1].j = ey - 1;
641: col[1].k = ez;
642: col[1].loc = DOWN;
643: col[1].c = 0;
644: valA[1] = 1.0 / (hy * hy);
645: col[2].i = ex;
646: col[2].j = ey + 1;
647: col[2].k = ez;
648: col[2].loc = DOWN;
649: col[2].c = 0;
650: valA[2] = 1.0 / (hy * hy);
651: /* Left term missing */
652: col[3].i = ex + 1;
653: col[3].j = ey;
654: col[3].k = ez;
655: col[3].loc = DOWN;
656: col[3].c = 0;
657: valA[3] = 1.0 / (hx * hx);
658: /* Back term missing */
659: col[4].i = ex;
660: col[4].j = ey;
661: col[4].k = ez + 1;
662: col[4].loc = DOWN;
663: col[4].c = 0;
664: valA[4] = 1.0 / (hz * hz);
665: col[5].i = ex;
666: col[5].j = ey - 1;
667: col[5].k = ez;
668: col[5].loc = ELEMENT;
669: col[5].c = 0;
670: valA[5] = 1.0 / hy;
671: col[6].i = ex;
672: col[6].j = ey;
673: col[6].k = ez;
674: col[6].loc = ELEMENT;
675: col[6].c = 0;
676: valA[6] = -1.0 / hy;
677: } else if (ez == N[2] - 1) {
678: nEntries = 7;
679: col[0].i = ex;
680: col[0].j = ey;
681: col[0].k = ez;
682: col[0].loc = DOWN;
683: col[0].c = 0;
684: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
685: col[1].i = ex;
686: col[1].j = ey - 1;
687: col[1].k = ez;
688: col[1].loc = DOWN;
689: col[1].c = 0;
690: valA[1] = 1.0 / (hy * hy);
691: col[2].i = ex;
692: col[2].j = ey + 1;
693: col[2].k = ez;
694: col[2].loc = DOWN;
695: col[2].c = 0;
696: valA[2] = 1.0 / (hy * hy);
697: /* Left term missing */
698: col[3].i = ex + 1;
699: col[3].j = ey;
700: col[3].k = ez;
701: col[3].loc = DOWN;
702: col[3].c = 0;
703: valA[3] = 1.0 / (hx * hx);
704: col[4].i = ex;
705: col[4].j = ey;
706: col[4].k = ez - 1;
707: col[4].loc = DOWN;
708: col[4].c = 0;
709: valA[4] = 1.0 / (hz * hz);
710: /* Front term missing */
711: col[5].i = ex;
712: col[5].j = ey - 1;
713: col[5].k = ez;
714: col[5].loc = ELEMENT;
715: col[5].c = 0;
716: valA[5] = 1.0 / hy;
717: col[6].i = ex;
718: col[6].j = ey;
719: col[6].k = ez;
720: col[6].loc = ELEMENT;
721: col[6].c = 0;
722: valA[6] = -1.0 / hy;
723: } else {
724: nEntries = 8;
725: col[0].i = ex;
726: col[0].j = ey;
727: col[0].k = ez;
728: col[0].loc = DOWN;
729: col[0].c = 0;
730: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
731: col[1].i = ex;
732: col[1].j = ey - 1;
733: col[1].k = ez;
734: col[1].loc = DOWN;
735: col[1].c = 0;
736: valA[1] = 1.0 / (hy * hy);
737: col[2].i = ex;
738: col[2].j = ey + 1;
739: col[2].k = ez;
740: col[2].loc = DOWN;
741: col[2].c = 0;
742: valA[2] = 1.0 / (hy * hy);
743: /* Left term missing */
744: col[3].i = ex + 1;
745: col[3].j = ey;
746: col[3].k = ez;
747: col[3].loc = DOWN;
748: col[3].c = 0;
749: valA[3] = 1.0 / (hx * hx);
750: col[4].i = ex;
751: col[4].j = ey;
752: col[4].k = ez - 1;
753: col[4].loc = DOWN;
754: col[4].c = 0;
755: valA[4] = 1.0 / (hz * hz);
756: col[5].i = ex;
757: col[5].j = ey;
758: col[5].k = ez + 1;
759: col[5].loc = DOWN;
760: col[5].c = 0;
761: valA[5] = 1.0 / (hz * hz);
762: col[6].i = ex;
763: col[6].j = ey - 1;
764: col[6].k = ez;
765: col[6].loc = ELEMENT;
766: col[6].c = 0;
767: valA[6] = 1.0 / hy;
768: col[7].i = ex;
769: col[7].j = ey;
770: col[7].k = ez;
771: col[7].loc = ELEMENT;
772: col[7].c = 0;
773: valA[7] = -1.0 / hy;
774: }
775: } else if (ex == N[0] - 1) {
776: if (ez == 0) {
777: nEntries = 7;
778: col[0].i = ex;
779: col[0].j = ey;
780: col[0].k = ez;
781: col[0].loc = DOWN;
782: col[0].c = 0;
783: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
784: col[1].i = ex;
785: col[1].j = ey - 1;
786: col[1].k = ez;
787: col[1].loc = DOWN;
788: col[1].c = 0;
789: valA[1] = 1.0 / (hy * hy);
790: col[2].i = ex;
791: col[2].j = ey + 1;
792: col[2].k = ez;
793: col[2].loc = DOWN;
794: col[2].c = 0;
795: valA[2] = 1.0 / (hy * hy);
796: col[3].i = ex - 1;
797: col[3].j = ey;
798: col[3].k = ez;
799: col[3].loc = DOWN;
800: col[3].c = 0;
801: valA[3] = 1.0 / (hx * hx);
802: /* Right term missing */
803: /* Back term missing */
804: col[4].i = ex;
805: col[4].j = ey;
806: col[4].k = ez + 1;
807: col[4].loc = DOWN;
808: col[4].c = 0;
809: valA[4] = 1.0 / (hz * hz);
810: col[5].i = ex;
811: col[5].j = ey - 1;
812: col[5].k = ez;
813: col[5].loc = ELEMENT;
814: col[5].c = 0;
815: valA[5] = 1.0 / hy;
816: col[6].i = ex;
817: col[6].j = ey;
818: col[6].k = ez;
819: col[6].loc = ELEMENT;
820: col[6].c = 0;
821: valA[6] = -1.0 / hy;
822: } else if (ez == N[2] - 1) {
823: nEntries = 7;
824: col[0].i = ex;
825: col[0].j = ey;
826: col[0].k = ez;
827: col[0].loc = DOWN;
828: col[0].c = 0;
829: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
830: col[1].i = ex;
831: col[1].j = ey - 1;
832: col[1].k = ez;
833: col[1].loc = DOWN;
834: col[1].c = 0;
835: valA[1] = 1.0 / (hy * hy);
836: col[2].i = ex;
837: col[2].j = ey + 1;
838: col[2].k = ez;
839: col[2].loc = DOWN;
840: col[2].c = 0;
841: valA[2] = 1.0 / (hy * hy);
842: col[3].i = ex - 1;
843: col[3].j = ey;
844: col[3].k = ez;
845: col[3].loc = DOWN;
846: col[3].c = 0;
847: valA[3] = 1.0 / (hx * hx);
848: /* Right term missing */
849: col[4].i = ex;
850: col[4].j = ey;
851: col[4].k = ez - 1;
852: col[4].loc = DOWN;
853: col[4].c = 0;
854: valA[4] = 1.0 / (hz * hz);
855: /* Front term missing */
856: col[5].i = ex;
857: col[5].j = ey - 1;
858: col[5].k = ez;
859: col[5].loc = ELEMENT;
860: col[5].c = 0;
861: valA[5] = 1.0 / hy;
862: col[6].i = ex;
863: col[6].j = ey;
864: col[6].k = ez;
865: col[6].loc = ELEMENT;
866: col[6].c = 0;
867: valA[6] = -1.0 / hy;
868: } else {
869: nEntries = 8;
870: col[0].i = ex;
871: col[0].j = ey;
872: col[0].k = ez;
873: col[0].loc = DOWN;
874: col[0].c = 0;
875: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
876: col[1].i = ex;
877: col[1].j = ey - 1;
878: col[1].k = ez;
879: col[1].loc = DOWN;
880: col[1].c = 0;
881: valA[1] = 1.0 / (hy * hy);
882: col[2].i = ex;
883: col[2].j = ey + 1;
884: col[2].k = ez;
885: col[2].loc = DOWN;
886: col[2].c = 0;
887: valA[2] = 1.0 / (hy * hy);
888: col[3].i = ex - 1;
889: col[3].j = ey;
890: col[3].k = ez;
891: col[3].loc = DOWN;
892: col[3].c = 0;
893: valA[3] = 1.0 / (hx * hx);
894: /* Right term missing */
895: col[4].i = ex;
896: col[4].j = ey;
897: col[4].k = ez - 1;
898: col[4].loc = DOWN;
899: col[4].c = 0;
900: valA[4] = 1.0 / (hz * hz);
901: col[5].i = ex;
902: col[5].j = ey;
903: col[5].k = ez + 1;
904: col[5].loc = DOWN;
905: col[5].c = 0;
906: valA[5] = 1.0 / (hz * hz);
907: col[6].i = ex;
908: col[6].j = ey - 1;
909: col[6].k = ez;
910: col[6].loc = ELEMENT;
911: col[6].c = 0;
912: valA[6] = 1.0 / hy;
913: col[7].i = ex;
914: col[7].j = ey;
915: col[7].k = ez;
916: col[7].loc = ELEMENT;
917: col[7].c = 0;
918: valA[7] = -1.0 / hy;
919: }
920: } else if (ez == 0) {
921: nEntries = 8;
922: col[0].i = ex;
923: col[0].j = ey;
924: col[0].k = ez;
925: col[0].loc = DOWN;
926: col[0].c = 0;
927: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
928: col[1].i = ex;
929: col[1].j = ey - 1;
930: col[1].k = ez;
931: col[1].loc = DOWN;
932: col[1].c = 0;
933: valA[1] = 1.0 / (hy * hy);
934: col[2].i = ex;
935: col[2].j = ey + 1;
936: col[2].k = ez;
937: col[2].loc = DOWN;
938: col[2].c = 0;
939: valA[2] = 1.0 / (hy * hy);
940: col[3].i = ex - 1;
941: col[3].j = ey;
942: col[3].k = ez;
943: col[3].loc = DOWN;
944: col[3].c = 0;
945: valA[3] = 1.0 / (hx * hx);
946: col[4].i = ex + 1;
947: col[4].j = ey;
948: col[4].k = ez;
949: col[4].loc = DOWN;
950: col[4].c = 0;
951: valA[4] = 1.0 / (hx * hx);
952: /* Back term missing */
953: col[5].i = ex;
954: col[5].j = ey;
955: col[5].k = ez + 1;
956: col[5].loc = DOWN;
957: col[5].c = 0;
958: valA[5] = 1.0 / (hz * hz);
959: col[6].i = ex;
960: col[6].j = ey - 1;
961: col[6].k = ez;
962: col[6].loc = ELEMENT;
963: col[6].c = 0;
964: valA[6] = 1.0 / hy;
965: col[7].i = ex;
966: col[7].j = ey;
967: col[7].k = ez;
968: col[7].loc = ELEMENT;
969: col[7].c = 0;
970: valA[7] = -1.0 / hy;
971: } else if (ez == N[2] - 1) {
972: nEntries = 8;
973: col[0].i = ex;
974: col[0].j = ey;
975: col[0].k = ez;
976: col[0].loc = DOWN;
977: col[0].c = 0;
978: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
979: col[1].i = ex;
980: col[1].j = ey - 1;
981: col[1].k = ez;
982: col[1].loc = DOWN;
983: col[1].c = 0;
984: valA[1] = 1.0 / (hy * hy);
985: col[2].i = ex;
986: col[2].j = ey + 1;
987: col[2].k = ez;
988: col[2].loc = DOWN;
989: col[2].c = 0;
990: valA[2] = 1.0 / (hy * hy);
991: col[3].i = ex - 1;
992: col[3].j = ey;
993: col[3].k = ez;
994: col[3].loc = DOWN;
995: col[3].c = 0;
996: valA[3] = 1.0 / (hx * hx);
997: col[4].i = ex + 1;
998: col[4].j = ey;
999: col[4].k = ez;
1000: col[4].loc = DOWN;
1001: col[4].c = 0;
1002: valA[4] = 1.0 / (hx * hx);
1003: col[5].i = ex;
1004: col[5].j = ey;
1005: col[5].k = ez - 1;
1006: col[5].loc = DOWN;
1007: col[5].c = 0;
1008: valA[5] = 1.0 / (hz * hz);
1009: /* Front term missing */
1010: col[6].i = ex;
1011: col[6].j = ey - 1;
1012: col[6].k = ez;
1013: col[6].loc = ELEMENT;
1014: col[6].c = 0;
1015: valA[6] = 1.0 / hy;
1016: col[7].i = ex;
1017: col[7].j = ey;
1018: col[7].k = ez;
1019: col[7].loc = ELEMENT;
1020: col[7].c = 0;
1021: valA[7] = -1.0 / hy;
1022: } else {
1023: nEntries = 9;
1024: col[0].i = ex;
1025: col[0].j = ey;
1026: col[0].k = ez;
1027: col[0].loc = DOWN;
1028: col[0].c = 0;
1029: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1030: col[1].i = ex;
1031: col[1].j = ey - 1;
1032: col[1].k = ez;
1033: col[1].loc = DOWN;
1034: col[1].c = 0;
1035: valA[1] = 1.0 / (hy * hy);
1036: col[2].i = ex;
1037: col[2].j = ey + 1;
1038: col[2].k = ez;
1039: col[2].loc = DOWN;
1040: col[2].c = 0;
1041: valA[2] = 1.0 / (hy * hy);
1042: col[3].i = ex - 1;
1043: col[3].j = ey;
1044: col[3].k = ez;
1045: col[3].loc = DOWN;
1046: col[3].c = 0;
1047: valA[3] = 1.0 / (hx * hx);
1048: col[4].i = ex + 1;
1049: col[4].j = ey;
1050: col[4].k = ez;
1051: col[4].loc = DOWN;
1052: col[4].c = 0;
1053: valA[4] = 1.0 / (hx * hx);
1054: col[5].i = ex;
1055: col[5].j = ey;
1056: col[5].k = ez - 1;
1057: col[5].loc = DOWN;
1058: col[5].c = 0;
1059: valA[5] = 1.0 / (hz * hz);
1060: col[6].i = ex;
1061: col[6].j = ey;
1062: col[6].k = ez + 1;
1063: col[6].loc = DOWN;
1064: col[6].c = 0;
1065: valA[6] = 1.0 / (hz * hz);
1066: col[7].i = ex;
1067: col[7].j = ey - 1;
1068: col[7].k = ez;
1069: col[7].loc = ELEMENT;
1070: col[7].c = 0;
1071: valA[7] = 1.0 / hy;
1072: col[8].i = ex;
1073: col[8].j = ey;
1074: col[8].k = ez;
1075: col[8].loc = ELEMENT;
1076: col[8].c = 0;
1077: valA[8] = -1.0 / hy;
1078: }
1079: DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES);
1080: }
1082: /* Equation on back face of this element */
1083: if (ez == 0) {
1084: /* Back boundary velocity Dirichlet */
1085: DMStagStencil row;
1086: const PetscScalar valA = 1.0;
1087: row.i = ex;
1088: row.j = ey;
1089: row.k = ez;
1090: row.loc = BACK;
1091: row.c = 0;
1092: DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
1093: } else {
1094: /* Z-momentum equation, (w_xx + w_yy + w_zz) - p_z = f^z */
1095: DMStagStencil row, col[9];
1096: PetscScalar valA[9];
1097: PetscInt nEntries;
1099: row.i = ex;
1100: row.j = ey;
1101: row.k = ez;
1102: row.loc = BACK;
1103: row.c = 0;
1104: if (ex == 0) {
1105: if (ey == 0) {
1106: nEntries = 7;
1107: col[0].i = ex;
1108: col[0].j = ey;
1109: col[0].k = ez;
1110: col[0].loc = BACK;
1111: col[0].c = 0;
1112: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1113: /* Down term missing */
1114: col[1].i = ex;
1115: col[1].j = ey + 1;
1116: col[1].k = ez;
1117: col[1].loc = BACK;
1118: col[1].c = 0;
1119: valA[1] = 1.0 / (hy * hy);
1120: /* Left term missing */
1121: col[2].i = ex + 1;
1122: col[2].j = ey;
1123: col[2].k = ez;
1124: col[2].loc = BACK;
1125: col[2].c = 0;
1126: valA[2] = 1.0 / (hx * hx);
1127: col[3].i = ex;
1128: col[3].j = ey;
1129: col[3].k = ez - 1;
1130: col[3].loc = BACK;
1131: col[3].c = 0;
1132: valA[3] = 1.0 / (hz * hz);
1133: col[4].i = ex;
1134: col[4].j = ey;
1135: col[4].k = ez + 1;
1136: col[4].loc = BACK;
1137: col[4].c = 0;
1138: valA[4] = 1.0 / (hz * hz);
1139: col[5].i = ex;
1140: col[5].j = ey;
1141: col[5].k = ez - 1;
1142: col[5].loc = ELEMENT;
1143: col[5].c = 0;
1144: valA[5] = 1.0 / hz;
1145: col[6].i = ex;
1146: col[6].j = ey;
1147: col[6].k = ez;
1148: col[6].loc = ELEMENT;
1149: col[6].c = 0;
1150: valA[6] = -1.0 / hz;
1151: } else if (ey == N[1] - 1) {
1152: nEntries = 7;
1153: col[0].i = ex;
1154: col[0].j = ey;
1155: col[0].k = ez;
1156: col[0].loc = BACK;
1157: col[0].c = 0;
1158: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1159: col[1].i = ex;
1160: col[1].j = ey - 1;
1161: col[1].k = ez;
1162: col[1].loc = BACK;
1163: col[1].c = 0;
1164: valA[1] = 1.0 / (hy * hy);
1165: /* Up term missing */
1166: /* Left term missing */
1167: col[2].i = ex + 1;
1168: col[2].j = ey;
1169: col[2].k = ez;
1170: col[2].loc = BACK;
1171: col[2].c = 0;
1172: valA[2] = 1.0 / (hx * hx);
1173: col[3].i = ex;
1174: col[3].j = ey;
1175: col[3].k = ez - 1;
1176: col[3].loc = BACK;
1177: col[3].c = 0;
1178: valA[3] = 1.0 / (hz * hz);
1179: col[4].i = ex;
1180: col[4].j = ey;
1181: col[4].k = ez + 1;
1182: col[4].loc = BACK;
1183: col[4].c = 0;
1184: valA[4] = 1.0 / (hz * hz);
1185: col[5].i = ex;
1186: col[5].j = ey;
1187: col[5].k = ez - 1;
1188: col[5].loc = ELEMENT;
1189: col[5].c = 0;
1190: valA[5] = 1.0 / hz;
1191: col[6].i = ex;
1192: col[6].j = ey;
1193: col[6].k = ez;
1194: col[6].loc = ELEMENT;
1195: col[6].c = 0;
1196: valA[6] = -1.0 / hz;
1197: } else {
1198: nEntries = 8;
1199: col[0].i = ex;
1200: col[0].j = ey;
1201: col[0].k = ez;
1202: col[0].loc = BACK;
1203: col[0].c = 0;
1204: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1205: col[1].i = ex;
1206: col[1].j = ey - 1;
1207: col[1].k = ez;
1208: col[1].loc = BACK;
1209: col[1].c = 0;
1210: valA[1] = 1.0 / (hy * hy);
1211: col[2].i = ex;
1212: col[2].j = ey + 1;
1213: col[2].k = ez;
1214: col[2].loc = BACK;
1215: col[2].c = 0;
1216: valA[2] = 1.0 / (hy * hy);
1217: /* Left term missing */
1218: col[3].i = ex + 1;
1219: col[3].j = ey;
1220: col[3].k = ez;
1221: col[3].loc = BACK;
1222: col[3].c = 0;
1223: valA[3] = 1.0 / (hx * hx);
1224: col[4].i = ex;
1225: col[4].j = ey;
1226: col[4].k = ez - 1;
1227: col[4].loc = BACK;
1228: col[4].c = 0;
1229: valA[4] = 1.0 / (hz * hz);
1230: col[5].i = ex;
1231: col[5].j = ey;
1232: col[5].k = ez + 1;
1233: col[5].loc = BACK;
1234: col[5].c = 0;
1235: valA[5] = 1.0 / (hz * hz);
1236: col[6].i = ex;
1237: col[6].j = ey;
1238: col[6].k = ez - 1;
1239: col[6].loc = ELEMENT;
1240: col[6].c = 0;
1241: valA[6] = 1.0 / hz;
1242: col[7].i = ex;
1243: col[7].j = ey;
1244: col[7].k = ez;
1245: col[7].loc = ELEMENT;
1246: col[7].c = 0;
1247: valA[7] = -1.0 / hz;
1248: }
1249: } else if (ex == N[0] - 1) {
1250: if (ey == 0) {
1251: nEntries = 7;
1252: col[0].i = ex;
1253: col[0].j = ey;
1254: col[0].k = ez;
1255: col[0].loc = BACK;
1256: col[0].c = 0;
1257: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1258: /* Down term missing */
1259: col[1].i = ex;
1260: col[1].j = ey + 1;
1261: col[1].k = ez;
1262: col[1].loc = BACK;
1263: col[1].c = 0;
1264: valA[1] = 1.0 / (hy * hy);
1265: col[2].i = ex - 1;
1266: col[2].j = ey;
1267: col[2].k = ez;
1268: col[2].loc = BACK;
1269: col[2].c = 0;
1270: valA[2] = 1.0 / (hx * hx);
1271: /* Right term missing */
1272: col[3].i = ex;
1273: col[3].j = ey;
1274: col[3].k = ez - 1;
1275: col[3].loc = BACK;
1276: col[3].c = 0;
1277: valA[3] = 1.0 / (hz * hz);
1278: col[4].i = ex;
1279: col[4].j = ey;
1280: col[4].k = ez + 1;
1281: col[4].loc = BACK;
1282: col[4].c = 0;
1283: valA[4] = 1.0 / (hz * hz);
1284: col[5].i = ex;
1285: col[5].j = ey;
1286: col[5].k = ez - 1;
1287: col[5].loc = ELEMENT;
1288: col[5].c = 0;
1289: valA[5] = 1.0 / hz;
1290: col[6].i = ex;
1291: col[6].j = ey;
1292: col[6].k = ez;
1293: col[6].loc = ELEMENT;
1294: col[6].c = 0;
1295: valA[6] = -1.0 / hz;
1296: } else if (ey == N[1] - 1) {
1297: nEntries = 7;
1298: col[0].i = ex;
1299: col[0].j = ey;
1300: col[0].k = ez;
1301: col[0].loc = BACK;
1302: col[0].c = 0;
1303: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1304: col[1].i = ex;
1305: col[1].j = ey - 1;
1306: col[1].k = ez;
1307: col[1].loc = BACK;
1308: col[1].c = 0;
1309: valA[1] = 1.0 / (hy * hy);
1310: /* Up term missing */
1311: col[2].i = ex - 1;
1312: col[2].j = ey;
1313: col[2].k = ez;
1314: col[2].loc = BACK;
1315: col[2].c = 0;
1316: valA[2] = 1.0 / (hx * hx);
1317: /* Right term missing */
1318: col[3].i = ex;
1319: col[3].j = ey;
1320: col[3].k = ez - 1;
1321: col[3].loc = BACK;
1322: col[3].c = 0;
1323: valA[3] = 1.0 / (hz * hz);
1324: col[4].i = ex;
1325: col[4].j = ey;
1326: col[4].k = ez + 1;
1327: col[4].loc = BACK;
1328: col[4].c = 0;
1329: valA[4] = 1.0 / (hz * hz);
1330: col[5].i = ex;
1331: col[5].j = ey;
1332: col[5].k = ez - 1;
1333: col[5].loc = ELEMENT;
1334: col[5].c = 0;
1335: valA[5] = 1.0 / hz;
1336: col[6].i = ex;
1337: col[6].j = ey;
1338: col[6].k = ez;
1339: col[6].loc = ELEMENT;
1340: col[6].c = 0;
1341: valA[6] = -1.0 / hz;
1342: } else {
1343: nEntries = 8;
1344: col[0].i = ex;
1345: col[0].j = ey;
1346: col[0].k = ez;
1347: col[0].loc = BACK;
1348: col[0].c = 0;
1349: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1350: col[1].i = ex;
1351: col[1].j = ey - 1;
1352: col[1].k = ez;
1353: col[1].loc = BACK;
1354: col[1].c = 0;
1355: valA[1] = 1.0 / (hy * hy);
1356: col[2].i = ex;
1357: col[2].j = ey + 1;
1358: col[2].k = ez;
1359: col[2].loc = BACK;
1360: col[2].c = 0;
1361: valA[2] = 1.0 / (hy * hy);
1362: col[3].i = ex - 1;
1363: col[3].j = ey;
1364: col[3].k = ez;
1365: col[3].loc = BACK;
1366: col[3].c = 0;
1367: valA[3] = 1.0 / (hx * hx);
1368: /* Right term missing */
1369: col[4].i = ex;
1370: col[4].j = ey;
1371: col[4].k = ez - 1;
1372: col[4].loc = BACK;
1373: col[4].c = 0;
1374: valA[4] = 1.0 / (hz * hz);
1375: col[5].i = ex;
1376: col[5].j = ey;
1377: col[5].k = ez + 1;
1378: col[5].loc = BACK;
1379: col[5].c = 0;
1380: valA[5] = 1.0 / (hz * hz);
1381: col[6].i = ex;
1382: col[6].j = ey;
1383: col[6].k = ez - 1;
1384: col[6].loc = ELEMENT;
1385: col[6].c = 0;
1386: valA[6] = 1.0 / hz;
1387: col[7].i = ex;
1388: col[7].j = ey;
1389: col[7].k = ez;
1390: col[7].loc = ELEMENT;
1391: col[7].c = 0;
1392: valA[7] = -1.0 / hz;
1393: }
1394: } else if (ey == 0) {
1395: nEntries = 8;
1396: col[0].i = ex;
1397: col[0].j = ey;
1398: col[0].k = ez;
1399: col[0].loc = BACK;
1400: col[0].c = 0;
1401: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1402: /* Down term missing */
1403: col[1].i = ex;
1404: col[1].j = ey + 1;
1405: col[1].k = ez;
1406: col[1].loc = BACK;
1407: col[1].c = 0;
1408: valA[1] = 1.0 / (hy * hy);
1409: col[2].i = ex - 1;
1410: col[2].j = ey;
1411: col[2].k = ez;
1412: col[2].loc = BACK;
1413: col[2].c = 0;
1414: valA[2] = 1.0 / (hx * hx);
1415: col[3].i = ex + 1;
1416: col[3].j = ey;
1417: col[3].k = ez;
1418: col[3].loc = BACK;
1419: col[3].c = 0;
1420: valA[3] = 1.0 / (hx * hx);
1421: col[4].i = ex;
1422: col[4].j = ey;
1423: col[4].k = ez - 1;
1424: col[4].loc = BACK;
1425: col[4].c = 0;
1426: valA[4] = 1.0 / (hz * hz);
1427: col[5].i = ex;
1428: col[5].j = ey;
1429: col[5].k = ez + 1;
1430: col[5].loc = BACK;
1431: col[5].c = 0;
1432: valA[5] = 1.0 / (hz * hz);
1433: col[6].i = ex;
1434: col[6].j = ey;
1435: col[6].k = ez - 1;
1436: col[6].loc = ELEMENT;
1437: col[6].c = 0;
1438: valA[6] = 1.0 / hz;
1439: col[7].i = ex;
1440: col[7].j = ey;
1441: col[7].k = ez;
1442: col[7].loc = ELEMENT;
1443: col[7].c = 0;
1444: valA[7] = -1.0 / hz;
1445: } else if (ey == N[1] - 1) {
1446: nEntries = 8;
1447: col[0].i = ex;
1448: col[0].j = ey;
1449: col[0].k = ez;
1450: col[0].loc = BACK;
1451: col[0].c = 0;
1452: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1453: col[1].i = ex;
1454: col[1].j = ey - 1;
1455: col[1].k = ez;
1456: col[1].loc = BACK;
1457: col[1].c = 0;
1458: valA[1] = 1.0 / (hy * hy);
1459: /* Up term missing */
1460: col[2].i = ex - 1;
1461: col[2].j = ey;
1462: col[2].k = ez;
1463: col[2].loc = BACK;
1464: col[2].c = 0;
1465: valA[2] = 1.0 / (hx * hx);
1466: col[3].i = ex + 1;
1467: col[3].j = ey;
1468: col[3].k = ez;
1469: col[3].loc = BACK;
1470: col[3].c = 0;
1471: valA[3] = 1.0 / (hx * hx);
1472: col[4].i = ex;
1473: col[4].j = ey;
1474: col[4].k = ez - 1;
1475: col[4].loc = BACK;
1476: col[4].c = 0;
1477: valA[4] = 1.0 / (hz * hz);
1478: col[5].i = ex;
1479: col[5].j = ey;
1480: col[5].k = ez + 1;
1481: col[5].loc = BACK;
1482: col[5].c = 0;
1483: valA[5] = 1.0 / (hz * hz);
1484: col[6].i = ex;
1485: col[6].j = ey;
1486: col[6].k = ez - 1;
1487: col[6].loc = ELEMENT;
1488: col[6].c = 0;
1489: valA[6] = 1.0 / hz;
1490: col[7].i = ex;
1491: col[7].j = ey;
1492: col[7].k = ez;
1493: col[7].loc = ELEMENT;
1494: col[7].c = 0;
1495: valA[7] = -1.0 / hz;
1496: } else {
1497: nEntries = 9;
1498: col[0].i = ex;
1499: col[0].j = ey;
1500: col[0].k = ez;
1501: col[0].loc = BACK;
1502: col[0].c = 0;
1503: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1504: col[1].i = ex;
1505: col[1].j = ey - 1;
1506: col[1].k = ez;
1507: col[1].loc = BACK;
1508: col[1].c = 0;
1509: valA[1] = 1.0 / (hy * hy);
1510: col[2].i = ex;
1511: col[2].j = ey + 1;
1512: col[2].k = ez;
1513: col[2].loc = BACK;
1514: col[2].c = 0;
1515: valA[2] = 1.0 / (hy * hy);
1516: col[3].i = ex - 1;
1517: col[3].j = ey;
1518: col[3].k = ez;
1519: col[3].loc = BACK;
1520: col[3].c = 0;
1521: valA[3] = 1.0 / (hx * hx);
1522: col[4].i = ex + 1;
1523: col[4].j = ey;
1524: col[4].k = ez;
1525: col[4].loc = BACK;
1526: col[4].c = 0;
1527: valA[4] = 1.0 / (hx * hx);
1528: col[5].i = ex;
1529: col[5].j = ey;
1530: col[5].k = ez - 1;
1531: col[5].loc = BACK;
1532: col[5].c = 0;
1533: valA[5] = 1.0 / (hz * hz);
1534: col[6].i = ex;
1535: col[6].j = ey;
1536: col[6].k = ez + 1;
1537: col[6].loc = BACK;
1538: col[6].c = 0;
1539: valA[6] = 1.0 / (hz * hz);
1540: col[7].i = ex;
1541: col[7].j = ey;
1542: col[7].k = ez - 1;
1543: col[7].loc = ELEMENT;
1544: col[7].c = 0;
1545: valA[7] = 1.0 / hz;
1546: col[8].i = ex;
1547: col[8].j = ey;
1548: col[8].k = ez;
1549: col[8].loc = ELEMENT;
1550: col[8].c = 0;
1551: valA[8] = -1.0 / hz;
1552: }
1553: DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES);
1554: }
1556: /* P equation : u_x + v_y + w_z = g
1557: Note that this includes an explicit zero on the diagonal. This is only needed for
1558: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
1559: {
1560: DMStagStencil row, col[7];
1561: PetscScalar valA[7];
1563: row.i = ex;
1564: row.j = ey;
1565: row.k = ez;
1566: row.loc = ELEMENT;
1567: row.c = 0;
1568: col[0].i = ex;
1569: col[0].j = ey;
1570: col[0].k = ez;
1571: col[0].loc = LEFT;
1572: col[0].c = 0;
1573: valA[0] = -1.0 / hx;
1574: col[1].i = ex;
1575: col[1].j = ey;
1576: col[1].k = ez;
1577: col[1].loc = RIGHT;
1578: col[1].c = 0;
1579: valA[1] = 1.0 / hx;
1580: col[2].i = ex;
1581: col[2].j = ey;
1582: col[2].k = ez;
1583: col[2].loc = DOWN;
1584: col[2].c = 0;
1585: valA[2] = -1.0 / hy;
1586: col[3].i = ex;
1587: col[3].j = ey;
1588: col[3].k = ez;
1589: col[3].loc = UP;
1590: col[3].c = 0;
1591: valA[3] = 1.0 / hy;
1592: col[4].i = ex;
1593: col[4].j = ey;
1594: col[4].k = ez;
1595: col[4].loc = BACK;
1596: col[4].c = 0;
1597: valA[4] = -1.0 / hz;
1598: col[5].i = ex;
1599: col[5].j = ey;
1600: col[5].k = ez;
1601: col[5].loc = FRONT;
1602: col[5].c = 0;
1603: valA[5] = 1.0 / hz;
1604: col[6] = row;
1605: valA[6] = 0.0;
1606: DMStagMatSetValuesStencil(dmSol, A, 1, &row, 7, col, valA, INSERT_VALUES);
1607: }
1608: }
1609: }
1610: }
1611: DMStagVecRestoreArrayRead(dmCoord, coordLocal, &arrCoord);
1612: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1613: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1615: return 0;
1616: }
1618: /* A helper function to check values */
1619: static PetscErrorCode check_vals(PetscInt ex, PetscInt ey, PetscInt ez, PetscInt n, const PetscScalar *ref, const PetscScalar *computed)
1620: {
1621: PetscInt i;
1624: for (i = 0; i < n; ++i) {
1626: }
1627: return 0;
1628: }
1630: /* The same function as above, but getting and checking values, instead of setting them */
1631: static PetscErrorCode CheckMat(DM dmSol, Mat A)
1632: {
1633: Vec coordLocal;
1634: PetscInt startx, starty, startz, N[3], nx, ny, nz, ex, ey, ez, d;
1635: PetscInt icp[3], icux[3], icuy[3], icuz[3], icux_right[3], icuy_up[3], icuz_front[3];
1636: PetscBool isLastRankx, isLastRanky, isLastRankz, isFirstRankx, isFirstRanky, isFirstRankz;
1637: PetscReal hx, hy, hz;
1638: DM dmCoord;
1639: PetscScalar ****arrCoord;
1640: PetscScalar computed[1024];
1643: DMStagGetCorners(dmSol, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL);
1644: DMStagGetGlobalSizes(dmSol, &N[0], &N[1], &N[2]);
1646: DMStagGetIsLastRank(dmSol, &isLastRankx, &isLastRanky, &isLastRankz);
1647: DMStagGetIsFirstRank(dmSol, &isFirstRankx, &isFirstRanky, &isFirstRankz);
1648: hx = 1.0 / N[0];
1649: hy = 1.0 / N[1];
1650: hz = 1.0 / N[2];
1651: DMGetCoordinateDM(dmSol, &dmCoord);
1652: DMGetCoordinatesLocal(dmSol, &coordLocal);
1653: DMStagVecGetArrayRead(dmCoord, coordLocal, &arrCoord);
1654: for (d = 0; d < 3; ++d) {
1655: DMStagGetLocationSlot(dmCoord, ELEMENT, d, &icp[d]);
1656: DMStagGetLocationSlot(dmCoord, LEFT, d, &icux[d]);
1657: DMStagGetLocationSlot(dmCoord, DOWN, d, &icuy[d]);
1658: DMStagGetLocationSlot(dmCoord, BACK, d, &icuz[d]);
1659: DMStagGetLocationSlot(dmCoord, RIGHT, d, &icux_right[d]);
1660: DMStagGetLocationSlot(dmCoord, UP, d, &icuy_up[d]);
1661: DMStagGetLocationSlot(dmCoord, FRONT, d, &icuz_front[d]);
1662: }
1664: for (ez = startz; ez < startz + nz; ++ez) {
1665: for (ey = starty; ey < starty + ny; ++ey) {
1666: for (ex = startx; ex < startx + nx; ++ex) {
1667: if (ex == N[0] - 1) {
1668: /* Right Boundary velocity Dirichlet */
1669: DMStagStencil row;
1670: const PetscScalar valA = 1.0;
1671: row.i = ex;
1672: row.j = ey;
1673: row.k = ez;
1674: row.loc = RIGHT;
1675: row.c = 0;
1676: DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed);
1677: check_vals(ex, ey, ez, 1, &valA, computed);
1678: }
1679: if (ey == N[1] - 1) {
1680: /* Top boundary velocity Dirichlet */
1681: DMStagStencil row;
1682: const PetscScalar valA = 1.0;
1683: row.i = ex;
1684: row.j = ey;
1685: row.k = ez;
1686: row.loc = UP;
1687: row.c = 0;
1688: DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed);
1689: check_vals(ex, ey, ez, 1, &valA, computed);
1690: }
1691: if (ez == N[2] - 1) {
1692: /* Top boundary velocity Dirichlet */
1693: DMStagStencil row;
1694: const PetscScalar valA = 1.0;
1695: row.i = ex;
1696: row.j = ey;
1697: row.k = ez;
1698: row.loc = FRONT;
1699: row.c = 0;
1700: DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed);
1701: check_vals(ex, ey, ez, 1, &valA, computed);
1702: }
1704: /* Equation on left face of this element */
1705: if (ex == 0) {
1706: /* Left velocity Dirichlet */
1707: DMStagStencil row;
1708: const PetscScalar valA = 1.0;
1709: row.i = ex;
1710: row.j = ey;
1711: row.k = ez;
1712: row.loc = LEFT;
1713: row.c = 0;
1714: DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed);
1715: check_vals(ex, ey, ez, 1, &valA, computed);
1716: } else {
1717: /* X-momentum interior equation : (u_xx + u_yy + u_zz) - p_x = f^x */
1718: DMStagStencil row, col[9];
1719: PetscScalar valA[9];
1720: PetscInt nEntries;
1722: row.i = ex;
1723: row.j = ey;
1724: row.k = ez;
1725: row.loc = LEFT;
1726: row.c = 0;
1727: if (ey == 0) {
1728: if (ez == 0) {
1729: nEntries = 7;
1730: col[0].i = ex;
1731: col[0].j = ey;
1732: col[0].k = ez;
1733: col[0].loc = LEFT;
1734: col[0].c = 0;
1735: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
1736: /* Missing down term */
1737: col[1].i = ex;
1738: col[1].j = ey + 1;
1739: col[1].k = ez;
1740: col[1].loc = LEFT;
1741: col[1].c = 0;
1742: valA[1] = 1.0 / (hy * hy);
1743: col[2].i = ex - 1;
1744: col[2].j = ey;
1745: col[2].k = ez;
1746: col[2].loc = LEFT;
1747: col[2].c = 0;
1748: valA[2] = 1.0 / (hx * hx);
1749: col[3].i = ex + 1;
1750: col[3].j = ey;
1751: col[3].k = ez;
1752: col[3].loc = LEFT;
1753: col[3].c = 0;
1754: valA[3] = 1.0 / (hx * hx);
1755: /* Missing back term */
1756: col[4].i = ex;
1757: col[4].j = ey;
1758: col[4].k = ez + 1;
1759: col[4].loc = LEFT;
1760: col[4].c = 0;
1761: valA[4] = 1.0 / (hz * hz);
1762: col[5].i = ex - 1;
1763: col[5].j = ey;
1764: col[5].k = ez;
1765: col[5].loc = ELEMENT;
1766: col[5].c = 0;
1767: valA[5] = 1.0 / hx;
1768: col[6].i = ex;
1769: col[6].j = ey;
1770: col[6].k = ez;
1771: col[6].loc = ELEMENT;
1772: col[6].c = 0;
1773: valA[6] = -1.0 / hx;
1774: } else if (ez == N[2] - 1) {
1775: nEntries = 7;
1776: col[0].i = ex;
1777: col[0].j = ey;
1778: col[0].k = ez;
1779: col[0].loc = LEFT;
1780: col[0].c = 0;
1781: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
1782: /* Missing down term */
1783: col[1].i = ex;
1784: col[1].j = ey + 1;
1785: col[1].k = ez;
1786: col[1].loc = LEFT;
1787: col[1].c = 0;
1788: valA[1] = 1.0 / (hy * hy);
1789: col[2].i = ex - 1;
1790: col[2].j = ey;
1791: col[2].k = ez;
1792: col[2].loc = LEFT;
1793: col[2].c = 0;
1794: valA[2] = 1.0 / (hx * hx);
1795: col[3].i = ex + 1;
1796: col[3].j = ey;
1797: col[3].k = ez;
1798: col[3].loc = LEFT;
1799: col[3].c = 0;
1800: valA[3] = 1.0 / (hx * hx);
1801: col[4].i = ex;
1802: col[4].j = ey;
1803: col[4].k = ez - 1;
1804: col[4].loc = LEFT;
1805: col[4].c = 0;
1806: valA[4] = 1.0 / (hz * hz);
1807: /* Missing front term */
1808: col[5].i = ex - 1;
1809: col[5].j = ey;
1810: col[5].k = ez;
1811: col[5].loc = ELEMENT;
1812: col[5].c = 0;
1813: valA[5] = 1.0 / hx;
1814: col[6].i = ex;
1815: col[6].j = ey;
1816: col[6].k = ez;
1817: col[6].loc = ELEMENT;
1818: col[6].c = 0;
1819: valA[6] = -1.0 / hx;
1820: } else {
1821: nEntries = 8;
1822: col[0].i = ex;
1823: col[0].j = ey;
1824: col[0].k = ez;
1825: col[0].loc = LEFT;
1826: col[0].c = 0;
1827: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1828: /* Missing down term */
1829: col[1].i = ex;
1830: col[1].j = ey + 1;
1831: col[1].k = ez;
1832: col[1].loc = LEFT;
1833: col[1].c = 0;
1834: valA[1] = 1.0 / (hy * hy);
1835: col[2].i = ex - 1;
1836: col[2].j = ey;
1837: col[2].k = ez;
1838: col[2].loc = LEFT;
1839: col[2].c = 0;
1840: valA[2] = 1.0 / (hx * hx);
1841: col[3].i = ex + 1;
1842: col[3].j = ey;
1843: col[3].k = ez;
1844: col[3].loc = LEFT;
1845: col[3].c = 0;
1846: valA[3] = 1.0 / (hx * hx);
1847: col[4].i = ex;
1848: col[4].j = ey;
1849: col[4].k = ez - 1;
1850: col[4].loc = LEFT;
1851: col[4].c = 0;
1852: valA[4] = 1.0 / (hz * hz);
1853: col[5].i = ex;
1854: col[5].j = ey;
1855: col[5].k = ez + 1;
1856: col[5].loc = LEFT;
1857: col[5].c = 0;
1858: valA[5] = 1.0 / (hz * hz);
1859: col[6].i = ex - 1;
1860: col[6].j = ey;
1861: col[6].k = ez;
1862: col[6].loc = ELEMENT;
1863: col[6].c = 0;
1864: valA[6] = 1.0 / hx;
1865: col[7].i = ex;
1866: col[7].j = ey;
1867: col[7].k = ez;
1868: col[7].loc = ELEMENT;
1869: col[7].c = 0;
1870: valA[7] = -1.0 / hx;
1871: }
1872: } else if (ey == N[1] - 1) {
1873: if (ez == 0) {
1874: nEntries = 7;
1875: col[0].i = ex;
1876: col[0].j = ey;
1877: col[0].k = ez;
1878: col[0].loc = LEFT;
1879: col[0].c = 0;
1880: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
1881: col[1].i = ex;
1882: col[1].j = ey - 1;
1883: col[1].k = ez;
1884: col[1].loc = LEFT;
1885: col[1].c = 0;
1886: valA[1] = 1.0 / (hy * hy);
1887: /* Missing up term */
1888: col[2].i = ex - 1;
1889: col[2].j = ey;
1890: col[2].k = ez;
1891: col[2].loc = LEFT;
1892: col[2].c = 0;
1893: valA[2] = 1.0 / (hx * hx);
1894: col[3].i = ex + 1;
1895: col[3].j = ey;
1896: col[3].k = ez;
1897: col[3].loc = LEFT;
1898: col[3].c = 0;
1899: valA[3] = 1.0 / (hx * hx);
1900: /* Missing back entry */
1901: col[4].i = ex;
1902: col[4].j = ey;
1903: col[4].k = ez + 1;
1904: col[4].loc = LEFT;
1905: col[4].c = 0;
1906: valA[4] = 1.0 / (hz * hz);
1907: col[5].i = ex - 1;
1908: col[5].j = ey;
1909: col[5].k = ez;
1910: col[5].loc = ELEMENT;
1911: col[5].c = 0;
1912: valA[5] = 1.0 / hx;
1913: col[6].i = ex;
1914: col[6].j = ey;
1915: col[6].k = ez;
1916: col[6].loc = ELEMENT;
1917: col[6].c = 0;
1918: valA[6] = -1.0 / hx;
1919: } else if (ez == N[2] - 1) {
1920: nEntries = 7;
1921: col[0].i = ex;
1922: col[0].j = ey;
1923: col[0].k = ez;
1924: col[0].loc = LEFT;
1925: col[0].c = 0;
1926: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
1927: col[1].i = ex;
1928: col[1].j = ey - 1;
1929: col[1].k = ez;
1930: col[1].loc = LEFT;
1931: col[1].c = 0;
1932: valA[1] = 1.0 / (hy * hy);
1933: /* Missing up term */
1934: col[2].i = ex - 1;
1935: col[2].j = ey;
1936: col[2].k = ez;
1937: col[2].loc = LEFT;
1938: col[2].c = 0;
1939: valA[2] = 1.0 / (hx * hx);
1940: col[3].i = ex + 1;
1941: col[3].j = ey;
1942: col[3].k = ez;
1943: col[3].loc = LEFT;
1944: col[3].c = 0;
1945: valA[3] = 1.0 / (hx * hx);
1946: col[4].i = ex;
1947: col[4].j = ey;
1948: col[4].k = ez - 1;
1949: col[4].loc = LEFT;
1950: col[4].c = 0;
1951: valA[4] = 1.0 / (hz * hz);
1952: /* Missing front term */
1953: col[5].i = ex - 1;
1954: col[5].j = ey;
1955: col[5].k = ez;
1956: col[5].loc = ELEMENT;
1957: col[5].c = 0;
1958: valA[5] = 1.0 / hx;
1959: col[6].i = ex;
1960: col[6].j = ey;
1961: col[6].k = ez;
1962: col[6].loc = ELEMENT;
1963: col[6].c = 0;
1964: valA[6] = -1.0 / hx;
1965: } else {
1966: nEntries = 8;
1967: col[0].i = ex;
1968: col[0].j = ey;
1969: col[0].k = ez;
1970: col[0].loc = LEFT;
1971: col[0].c = 0;
1972: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1973: col[1].i = ex;
1974: col[1].j = ey - 1;
1975: col[1].k = ez;
1976: col[1].loc = LEFT;
1977: col[1].c = 0;
1978: valA[1] = 1.0 / (hy * hy);
1979: /* Missing up term */
1980: col[2].i = ex - 1;
1981: col[2].j = ey;
1982: col[2].k = ez;
1983: col[2].loc = LEFT;
1984: col[2].c = 0;
1985: valA[2] = 1.0 / (hx * hx);
1986: col[3].i = ex + 1;
1987: col[3].j = ey;
1988: col[3].k = ez;
1989: col[3].loc = LEFT;
1990: col[3].c = 0;
1991: valA[3] = 1.0 / (hx * hx);
1992: col[4].i = ex;
1993: col[4].j = ey;
1994: col[4].k = ez - 1;
1995: col[4].loc = LEFT;
1996: col[4].c = 0;
1997: valA[4] = 1.0 / (hz * hz);
1998: col[5].i = ex;
1999: col[5].j = ey;
2000: col[5].k = ez + 1;
2001: col[5].loc = LEFT;
2002: col[5].c = 0;
2003: valA[5] = 1.0 / (hz * hz);
2004: col[6].i = ex - 1;
2005: col[6].j = ey;
2006: col[6].k = ez;
2007: col[6].loc = ELEMENT;
2008: col[6].c = 0;
2009: valA[6] = 1.0 / hx;
2010: col[7].i = ex;
2011: col[7].j = ey;
2012: col[7].k = ez;
2013: col[7].loc = ELEMENT;
2014: col[7].c = 0;
2015: valA[7] = -1.0 / hx;
2016: }
2017: } else if (ez == 0) {
2018: nEntries = 8;
2019: col[0].i = ex;
2020: col[0].j = ey;
2021: col[0].k = ez;
2022: col[0].loc = LEFT;
2023: col[0].c = 0;
2024: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2025: col[1].i = ex;
2026: col[1].j = ey - 1;
2027: col[1].k = ez;
2028: col[1].loc = LEFT;
2029: col[1].c = 0;
2030: valA[1] = 1.0 / (hy * hy);
2031: col[2].i = ex;
2032: col[2].j = ey + 1;
2033: col[2].k = ez;
2034: col[2].loc = LEFT;
2035: col[2].c = 0;
2036: valA[2] = 1.0 / (hy * hy);
2037: col[3].i = ex - 1;
2038: col[3].j = ey;
2039: col[3].k = ez;
2040: col[3].loc = LEFT;
2041: col[3].c = 0;
2042: valA[3] = 1.0 / (hx * hx);
2043: col[4].i = ex + 1;
2044: col[4].j = ey;
2045: col[4].k = ez;
2046: col[4].loc = LEFT;
2047: col[4].c = 0;
2048: valA[4] = 1.0 / (hx * hx);
2049: /* Missing back term */
2050: col[5].i = ex;
2051: col[5].j = ey;
2052: col[5].k = ez + 1;
2053: col[5].loc = LEFT;
2054: col[5].c = 0;
2055: valA[5] = 1.0 / (hz * hz);
2056: col[6].i = ex - 1;
2057: col[6].j = ey;
2058: col[6].k = ez;
2059: col[6].loc = ELEMENT;
2060: col[6].c = 0;
2061: valA[6] = 1.0 / hx;
2062: col[7].i = ex;
2063: col[7].j = ey;
2064: col[7].k = ez;
2065: col[7].loc = ELEMENT;
2066: col[7].c = 0;
2067: valA[7] = -1.0 / hx;
2068: } else if (ez == N[2] - 1) {
2069: nEntries = 8;
2070: col[0].i = ex;
2071: col[0].j = ey;
2072: col[0].k = ez;
2073: col[0].loc = LEFT;
2074: col[0].c = 0;
2075: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2076: col[1].i = ex;
2077: col[1].j = ey - 1;
2078: col[1].k = ez;
2079: col[1].loc = LEFT;
2080: col[1].c = 0;
2081: valA[1] = 1.0 / (hy * hy);
2082: col[2].i = ex;
2083: col[2].j = ey + 1;
2084: col[2].k = ez;
2085: col[2].loc = LEFT;
2086: col[2].c = 0;
2087: valA[2] = 1.0 / (hy * hy);
2088: col[3].i = ex - 1;
2089: col[3].j = ey;
2090: col[3].k = ez;
2091: col[3].loc = LEFT;
2092: col[3].c = 0;
2093: valA[3] = 1.0 / (hx * hx);
2094: col[4].i = ex + 1;
2095: col[4].j = ey;
2096: col[4].k = ez;
2097: col[4].loc = LEFT;
2098: col[4].c = 0;
2099: valA[4] = 1.0 / (hx * hx);
2100: col[5].i = ex;
2101: col[5].j = ey;
2102: col[5].k = ez - 1;
2103: col[5].loc = LEFT;
2104: col[5].c = 0;
2105: valA[5] = 1.0 / (hz * hz);
2106: /* Missing front term */
2107: col[6].i = ex - 1;
2108: col[6].j = ey;
2109: col[6].k = ez;
2110: col[6].loc = ELEMENT;
2111: col[6].c = 0;
2112: valA[6] = 1.0 / hx;
2113: col[7].i = ex;
2114: col[7].j = ey;
2115: col[7].k = ez;
2116: col[7].loc = ELEMENT;
2117: col[7].c = 0;
2118: valA[7] = -1.0 / hx;
2119: } else {
2120: nEntries = 9;
2121: col[0].i = ex;
2122: col[0].j = ey;
2123: col[0].k = ez;
2124: col[0].loc = LEFT;
2125: col[0].c = 0;
2126: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2127: col[1].i = ex;
2128: col[1].j = ey - 1;
2129: col[1].k = ez;
2130: col[1].loc = LEFT;
2131: col[1].c = 0;
2132: valA[1] = 1.0 / (hy * hy);
2133: col[2].i = ex;
2134: col[2].j = ey + 1;
2135: col[2].k = ez;
2136: col[2].loc = LEFT;
2137: col[2].c = 0;
2138: valA[2] = 1.0 / (hy * hy);
2139: col[3].i = ex - 1;
2140: col[3].j = ey;
2141: col[3].k = ez;
2142: col[3].loc = LEFT;
2143: col[3].c = 0;
2144: valA[3] = 1.0 / (hx * hx);
2145: col[4].i = ex + 1;
2146: col[4].j = ey;
2147: col[4].k = ez;
2148: col[4].loc = LEFT;
2149: col[4].c = 0;
2150: valA[4] = 1.0 / (hx * hx);
2151: col[5].i = ex;
2152: col[5].j = ey;
2153: col[5].k = ez - 1;
2154: col[5].loc = LEFT;
2155: col[5].c = 0;
2156: valA[5] = 1.0 / (hz * hz);
2157: col[6].i = ex;
2158: col[6].j = ey;
2159: col[6].k = ez + 1;
2160: col[6].loc = LEFT;
2161: col[6].c = 0;
2162: valA[6] = 1.0 / (hz * hz);
2163: col[7].i = ex - 1;
2164: col[7].j = ey;
2165: col[7].k = ez;
2166: col[7].loc = ELEMENT;
2167: col[7].c = 0;
2168: valA[7] = 1.0 / hx;
2169: col[8].i = ex;
2170: col[8].j = ey;
2171: col[8].k = ez;
2172: col[8].loc = ELEMENT;
2173: col[8].c = 0;
2174: valA[8] = -1.0 / hx;
2175: }
2176: DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed);
2177: check_vals(ex, ey, ez, nEntries, valA, computed);
2178: }
2180: /* Equation on bottom face of this element */
2181: if (ey == 0) {
2182: /* Bottom boundary velocity Dirichlet */
2183: DMStagStencil row;
2184: const PetscScalar valA = 1.0;
2185: row.i = ex;
2186: row.j = ey;
2187: row.k = ez;
2188: row.loc = DOWN;
2189: row.c = 0;
2190: DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed);
2191: check_vals(ex, ey, ez, 1, &valA, computed);
2192: } else {
2193: /* Y-momentum equation, (v_xx + v_yy + v_zz) - p_y = f^y */
2194: DMStagStencil row, col[9];
2195: PetscScalar valA[9];
2196: PetscInt nEntries;
2198: row.i = ex;
2199: row.j = ey;
2200: row.k = ez;
2201: row.loc = DOWN;
2202: row.c = 0;
2203: if (ex == 0) {
2204: if (ez == 0) {
2205: nEntries = 7;
2206: col[0].i = ex;
2207: col[0].j = ey;
2208: col[0].k = ez;
2209: col[0].loc = DOWN;
2210: col[0].c = 0;
2211: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2212: col[1].i = ex;
2213: col[1].j = ey - 1;
2214: col[1].k = ez;
2215: col[1].loc = DOWN;
2216: col[1].c = 0;
2217: valA[1] = 1.0 / (hy * hy);
2218: col[2].i = ex;
2219: col[2].j = ey + 1;
2220: col[2].k = ez;
2221: col[2].loc = DOWN;
2222: col[2].c = 0;
2223: valA[2] = 1.0 / (hy * hy);
2224: /* Left term missing */
2225: col[3].i = ex + 1;
2226: col[3].j = ey;
2227: col[3].k = ez;
2228: col[3].loc = DOWN;
2229: col[3].c = 0;
2230: valA[3] = 1.0 / (hx * hx);
2231: /* Back term missing */
2232: col[4].i = ex;
2233: col[4].j = ey;
2234: col[4].k = ez + 1;
2235: col[4].loc = DOWN;
2236: col[4].c = 0;
2237: valA[4] = 1.0 / (hz * hz);
2238: col[5].i = ex;
2239: col[5].j = ey - 1;
2240: col[5].k = ez;
2241: col[5].loc = ELEMENT;
2242: col[5].c = 0;
2243: valA[5] = 1.0 / hy;
2244: col[6].i = ex;
2245: col[6].j = ey;
2246: col[6].k = ez;
2247: col[6].loc = ELEMENT;
2248: col[6].c = 0;
2249: valA[6] = -1.0 / hy;
2250: } else if (ez == N[2] - 1) {
2251: nEntries = 7;
2252: col[0].i = ex;
2253: col[0].j = ey;
2254: col[0].k = ez;
2255: col[0].loc = DOWN;
2256: col[0].c = 0;
2257: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2258: col[1].i = ex;
2259: col[1].j = ey - 1;
2260: col[1].k = ez;
2261: col[1].loc = DOWN;
2262: col[1].c = 0;
2263: valA[1] = 1.0 / (hy * hy);
2264: col[2].i = ex;
2265: col[2].j = ey + 1;
2266: col[2].k = ez;
2267: col[2].loc = DOWN;
2268: col[2].c = 0;
2269: valA[2] = 1.0 / (hy * hy);
2270: /* Left term missing */
2271: col[3].i = ex + 1;
2272: col[3].j = ey;
2273: col[3].k = ez;
2274: col[3].loc = DOWN;
2275: col[3].c = 0;
2276: valA[3] = 1.0 / (hx * hx);
2277: col[4].i = ex;
2278: col[4].j = ey;
2279: col[4].k = ez - 1;
2280: col[4].loc = DOWN;
2281: col[4].c = 0;
2282: valA[4] = 1.0 / (hz * hz);
2283: /* Front term missing */
2284: col[5].i = ex;
2285: col[5].j = ey - 1;
2286: col[5].k = ez;
2287: col[5].loc = ELEMENT;
2288: col[5].c = 0;
2289: valA[5] = 1.0 / hy;
2290: col[6].i = ex;
2291: col[6].j = ey;
2292: col[6].k = ez;
2293: col[6].loc = ELEMENT;
2294: col[6].c = 0;
2295: valA[6] = -1.0 / hy;
2296: } else {
2297: nEntries = 8;
2298: col[0].i = ex;
2299: col[0].j = ey;
2300: col[0].k = ez;
2301: col[0].loc = DOWN;
2302: col[0].c = 0;
2303: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2304: col[1].i = ex;
2305: col[1].j = ey - 1;
2306: col[1].k = ez;
2307: col[1].loc = DOWN;
2308: col[1].c = 0;
2309: valA[1] = 1.0 / (hy * hy);
2310: col[2].i = ex;
2311: col[2].j = ey + 1;
2312: col[2].k = ez;
2313: col[2].loc = DOWN;
2314: col[2].c = 0;
2315: valA[2] = 1.0 / (hy * hy);
2316: /* Left term missing */
2317: col[3].i = ex + 1;
2318: col[3].j = ey;
2319: col[3].k = ez;
2320: col[3].loc = DOWN;
2321: col[3].c = 0;
2322: valA[3] = 1.0 / (hx * hx);
2323: col[4].i = ex;
2324: col[4].j = ey;
2325: col[4].k = ez - 1;
2326: col[4].loc = DOWN;
2327: col[4].c = 0;
2328: valA[4] = 1.0 / (hz * hz);
2329: col[5].i = ex;
2330: col[5].j = ey;
2331: col[5].k = ez + 1;
2332: col[5].loc = DOWN;
2333: col[5].c = 0;
2334: valA[5] = 1.0 / (hz * hz);
2335: col[6].i = ex;
2336: col[6].j = ey - 1;
2337: col[6].k = ez;
2338: col[6].loc = ELEMENT;
2339: col[6].c = 0;
2340: valA[6] = 1.0 / hy;
2341: col[7].i = ex;
2342: col[7].j = ey;
2343: col[7].k = ez;
2344: col[7].loc = ELEMENT;
2345: col[7].c = 0;
2346: valA[7] = -1.0 / hy;
2347: }
2348: } else if (ex == N[0] - 1) {
2349: if (ez == 0) {
2350: nEntries = 7;
2351: col[0].i = ex;
2352: col[0].j = ey;
2353: col[0].k = ez;
2354: col[0].loc = DOWN;
2355: col[0].c = 0;
2356: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2357: col[1].i = ex;
2358: col[1].j = ey - 1;
2359: col[1].k = ez;
2360: col[1].loc = DOWN;
2361: col[1].c = 0;
2362: valA[1] = 1.0 / (hy * hy);
2363: col[2].i = ex;
2364: col[2].j = ey + 1;
2365: col[2].k = ez;
2366: col[2].loc = DOWN;
2367: col[2].c = 0;
2368: valA[2] = 1.0 / (hy * hy);
2369: col[3].i = ex - 1;
2370: col[3].j = ey;
2371: col[3].k = ez;
2372: col[3].loc = DOWN;
2373: col[3].c = 0;
2374: valA[3] = 1.0 / (hx * hx);
2375: /* Right term missing */
2376: /* Back term missing */
2377: col[4].i = ex;
2378: col[4].j = ey;
2379: col[4].k = ez + 1;
2380: col[4].loc = DOWN;
2381: col[4].c = 0;
2382: valA[4] = 1.0 / (hz * hz);
2383: col[5].i = ex;
2384: col[5].j = ey - 1;
2385: col[5].k = ez;
2386: col[5].loc = ELEMENT;
2387: col[5].c = 0;
2388: valA[5] = 1.0 / hy;
2389: col[6].i = ex;
2390: col[6].j = ey;
2391: col[6].k = ez;
2392: col[6].loc = ELEMENT;
2393: col[6].c = 0;
2394: valA[6] = -1.0 / hy;
2395: } else if (ez == N[2] - 1) {
2396: nEntries = 7;
2397: col[0].i = ex;
2398: col[0].j = ey;
2399: col[0].k = ez;
2400: col[0].loc = DOWN;
2401: col[0].c = 0;
2402: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2403: col[1].i = ex;
2404: col[1].j = ey - 1;
2405: col[1].k = ez;
2406: col[1].loc = DOWN;
2407: col[1].c = 0;
2408: valA[1] = 1.0 / (hy * hy);
2409: col[2].i = ex;
2410: col[2].j = ey + 1;
2411: col[2].k = ez;
2412: col[2].loc = DOWN;
2413: col[2].c = 0;
2414: valA[2] = 1.0 / (hy * hy);
2415: col[3].i = ex - 1;
2416: col[3].j = ey;
2417: col[3].k = ez;
2418: col[3].loc = DOWN;
2419: col[3].c = 0;
2420: valA[3] = 1.0 / (hx * hx);
2421: /* Right term missing */
2422: col[4].i = ex;
2423: col[4].j = ey;
2424: col[4].k = ez - 1;
2425: col[4].loc = DOWN;
2426: col[4].c = 0;
2427: valA[4] = 1.0 / (hz * hz);
2428: /* Front term missing */
2429: col[5].i = ex;
2430: col[5].j = ey - 1;
2431: col[5].k = ez;
2432: col[5].loc = ELEMENT;
2433: col[5].c = 0;
2434: valA[5] = 1.0 / hy;
2435: col[6].i = ex;
2436: col[6].j = ey;
2437: col[6].k = ez;
2438: col[6].loc = ELEMENT;
2439: col[6].c = 0;
2440: valA[6] = -1.0 / hy;
2441: } else {
2442: nEntries = 8;
2443: col[0].i = ex;
2444: col[0].j = ey;
2445: col[0].k = ez;
2446: col[0].loc = DOWN;
2447: col[0].c = 0;
2448: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2449: col[1].i = ex;
2450: col[1].j = ey - 1;
2451: col[1].k = ez;
2452: col[1].loc = DOWN;
2453: col[1].c = 0;
2454: valA[1] = 1.0 / (hy * hy);
2455: col[2].i = ex;
2456: col[2].j = ey + 1;
2457: col[2].k = ez;
2458: col[2].loc = DOWN;
2459: col[2].c = 0;
2460: valA[2] = 1.0 / (hy * hy);
2461: col[3].i = ex - 1;
2462: col[3].j = ey;
2463: col[3].k = ez;
2464: col[3].loc = DOWN;
2465: col[3].c = 0;
2466: valA[3] = 1.0 / (hx * hx);
2467: /* Right term missing */
2468: col[4].i = ex;
2469: col[4].j = ey;
2470: col[4].k = ez - 1;
2471: col[4].loc = DOWN;
2472: col[4].c = 0;
2473: valA[4] = 1.0 / (hz * hz);
2474: col[5].i = ex;
2475: col[5].j = ey;
2476: col[5].k = ez + 1;
2477: col[5].loc = DOWN;
2478: col[5].c = 0;
2479: valA[5] = 1.0 / (hz * hz);
2480: col[6].i = ex;
2481: col[6].j = ey - 1;
2482: col[6].k = ez;
2483: col[6].loc = ELEMENT;
2484: col[6].c = 0;
2485: valA[6] = 1.0 / hy;
2486: col[7].i = ex;
2487: col[7].j = ey;
2488: col[7].k = ez;
2489: col[7].loc = ELEMENT;
2490: col[7].c = 0;
2491: valA[7] = -1.0 / hy;
2492: }
2493: } else if (ez == 0) {
2494: nEntries = 8;
2495: col[0].i = ex;
2496: col[0].j = ey;
2497: col[0].k = ez;
2498: col[0].loc = DOWN;
2499: col[0].c = 0;
2500: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2501: col[1].i = ex;
2502: col[1].j = ey - 1;
2503: col[1].k = ez;
2504: col[1].loc = DOWN;
2505: col[1].c = 0;
2506: valA[1] = 1.0 / (hy * hy);
2507: col[2].i = ex;
2508: col[2].j = ey + 1;
2509: col[2].k = ez;
2510: col[2].loc = DOWN;
2511: col[2].c = 0;
2512: valA[2] = 1.0 / (hy * hy);
2513: col[3].i = ex - 1;
2514: col[3].j = ey;
2515: col[3].k = ez;
2516: col[3].loc = DOWN;
2517: col[3].c = 0;
2518: valA[3] = 1.0 / (hx * hx);
2519: col[4].i = ex + 1;
2520: col[4].j = ey;
2521: col[4].k = ez;
2522: col[4].loc = DOWN;
2523: col[4].c = 0;
2524: valA[4] = 1.0 / (hx * hx);
2525: /* Back term missing */
2526: col[5].i = ex;
2527: col[5].j = ey;
2528: col[5].k = ez + 1;
2529: col[5].loc = DOWN;
2530: col[5].c = 0;
2531: valA[5] = 1.0 / (hz * hz);
2532: col[6].i = ex;
2533: col[6].j = ey - 1;
2534: col[6].k = ez;
2535: col[6].loc = ELEMENT;
2536: col[6].c = 0;
2537: valA[6] = 1.0 / hy;
2538: col[7].i = ex;
2539: col[7].j = ey;
2540: col[7].k = ez;
2541: col[7].loc = ELEMENT;
2542: col[7].c = 0;
2543: valA[7] = -1.0 / hy;
2544: } else if (ez == N[2] - 1) {
2545: nEntries = 8;
2546: col[0].i = ex;
2547: col[0].j = ey;
2548: col[0].k = ez;
2549: col[0].loc = DOWN;
2550: col[0].c = 0;
2551: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2552: col[1].i = ex;
2553: col[1].j = ey - 1;
2554: col[1].k = ez;
2555: col[1].loc = DOWN;
2556: col[1].c = 0;
2557: valA[1] = 1.0 / (hy * hy);
2558: col[2].i = ex;
2559: col[2].j = ey + 1;
2560: col[2].k = ez;
2561: col[2].loc = DOWN;
2562: col[2].c = 0;
2563: valA[2] = 1.0 / (hy * hy);
2564: col[3].i = ex - 1;
2565: col[3].j = ey;
2566: col[3].k = ez;
2567: col[3].loc = DOWN;
2568: col[3].c = 0;
2569: valA[3] = 1.0 / (hx * hx);
2570: col[4].i = ex + 1;
2571: col[4].j = ey;
2572: col[4].k = ez;
2573: col[4].loc = DOWN;
2574: col[4].c = 0;
2575: valA[4] = 1.0 / (hx * hx);
2576: col[5].i = ex;
2577: col[5].j = ey;
2578: col[5].k = ez - 1;
2579: col[5].loc = DOWN;
2580: col[5].c = 0;
2581: valA[5] = 1.0 / (hz * hz);
2582: /* Front term missing */
2583: col[6].i = ex;
2584: col[6].j = ey - 1;
2585: col[6].k = ez;
2586: col[6].loc = ELEMENT;
2587: col[6].c = 0;
2588: valA[6] = 1.0 / hy;
2589: col[7].i = ex;
2590: col[7].j = ey;
2591: col[7].k = ez;
2592: col[7].loc = ELEMENT;
2593: col[7].c = 0;
2594: valA[7] = -1.0 / hy;
2595: } else {
2596: nEntries = 9;
2597: col[0].i = ex;
2598: col[0].j = ey;
2599: col[0].k = ez;
2600: col[0].loc = DOWN;
2601: col[0].c = 0;
2602: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2603: col[1].i = ex;
2604: col[1].j = ey - 1;
2605: col[1].k = ez;
2606: col[1].loc = DOWN;
2607: col[1].c = 0;
2608: valA[1] = 1.0 / (hy * hy);
2609: col[2].i = ex;
2610: col[2].j = ey + 1;
2611: col[2].k = ez;
2612: col[2].loc = DOWN;
2613: col[2].c = 0;
2614: valA[2] = 1.0 / (hy * hy);
2615: col[3].i = ex - 1;
2616: col[3].j = ey;
2617: col[3].k = ez;
2618: col[3].loc = DOWN;
2619: col[3].c = 0;
2620: valA[3] = 1.0 / (hx * hx);
2621: col[4].i = ex + 1;
2622: col[4].j = ey;
2623: col[4].k = ez;
2624: col[4].loc = DOWN;
2625: col[4].c = 0;
2626: valA[4] = 1.0 / (hx * hx);
2627: col[5].i = ex;
2628: col[5].j = ey;
2629: col[5].k = ez - 1;
2630: col[5].loc = DOWN;
2631: col[5].c = 0;
2632: valA[5] = 1.0 / (hz * hz);
2633: col[6].i = ex;
2634: col[6].j = ey;
2635: col[6].k = ez + 1;
2636: col[6].loc = DOWN;
2637: col[6].c = 0;
2638: valA[6] = 1.0 / (hz * hz);
2639: col[7].i = ex;
2640: col[7].j = ey - 1;
2641: col[7].k = ez;
2642: col[7].loc = ELEMENT;
2643: col[7].c = 0;
2644: valA[7] = 1.0 / hy;
2645: col[8].i = ex;
2646: col[8].j = ey;
2647: col[8].k = ez;
2648: col[8].loc = ELEMENT;
2649: col[8].c = 0;
2650: valA[8] = -1.0 / hy;
2651: }
2652: DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed);
2653: check_vals(ex, ey, ez, nEntries, valA, computed);
2654: }
2656: /* Equation on back face of this element */
2657: if (ez == 0) {
2658: /* Back boundary velocity Dirichlet */
2659: DMStagStencil row;
2660: const PetscScalar valA = 1.0;
2661: row.i = ex;
2662: row.j = ey;
2663: row.k = ez;
2664: row.loc = BACK;
2665: row.c = 0;
2666: DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed);
2667: check_vals(ex, ey, ez, 1, &valA, computed);
2668: } else {
2669: /* Z-momentum equation, (w_xx + w_yy + w_zz) - p_z = f^z */
2670: DMStagStencil row, col[9];
2671: PetscScalar valA[9];
2672: PetscInt nEntries;
2674: row.i = ex;
2675: row.j = ey;
2676: row.k = ez;
2677: row.loc = BACK;
2678: row.c = 0;
2679: if (ex == 0) {
2680: if (ey == 0) {
2681: nEntries = 7;
2682: col[0].i = ex;
2683: col[0].j = ey;
2684: col[0].k = ez;
2685: col[0].loc = BACK;
2686: col[0].c = 0;
2687: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2688: /* Down term missing */
2689: col[1].i = ex;
2690: col[1].j = ey + 1;
2691: col[1].k = ez;
2692: col[1].loc = BACK;
2693: col[1].c = 0;
2694: valA[1] = 1.0 / (hy * hy);
2695: /* Left term missing */
2696: col[2].i = ex + 1;
2697: col[2].j = ey;
2698: col[2].k = ez;
2699: col[2].loc = BACK;
2700: col[2].c = 0;
2701: valA[2] = 1.0 / (hx * hx);
2702: col[3].i = ex;
2703: col[3].j = ey;
2704: col[3].k = ez - 1;
2705: col[3].loc = BACK;
2706: col[3].c = 0;
2707: valA[3] = 1.0 / (hz * hz);
2708: col[4].i = ex;
2709: col[4].j = ey;
2710: col[4].k = ez + 1;
2711: col[4].loc = BACK;
2712: col[4].c = 0;
2713: valA[4] = 1.0 / (hz * hz);
2714: col[5].i = ex;
2715: col[5].j = ey;
2716: col[5].k = ez - 1;
2717: col[5].loc = ELEMENT;
2718: col[5].c = 0;
2719: valA[5] = 1.0 / hz;
2720: col[6].i = ex;
2721: col[6].j = ey;
2722: col[6].k = ez;
2723: col[6].loc = ELEMENT;
2724: col[6].c = 0;
2725: valA[6] = -1.0 / hz;
2726: } else if (ey == N[1] - 1) {
2727: nEntries = 7;
2728: col[0].i = ex;
2729: col[0].j = ey;
2730: col[0].k = ez;
2731: col[0].loc = BACK;
2732: col[0].c = 0;
2733: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2734: col[1].i = ex;
2735: col[1].j = ey - 1;
2736: col[1].k = ez;
2737: col[1].loc = BACK;
2738: col[1].c = 0;
2739: valA[1] = 1.0 / (hy * hy);
2740: /* Up term missing */
2741: /* Left term missing */
2742: col[2].i = ex + 1;
2743: col[2].j = ey;
2744: col[2].k = ez;
2745: col[2].loc = BACK;
2746: col[2].c = 0;
2747: valA[2] = 1.0 / (hx * hx);
2748: col[3].i = ex;
2749: col[3].j = ey;
2750: col[3].k = ez - 1;
2751: col[3].loc = BACK;
2752: col[3].c = 0;
2753: valA[3] = 1.0 / (hz * hz);
2754: col[4].i = ex;
2755: col[4].j = ey;
2756: col[4].k = ez + 1;
2757: col[4].loc = BACK;
2758: col[4].c = 0;
2759: valA[4] = 1.0 / (hz * hz);
2760: col[5].i = ex;
2761: col[5].j = ey;
2762: col[5].k = ez - 1;
2763: col[5].loc = ELEMENT;
2764: col[5].c = 0;
2765: valA[5] = 1.0 / hz;
2766: col[6].i = ex;
2767: col[6].j = ey;
2768: col[6].k = ez;
2769: col[6].loc = ELEMENT;
2770: col[6].c = 0;
2771: valA[6] = -1.0 / hz;
2772: } else {
2773: nEntries = 8;
2774: col[0].i = ex;
2775: col[0].j = ey;
2776: col[0].k = ez;
2777: col[0].loc = BACK;
2778: col[0].c = 0;
2779: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2780: col[1].i = ex;
2781: col[1].j = ey - 1;
2782: col[1].k = ez;
2783: col[1].loc = BACK;
2784: col[1].c = 0;
2785: valA[1] = 1.0 / (hy * hy);
2786: col[2].i = ex;
2787: col[2].j = ey + 1;
2788: col[2].k = ez;
2789: col[2].loc = BACK;
2790: col[2].c = 0;
2791: valA[2] = 1.0 / (hy * hy);
2792: /* Left term missing */
2793: col[3].i = ex + 1;
2794: col[3].j = ey;
2795: col[3].k = ez;
2796: col[3].loc = BACK;
2797: col[3].c = 0;
2798: valA[3] = 1.0 / (hx * hx);
2799: col[4].i = ex;
2800: col[4].j = ey;
2801: col[4].k = ez - 1;
2802: col[4].loc = BACK;
2803: col[4].c = 0;
2804: valA[4] = 1.0 / (hz * hz);
2805: col[5].i = ex;
2806: col[5].j = ey;
2807: col[5].k = ez + 1;
2808: col[5].loc = BACK;
2809: col[5].c = 0;
2810: valA[5] = 1.0 / (hz * hz);
2811: col[6].i = ex;
2812: col[6].j = ey;
2813: col[6].k = ez - 1;
2814: col[6].loc = ELEMENT;
2815: col[6].c = 0;
2816: valA[6] = 1.0 / hz;
2817: col[7].i = ex;
2818: col[7].j = ey;
2819: col[7].k = ez;
2820: col[7].loc = ELEMENT;
2821: col[7].c = 0;
2822: valA[7] = -1.0 / hz;
2823: }
2824: } else if (ex == N[0] - 1) {
2825: if (ey == 0) {
2826: nEntries = 7;
2827: col[0].i = ex;
2828: col[0].j = ey;
2829: col[0].k = ez;
2830: col[0].loc = BACK;
2831: col[0].c = 0;
2832: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2833: /* Down term missing */
2834: col[1].i = ex;
2835: col[1].j = ey + 1;
2836: col[1].k = ez;
2837: col[1].loc = BACK;
2838: col[1].c = 0;
2839: valA[1] = 1.0 / (hy * hy);
2840: col[2].i = ex - 1;
2841: col[2].j = ey;
2842: col[2].k = ez;
2843: col[2].loc = BACK;
2844: col[2].c = 0;
2845: valA[2] = 1.0 / (hx * hx);
2846: /* Right term missing */
2847: col[3].i = ex;
2848: col[3].j = ey;
2849: col[3].k = ez - 1;
2850: col[3].loc = BACK;
2851: col[3].c = 0;
2852: valA[3] = 1.0 / (hz * hz);
2853: col[4].i = ex;
2854: col[4].j = ey;
2855: col[4].k = ez + 1;
2856: col[4].loc = BACK;
2857: col[4].c = 0;
2858: valA[4] = 1.0 / (hz * hz);
2859: col[5].i = ex;
2860: col[5].j = ey;
2861: col[5].k = ez - 1;
2862: col[5].loc = ELEMENT;
2863: col[5].c = 0;
2864: valA[5] = 1.0 / hz;
2865: col[6].i = ex;
2866: col[6].j = ey;
2867: col[6].k = ez;
2868: col[6].loc = ELEMENT;
2869: col[6].c = 0;
2870: valA[6] = -1.0 / hz;
2871: } else if (ey == N[1] - 1) {
2872: nEntries = 7;
2873: col[0].i = ex;
2874: col[0].j = ey;
2875: col[0].k = ez;
2876: col[0].loc = BACK;
2877: col[0].c = 0;
2878: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2879: col[1].i = ex;
2880: col[1].j = ey - 1;
2881: col[1].k = ez;
2882: col[1].loc = BACK;
2883: col[1].c = 0;
2884: valA[1] = 1.0 / (hy * hy);
2885: /* Up term missing */
2886: col[2].i = ex - 1;
2887: col[2].j = ey;
2888: col[2].k = ez;
2889: col[2].loc = BACK;
2890: col[2].c = 0;
2891: valA[2] = 1.0 / (hx * hx);
2892: /* Right term missing */
2893: col[3].i = ex;
2894: col[3].j = ey;
2895: col[3].k = ez - 1;
2896: col[3].loc = BACK;
2897: col[3].c = 0;
2898: valA[3] = 1.0 / (hz * hz);
2899: col[4].i = ex;
2900: col[4].j = ey;
2901: col[4].k = ez + 1;
2902: col[4].loc = BACK;
2903: col[4].c = 0;
2904: valA[4] = 1.0 / (hz * hz);
2905: col[5].i = ex;
2906: col[5].j = ey;
2907: col[5].k = ez - 1;
2908: col[5].loc = ELEMENT;
2909: col[5].c = 0;
2910: valA[5] = 1.0 / hz;
2911: col[6].i = ex;
2912: col[6].j = ey;
2913: col[6].k = ez;
2914: col[6].loc = ELEMENT;
2915: col[6].c = 0;
2916: valA[6] = -1.0 / hz;
2917: } else {
2918: nEntries = 8;
2919: col[0].i = ex;
2920: col[0].j = ey;
2921: col[0].k = ez;
2922: col[0].loc = BACK;
2923: col[0].c = 0;
2924: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2925: col[1].i = ex;
2926: col[1].j = ey - 1;
2927: col[1].k = ez;
2928: col[1].loc = BACK;
2929: col[1].c = 0;
2930: valA[1] = 1.0 / (hy * hy);
2931: col[2].i = ex;
2932: col[2].j = ey + 1;
2933: col[2].k = ez;
2934: col[2].loc = BACK;
2935: col[2].c = 0;
2936: valA[2] = 1.0 / (hy * hy);
2937: col[3].i = ex - 1;
2938: col[3].j = ey;
2939: col[3].k = ez;
2940: col[3].loc = BACK;
2941: col[3].c = 0;
2942: valA[3] = 1.0 / (hx * hx);
2943: /* Right term missing */
2944: col[4].i = ex;
2945: col[4].j = ey;
2946: col[4].k = ez - 1;
2947: col[4].loc = BACK;
2948: col[4].c = 0;
2949: valA[4] = 1.0 / (hz * hz);
2950: col[5].i = ex;
2951: col[5].j = ey;
2952: col[5].k = ez + 1;
2953: col[5].loc = BACK;
2954: col[5].c = 0;
2955: valA[5] = 1.0 / (hz * hz);
2956: col[6].i = ex;
2957: col[6].j = ey;
2958: col[6].k = ez - 1;
2959: col[6].loc = ELEMENT;
2960: col[6].c = 0;
2961: valA[6] = 1.0 / hz;
2962: col[7].i = ex;
2963: col[7].j = ey;
2964: col[7].k = ez;
2965: col[7].loc = ELEMENT;
2966: col[7].c = 0;
2967: valA[7] = -1.0 / hz;
2968: }
2969: } else if (ey == 0) {
2970: nEntries = 8;
2971: col[0].i = ex;
2972: col[0].j = ey;
2973: col[0].k = ez;
2974: col[0].loc = BACK;
2975: col[0].c = 0;
2976: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2977: /* Down term missing */
2978: col[1].i = ex;
2979: col[1].j = ey + 1;
2980: col[1].k = ez;
2981: col[1].loc = BACK;
2982: col[1].c = 0;
2983: valA[1] = 1.0 / (hy * hy);
2984: col[2].i = ex - 1;
2985: col[2].j = ey;
2986: col[2].k = ez;
2987: col[2].loc = BACK;
2988: col[2].c = 0;
2989: valA[2] = 1.0 / (hx * hx);
2990: col[3].i = ex + 1;
2991: col[3].j = ey;
2992: col[3].k = ez;
2993: col[3].loc = BACK;
2994: col[3].c = 0;
2995: valA[3] = 1.0 / (hx * hx);
2996: col[4].i = ex;
2997: col[4].j = ey;
2998: col[4].k = ez - 1;
2999: col[4].loc = BACK;
3000: col[4].c = 0;
3001: valA[4] = 1.0 / (hz * hz);
3002: col[5].i = ex;
3003: col[5].j = ey;
3004: col[5].k = ez + 1;
3005: col[5].loc = BACK;
3006: col[5].c = 0;
3007: valA[5] = 1.0 / (hz * hz);
3008: col[6].i = ex;
3009: col[6].j = ey;
3010: col[6].k = ez - 1;
3011: col[6].loc = ELEMENT;
3012: col[6].c = 0;
3013: valA[6] = 1.0 / hz;
3014: col[7].i = ex;
3015: col[7].j = ey;
3016: col[7].k = ez;
3017: col[7].loc = ELEMENT;
3018: col[7].c = 0;
3019: valA[7] = -1.0 / hz;
3020: } else if (ey == N[1] - 1) {
3021: nEntries = 8;
3022: col[0].i = ex;
3023: col[0].j = ey;
3024: col[0].k = ez;
3025: col[0].loc = BACK;
3026: col[0].c = 0;
3027: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
3028: col[1].i = ex;
3029: col[1].j = ey - 1;
3030: col[1].k = ez;
3031: col[1].loc = BACK;
3032: col[1].c = 0;
3033: valA[1] = 1.0 / (hy * hy);
3034: /* Up term missing */
3035: col[2].i = ex - 1;
3036: col[2].j = ey;
3037: col[2].k = ez;
3038: col[2].loc = BACK;
3039: col[2].c = 0;
3040: valA[2] = 1.0 / (hx * hx);
3041: col[3].i = ex + 1;
3042: col[3].j = ey;
3043: col[3].k = ez;
3044: col[3].loc = BACK;
3045: col[3].c = 0;
3046: valA[3] = 1.0 / (hx * hx);
3047: col[4].i = ex;
3048: col[4].j = ey;
3049: col[4].k = ez - 1;
3050: col[4].loc = BACK;
3051: col[4].c = 0;
3052: valA[4] = 1.0 / (hz * hz);
3053: col[5].i = ex;
3054: col[5].j = ey;
3055: col[5].k = ez + 1;
3056: col[5].loc = BACK;
3057: col[5].c = 0;
3058: valA[5] = 1.0 / (hz * hz);
3059: col[6].i = ex;
3060: col[6].j = ey;
3061: col[6].k = ez - 1;
3062: col[6].loc = ELEMENT;
3063: col[6].c = 0;
3064: valA[6] = 1.0 / hz;
3065: col[7].i = ex;
3066: col[7].j = ey;
3067: col[7].k = ez;
3068: col[7].loc = ELEMENT;
3069: col[7].c = 0;
3070: valA[7] = -1.0 / hz;
3071: } else {
3072: nEntries = 9;
3073: col[0].i = ex;
3074: col[0].j = ey;
3075: col[0].k = ez;
3076: col[0].loc = BACK;
3077: col[0].c = 0;
3078: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
3079: col[1].i = ex;
3080: col[1].j = ey - 1;
3081: col[1].k = ez;
3082: col[1].loc = BACK;
3083: col[1].c = 0;
3084: valA[1] = 1.0 / (hy * hy);
3085: col[2].i = ex;
3086: col[2].j = ey + 1;
3087: col[2].k = ez;
3088: col[2].loc = BACK;
3089: col[2].c = 0;
3090: valA[2] = 1.0 / (hy * hy);
3091: col[3].i = ex - 1;
3092: col[3].j = ey;
3093: col[3].k = ez;
3094: col[3].loc = BACK;
3095: col[3].c = 0;
3096: valA[3] = 1.0 / (hx * hx);
3097: col[4].i = ex + 1;
3098: col[4].j = ey;
3099: col[4].k = ez;
3100: col[4].loc = BACK;
3101: col[4].c = 0;
3102: valA[4] = 1.0 / (hx * hx);
3103: col[5].i = ex;
3104: col[5].j = ey;
3105: col[5].k = ez - 1;
3106: col[5].loc = BACK;
3107: col[5].c = 0;
3108: valA[5] = 1.0 / (hz * hz);
3109: col[6].i = ex;
3110: col[6].j = ey;
3111: col[6].k = ez + 1;
3112: col[6].loc = BACK;
3113: col[6].c = 0;
3114: valA[6] = 1.0 / (hz * hz);
3115: col[7].i = ex;
3116: col[7].j = ey;
3117: col[7].k = ez - 1;
3118: col[7].loc = ELEMENT;
3119: col[7].c = 0;
3120: valA[7] = 1.0 / hz;
3121: col[8].i = ex;
3122: col[8].j = ey;
3123: col[8].k = ez;
3124: col[8].loc = ELEMENT;
3125: col[8].c = 0;
3126: valA[8] = -1.0 / hz;
3127: }
3128: DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed);
3129: check_vals(ex, ey, ez, nEntries, valA, computed);
3130: }
3132: /* P equation : u_x + v_y + w_z = g
3133: Note that this includes an explicit zero on the diagonal. This is only needed for
3134: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
3135: {
3136: DMStagStencil row, col[7];
3137: PetscScalar valA[7];
3139: row.i = ex;
3140: row.j = ey;
3141: row.k = ez;
3142: row.loc = ELEMENT;
3143: row.c = 0;
3144: col[0].i = ex;
3145: col[0].j = ey;
3146: col[0].k = ez;
3147: col[0].loc = LEFT;
3148: col[0].c = 0;
3149: valA[0] = -1.0 / hx;
3150: col[1].i = ex;
3151: col[1].j = ey;
3152: col[1].k = ez;
3153: col[1].loc = RIGHT;
3154: col[1].c = 0;
3155: valA[1] = 1.0 / hx;
3156: col[2].i = ex;
3157: col[2].j = ey;
3158: col[2].k = ez;
3159: col[2].loc = DOWN;
3160: col[2].c = 0;
3161: valA[2] = -1.0 / hy;
3162: col[3].i = ex;
3163: col[3].j = ey;
3164: col[3].k = ez;
3165: col[3].loc = UP;
3166: col[3].c = 0;
3167: valA[3] = 1.0 / hy;
3168: col[4].i = ex;
3169: col[4].j = ey;
3170: col[4].k = ez;
3171: col[4].loc = BACK;
3172: col[4].c = 0;
3173: valA[4] = -1.0 / hz;
3174: col[5].i = ex;
3175: col[5].j = ey;
3176: col[5].k = ez;
3177: col[5].loc = FRONT;
3178: col[5].c = 0;
3179: valA[5] = 1.0 / hz;
3180: col[6] = row;
3181: valA[6] = 0.0;
3182: DMStagMatGetValuesStencil(dmSol, A, 1, &row, 7, col, computed);
3183: check_vals(ex, ey, ez, 7, valA, computed);
3184: }
3185: }
3186: }
3187: }
3188: DMStagVecRestoreArrayRead(dmCoord, coordLocal, &arrCoord);
3189: return 0;
3190: }
3192: /*TEST
3194: test:
3195: suffix: 1
3196: nsize: 1
3198: test:
3199: suffix: 2
3200: nsize: 4
3202: TEST*/