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;
44: PetscFunctionBeginUser;
45: PetscCall(PetscInitialize(&argc, &argv, NULL, 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: PetscCall(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: PetscCall(DMSetFromOptions(dmSol));
51: PetscCall(DMSetUp(dmSol));
52: PetscCall(DMStagSetUniformCoordinatesExplicit(dmSol, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
53: }
54: PetscCall(CreateMat(dmSol, &A));
55: PetscCall(CheckMat(dmSol, A));
56: PetscCall(MatDestroy(&A));
57: PetscCall(DMDestroy(&dmSol));
58: PetscCall(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;
73: PetscFunctionBeginUser;
74: PetscCall(DMCreateMatrix(dmSol, pA));
75: A = *pA;
76: PetscCall(DMStagGetCorners(dmSol, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
77: PetscCall(DMStagGetGlobalSizes(dmSol, &N[0], &N[1], &N[2]));
78: PetscCheck(N[0] >= 2 && N[1] >= 2 && N[2] >= 2, PetscObjectComm((PetscObject)dmSol), PETSC_ERR_ARG_SIZ, "This example requires at least two elements in each dimensions");
79: PetscCall(DMStagGetIsLastRank(dmSol, &isLastRankx, &isLastRanky, &isLastRankz));
80: PetscCall(DMStagGetIsFirstRank(dmSol, &isFirstRankx, &isFirstRanky, &isFirstRankz));
81: hx = 1.0 / N[0];
82: hy = 1.0 / N[1];
83: hz = 1.0 / N[2];
84: PetscCall(DMGetCoordinateDM(dmSol, &dmCoord));
85: PetscCall(DMGetCoordinatesLocal(dmSol, &coordLocal));
86: PetscCall(DMStagVecGetArrayRead(dmCoord, coordLocal, &arrCoord));
87: for (d = 0; d < 3; ++d) {
88: PetscCall(DMStagGetLocationSlot(dmCoord, ELEMENT, d, &icp[d]));
89: PetscCall(DMStagGetLocationSlot(dmCoord, LEFT, d, &icux[d]));
90: PetscCall(DMStagGetLocationSlot(dmCoord, DOWN, d, &icuy[d]));
91: PetscCall(DMStagGetLocationSlot(dmCoord, BACK, d, &icuz[d]));
92: PetscCall(DMStagGetLocationSlot(dmCoord, RIGHT, d, &icux_right[d]));
93: PetscCall(DMStagGetLocationSlot(dmCoord, UP, d, &icuy_up[d]));
94: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 7, col, valA, INSERT_VALUES));
1607: }
1608: }
1609: }
1610: }
1611: PetscCall(DMStagVecRestoreArrayRead(dmCoord, coordLocal, &arrCoord));
1612: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1613: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1614: PetscFunctionReturn(PETSC_SUCCESS);
1615: }
1617: /* A helper function to check values */
1618: static PetscErrorCode check_vals(PetscInt ex, PetscInt ey, PetscInt ez, PetscInt n, const PetscScalar *ref, const PetscScalar *computed)
1619: {
1620: PetscFunctionBeginUser;
1621: for (PetscInt i = 0; i < n; ++i) {
1622: PetscCheck(ref[i] == computed[i], PETSC_COMM_SELF, PETSC_ERR_PLIB, "(%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ") Assertion Failure. (ref[%" PetscInt_FMT "]) %g != %g (computed)[%" PetscInt_FMT "]", ex, ey, ez, i, (double)PetscRealPart(ref[i]), (double)PetscRealPart(computed[i]), i);
1623: }
1624: PetscFunctionReturn(PETSC_SUCCESS);
1625: }
1627: /* The same function as above, but getting and checking values, instead of setting them */
1628: static PetscErrorCode CheckMat(DM dmSol, Mat A)
1629: {
1630: Vec coordLocal;
1631: PetscInt startx, starty, startz, N[3], nx, ny, nz, ex, ey, ez, d;
1632: PetscInt icp[3], icux[3], icuy[3], icuz[3], icux_right[3], icuy_up[3], icuz_front[3];
1633: PetscBool isLastRankx, isLastRanky, isLastRankz, isFirstRankx, isFirstRanky, isFirstRankz;
1634: PetscReal hx, hy, hz;
1635: DM dmCoord;
1636: PetscScalar ****arrCoord;
1637: PetscScalar computed[1024];
1639: PetscFunctionBeginUser;
1640: PetscCall(DMStagGetCorners(dmSol, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
1641: PetscCall(DMStagGetGlobalSizes(dmSol, &N[0], &N[1], &N[2]));
1642: PetscCheck(N[0] >= 2 && N[1] >= 2 && N[2] >= 2, PetscObjectComm((PetscObject)dmSol), PETSC_ERR_ARG_SIZ, "This example requires at least two elements in each dimensions");
1643: PetscCall(DMStagGetIsLastRank(dmSol, &isLastRankx, &isLastRanky, &isLastRankz));
1644: PetscCall(DMStagGetIsFirstRank(dmSol, &isFirstRankx, &isFirstRanky, &isFirstRankz));
1645: hx = 1.0 / N[0];
1646: hy = 1.0 / N[1];
1647: hz = 1.0 / N[2];
1648: PetscCall(DMGetCoordinateDM(dmSol, &dmCoord));
1649: PetscCall(DMGetCoordinatesLocal(dmSol, &coordLocal));
1650: PetscCall(DMStagVecGetArrayRead(dmCoord, coordLocal, &arrCoord));
1651: for (d = 0; d < 3; ++d) {
1652: PetscCall(DMStagGetLocationSlot(dmCoord, ELEMENT, d, &icp[d]));
1653: PetscCall(DMStagGetLocationSlot(dmCoord, LEFT, d, &icux[d]));
1654: PetscCall(DMStagGetLocationSlot(dmCoord, DOWN, d, &icuy[d]));
1655: PetscCall(DMStagGetLocationSlot(dmCoord, BACK, d, &icuz[d]));
1656: PetscCall(DMStagGetLocationSlot(dmCoord, RIGHT, d, &icux_right[d]));
1657: PetscCall(DMStagGetLocationSlot(dmCoord, UP, d, &icuy_up[d]));
1658: PetscCall(DMStagGetLocationSlot(dmCoord, FRONT, d, &icuz_front[d]));
1659: }
1661: for (ez = startz; ez < startz + nz; ++ez) {
1662: for (ey = starty; ey < starty + ny; ++ey) {
1663: for (ex = startx; ex < startx + nx; ++ex) {
1664: if (ex == N[0] - 1) {
1665: /* Right Boundary velocity Dirichlet */
1666: DMStagStencil row;
1667: const PetscScalar valA = 1.0;
1668: row.i = ex;
1669: row.j = ey;
1670: row.k = ez;
1671: row.loc = RIGHT;
1672: row.c = 0;
1673: PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1674: PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
1675: }
1676: if (ey == N[1] - 1) {
1677: /* Top boundary velocity Dirichlet */
1678: DMStagStencil row;
1679: const PetscScalar valA = 1.0;
1680: row.i = ex;
1681: row.j = ey;
1682: row.k = ez;
1683: row.loc = UP;
1684: row.c = 0;
1685: PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1686: PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
1687: }
1688: if (ez == N[2] - 1) {
1689: /* Top boundary velocity Dirichlet */
1690: DMStagStencil row;
1691: const PetscScalar valA = 1.0;
1692: row.i = ex;
1693: row.j = ey;
1694: row.k = ez;
1695: row.loc = FRONT;
1696: row.c = 0;
1697: PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1698: PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
1699: }
1701: /* Equation on left face of this element */
1702: if (ex == 0) {
1703: /* Left velocity Dirichlet */
1704: DMStagStencil row;
1705: const PetscScalar valA = 1.0;
1706: row.i = ex;
1707: row.j = ey;
1708: row.k = ez;
1709: row.loc = LEFT;
1710: row.c = 0;
1711: PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1712: PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
1713: } else {
1714: /* X-momentum interior equation : (u_xx + u_yy + u_zz) - p_x = f^x */
1715: DMStagStencil row, col[9];
1716: PetscScalar valA[9];
1717: PetscInt nEntries;
1719: row.i = ex;
1720: row.j = ey;
1721: row.k = ez;
1722: row.loc = LEFT;
1723: row.c = 0;
1724: if (ey == 0) {
1725: if (ez == 0) {
1726: nEntries = 7;
1727: col[0].i = ex;
1728: col[0].j = ey;
1729: col[0].k = ez;
1730: col[0].loc = LEFT;
1731: col[0].c = 0;
1732: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
1733: /* Missing down term */
1734: col[1].i = ex;
1735: col[1].j = ey + 1;
1736: col[1].k = ez;
1737: col[1].loc = LEFT;
1738: col[1].c = 0;
1739: valA[1] = 1.0 / (hy * hy);
1740: col[2].i = ex - 1;
1741: col[2].j = ey;
1742: col[2].k = ez;
1743: col[2].loc = LEFT;
1744: col[2].c = 0;
1745: valA[2] = 1.0 / (hx * hx);
1746: col[3].i = ex + 1;
1747: col[3].j = ey;
1748: col[3].k = ez;
1749: col[3].loc = LEFT;
1750: col[3].c = 0;
1751: valA[3] = 1.0 / (hx * hx);
1752: /* Missing back term */
1753: col[4].i = ex;
1754: col[4].j = ey;
1755: col[4].k = ez + 1;
1756: col[4].loc = LEFT;
1757: col[4].c = 0;
1758: valA[4] = 1.0 / (hz * hz);
1759: col[5].i = ex - 1;
1760: col[5].j = ey;
1761: col[5].k = ez;
1762: col[5].loc = ELEMENT;
1763: col[5].c = 0;
1764: valA[5] = 1.0 / hx;
1765: col[6].i = ex;
1766: col[6].j = ey;
1767: col[6].k = ez;
1768: col[6].loc = ELEMENT;
1769: col[6].c = 0;
1770: valA[6] = -1.0 / hx;
1771: } else if (ez == N[2] - 1) {
1772: nEntries = 7;
1773: col[0].i = ex;
1774: col[0].j = ey;
1775: col[0].k = ez;
1776: col[0].loc = LEFT;
1777: col[0].c = 0;
1778: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
1779: /* Missing down term */
1780: col[1].i = ex;
1781: col[1].j = ey + 1;
1782: col[1].k = ez;
1783: col[1].loc = LEFT;
1784: col[1].c = 0;
1785: valA[1] = 1.0 / (hy * hy);
1786: col[2].i = ex - 1;
1787: col[2].j = ey;
1788: col[2].k = ez;
1789: col[2].loc = LEFT;
1790: col[2].c = 0;
1791: valA[2] = 1.0 / (hx * hx);
1792: col[3].i = ex + 1;
1793: col[3].j = ey;
1794: col[3].k = ez;
1795: col[3].loc = LEFT;
1796: col[3].c = 0;
1797: valA[3] = 1.0 / (hx * hx);
1798: col[4].i = ex;
1799: col[4].j = ey;
1800: col[4].k = ez - 1;
1801: col[4].loc = LEFT;
1802: col[4].c = 0;
1803: valA[4] = 1.0 / (hz * hz);
1804: /* Missing front term */
1805: col[5].i = ex - 1;
1806: col[5].j = ey;
1807: col[5].k = ez;
1808: col[5].loc = ELEMENT;
1809: col[5].c = 0;
1810: valA[5] = 1.0 / hx;
1811: col[6].i = ex;
1812: col[6].j = ey;
1813: col[6].k = ez;
1814: col[6].loc = ELEMENT;
1815: col[6].c = 0;
1816: valA[6] = -1.0 / hx;
1817: } else {
1818: nEntries = 8;
1819: col[0].i = ex;
1820: col[0].j = ey;
1821: col[0].k = ez;
1822: col[0].loc = LEFT;
1823: col[0].c = 0;
1824: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1825: /* Missing down term */
1826: col[1].i = ex;
1827: col[1].j = ey + 1;
1828: col[1].k = ez;
1829: col[1].loc = LEFT;
1830: col[1].c = 0;
1831: valA[1] = 1.0 / (hy * hy);
1832: col[2].i = ex - 1;
1833: col[2].j = ey;
1834: col[2].k = ez;
1835: col[2].loc = LEFT;
1836: col[2].c = 0;
1837: valA[2] = 1.0 / (hx * hx);
1838: col[3].i = ex + 1;
1839: col[3].j = ey;
1840: col[3].k = ez;
1841: col[3].loc = LEFT;
1842: col[3].c = 0;
1843: valA[3] = 1.0 / (hx * hx);
1844: col[4].i = ex;
1845: col[4].j = ey;
1846: col[4].k = ez - 1;
1847: col[4].loc = LEFT;
1848: col[4].c = 0;
1849: valA[4] = 1.0 / (hz * hz);
1850: col[5].i = ex;
1851: col[5].j = ey;
1852: col[5].k = ez + 1;
1853: col[5].loc = LEFT;
1854: col[5].c = 0;
1855: valA[5] = 1.0 / (hz * hz);
1856: col[6].i = ex - 1;
1857: col[6].j = ey;
1858: col[6].k = ez;
1859: col[6].loc = ELEMENT;
1860: col[6].c = 0;
1861: valA[6] = 1.0 / hx;
1862: col[7].i = ex;
1863: col[7].j = ey;
1864: col[7].k = ez;
1865: col[7].loc = ELEMENT;
1866: col[7].c = 0;
1867: valA[7] = -1.0 / hx;
1868: }
1869: } else if (ey == N[1] - 1) {
1870: if (ez == 0) {
1871: nEntries = 7;
1872: col[0].i = ex;
1873: col[0].j = ey;
1874: col[0].k = ez;
1875: col[0].loc = LEFT;
1876: col[0].c = 0;
1877: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
1878: col[1].i = ex;
1879: col[1].j = ey - 1;
1880: col[1].k = ez;
1881: col[1].loc = LEFT;
1882: col[1].c = 0;
1883: valA[1] = 1.0 / (hy * hy);
1884: /* Missing up term */
1885: col[2].i = ex - 1;
1886: col[2].j = ey;
1887: col[2].k = ez;
1888: col[2].loc = LEFT;
1889: col[2].c = 0;
1890: valA[2] = 1.0 / (hx * hx);
1891: col[3].i = ex + 1;
1892: col[3].j = ey;
1893: col[3].k = ez;
1894: col[3].loc = LEFT;
1895: col[3].c = 0;
1896: valA[3] = 1.0 / (hx * hx);
1897: /* Missing back entry */
1898: col[4].i = ex;
1899: col[4].j = ey;
1900: col[4].k = ez + 1;
1901: col[4].loc = LEFT;
1902: col[4].c = 0;
1903: valA[4] = 1.0 / (hz * hz);
1904: col[5].i = ex - 1;
1905: col[5].j = ey;
1906: col[5].k = ez;
1907: col[5].loc = ELEMENT;
1908: col[5].c = 0;
1909: valA[5] = 1.0 / hx;
1910: col[6].i = ex;
1911: col[6].j = ey;
1912: col[6].k = ez;
1913: col[6].loc = ELEMENT;
1914: col[6].c = 0;
1915: valA[6] = -1.0 / hx;
1916: } else if (ez == N[2] - 1) {
1917: nEntries = 7;
1918: col[0].i = ex;
1919: col[0].j = ey;
1920: col[0].k = ez;
1921: col[0].loc = LEFT;
1922: col[0].c = 0;
1923: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
1924: col[1].i = ex;
1925: col[1].j = ey - 1;
1926: col[1].k = ez;
1927: col[1].loc = LEFT;
1928: col[1].c = 0;
1929: valA[1] = 1.0 / (hy * hy);
1930: /* Missing up term */
1931: col[2].i = ex - 1;
1932: col[2].j = ey;
1933: col[2].k = ez;
1934: col[2].loc = LEFT;
1935: col[2].c = 0;
1936: valA[2] = 1.0 / (hx * hx);
1937: col[3].i = ex + 1;
1938: col[3].j = ey;
1939: col[3].k = ez;
1940: col[3].loc = LEFT;
1941: col[3].c = 0;
1942: valA[3] = 1.0 / (hx * hx);
1943: col[4].i = ex;
1944: col[4].j = ey;
1945: col[4].k = ez - 1;
1946: col[4].loc = LEFT;
1947: col[4].c = 0;
1948: valA[4] = 1.0 / (hz * hz);
1949: /* Missing front term */
1950: col[5].i = ex - 1;
1951: col[5].j = ey;
1952: col[5].k = ez;
1953: col[5].loc = ELEMENT;
1954: col[5].c = 0;
1955: valA[5] = 1.0 / hx;
1956: col[6].i = ex;
1957: col[6].j = ey;
1958: col[6].k = ez;
1959: col[6].loc = ELEMENT;
1960: col[6].c = 0;
1961: valA[6] = -1.0 / hx;
1962: } else {
1963: nEntries = 8;
1964: col[0].i = ex;
1965: col[0].j = ey;
1966: col[0].k = ez;
1967: col[0].loc = LEFT;
1968: col[0].c = 0;
1969: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1970: col[1].i = ex;
1971: col[1].j = ey - 1;
1972: col[1].k = ez;
1973: col[1].loc = LEFT;
1974: col[1].c = 0;
1975: valA[1] = 1.0 / (hy * hy);
1976: /* Missing up term */
1977: col[2].i = ex - 1;
1978: col[2].j = ey;
1979: col[2].k = ez;
1980: col[2].loc = LEFT;
1981: col[2].c = 0;
1982: valA[2] = 1.0 / (hx * hx);
1983: col[3].i = ex + 1;
1984: col[3].j = ey;
1985: col[3].k = ez;
1986: col[3].loc = LEFT;
1987: col[3].c = 0;
1988: valA[3] = 1.0 / (hx * hx);
1989: col[4].i = ex;
1990: col[4].j = ey;
1991: col[4].k = ez - 1;
1992: col[4].loc = LEFT;
1993: col[4].c = 0;
1994: valA[4] = 1.0 / (hz * hz);
1995: col[5].i = ex;
1996: col[5].j = ey;
1997: col[5].k = ez + 1;
1998: col[5].loc = LEFT;
1999: col[5].c = 0;
2000: valA[5] = 1.0 / (hz * hz);
2001: col[6].i = ex - 1;
2002: col[6].j = ey;
2003: col[6].k = ez;
2004: col[6].loc = ELEMENT;
2005: col[6].c = 0;
2006: valA[6] = 1.0 / hx;
2007: col[7].i = ex;
2008: col[7].j = ey;
2009: col[7].k = ez;
2010: col[7].loc = ELEMENT;
2011: col[7].c = 0;
2012: valA[7] = -1.0 / hx;
2013: }
2014: } else if (ez == 0) {
2015: nEntries = 8;
2016: col[0].i = ex;
2017: col[0].j = ey;
2018: col[0].k = ez;
2019: col[0].loc = LEFT;
2020: col[0].c = 0;
2021: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2022: col[1].i = ex;
2023: col[1].j = ey - 1;
2024: col[1].k = ez;
2025: col[1].loc = LEFT;
2026: col[1].c = 0;
2027: valA[1] = 1.0 / (hy * hy);
2028: col[2].i = ex;
2029: col[2].j = ey + 1;
2030: col[2].k = ez;
2031: col[2].loc = LEFT;
2032: col[2].c = 0;
2033: valA[2] = 1.0 / (hy * hy);
2034: col[3].i = ex - 1;
2035: col[3].j = ey;
2036: col[3].k = ez;
2037: col[3].loc = LEFT;
2038: col[3].c = 0;
2039: valA[3] = 1.0 / (hx * hx);
2040: col[4].i = ex + 1;
2041: col[4].j = ey;
2042: col[4].k = ez;
2043: col[4].loc = LEFT;
2044: col[4].c = 0;
2045: valA[4] = 1.0 / (hx * hx);
2046: /* Missing back term */
2047: col[5].i = ex;
2048: col[5].j = ey;
2049: col[5].k = ez + 1;
2050: col[5].loc = LEFT;
2051: col[5].c = 0;
2052: valA[5] = 1.0 / (hz * hz);
2053: col[6].i = ex - 1;
2054: col[6].j = ey;
2055: col[6].k = ez;
2056: col[6].loc = ELEMENT;
2057: col[6].c = 0;
2058: valA[6] = 1.0 / hx;
2059: col[7].i = ex;
2060: col[7].j = ey;
2061: col[7].k = ez;
2062: col[7].loc = ELEMENT;
2063: col[7].c = 0;
2064: valA[7] = -1.0 / hx;
2065: } else if (ez == N[2] - 1) {
2066: nEntries = 8;
2067: col[0].i = ex;
2068: col[0].j = ey;
2069: col[0].k = ez;
2070: col[0].loc = LEFT;
2071: col[0].c = 0;
2072: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2073: col[1].i = ex;
2074: col[1].j = ey - 1;
2075: col[1].k = ez;
2076: col[1].loc = LEFT;
2077: col[1].c = 0;
2078: valA[1] = 1.0 / (hy * hy);
2079: col[2].i = ex;
2080: col[2].j = ey + 1;
2081: col[2].k = ez;
2082: col[2].loc = LEFT;
2083: col[2].c = 0;
2084: valA[2] = 1.0 / (hy * hy);
2085: col[3].i = ex - 1;
2086: col[3].j = ey;
2087: col[3].k = ez;
2088: col[3].loc = LEFT;
2089: col[3].c = 0;
2090: valA[3] = 1.0 / (hx * hx);
2091: col[4].i = ex + 1;
2092: col[4].j = ey;
2093: col[4].k = ez;
2094: col[4].loc = LEFT;
2095: col[4].c = 0;
2096: valA[4] = 1.0 / (hx * hx);
2097: col[5].i = ex;
2098: col[5].j = ey;
2099: col[5].k = ez - 1;
2100: col[5].loc = LEFT;
2101: col[5].c = 0;
2102: valA[5] = 1.0 / (hz * hz);
2103: /* Missing front term */
2104: col[6].i = ex - 1;
2105: col[6].j = ey;
2106: col[6].k = ez;
2107: col[6].loc = ELEMENT;
2108: col[6].c = 0;
2109: valA[6] = 1.0 / hx;
2110: col[7].i = ex;
2111: col[7].j = ey;
2112: col[7].k = ez;
2113: col[7].loc = ELEMENT;
2114: col[7].c = 0;
2115: valA[7] = -1.0 / hx;
2116: } else {
2117: nEntries = 9;
2118: col[0].i = ex;
2119: col[0].j = ey;
2120: col[0].k = ez;
2121: col[0].loc = LEFT;
2122: col[0].c = 0;
2123: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2124: col[1].i = ex;
2125: col[1].j = ey - 1;
2126: col[1].k = ez;
2127: col[1].loc = LEFT;
2128: col[1].c = 0;
2129: valA[1] = 1.0 / (hy * hy);
2130: col[2].i = ex;
2131: col[2].j = ey + 1;
2132: col[2].k = ez;
2133: col[2].loc = LEFT;
2134: col[2].c = 0;
2135: valA[2] = 1.0 / (hy * hy);
2136: col[3].i = ex - 1;
2137: col[3].j = ey;
2138: col[3].k = ez;
2139: col[3].loc = LEFT;
2140: col[3].c = 0;
2141: valA[3] = 1.0 / (hx * hx);
2142: col[4].i = ex + 1;
2143: col[4].j = ey;
2144: col[4].k = ez;
2145: col[4].loc = LEFT;
2146: col[4].c = 0;
2147: valA[4] = 1.0 / (hx * hx);
2148: col[5].i = ex;
2149: col[5].j = ey;
2150: col[5].k = ez - 1;
2151: col[5].loc = LEFT;
2152: col[5].c = 0;
2153: valA[5] = 1.0 / (hz * hz);
2154: col[6].i = ex;
2155: col[6].j = ey;
2156: col[6].k = ez + 1;
2157: col[6].loc = LEFT;
2158: col[6].c = 0;
2159: valA[6] = 1.0 / (hz * hz);
2160: col[7].i = ex - 1;
2161: col[7].j = ey;
2162: col[7].k = ez;
2163: col[7].loc = ELEMENT;
2164: col[7].c = 0;
2165: valA[7] = 1.0 / hx;
2166: col[8].i = ex;
2167: col[8].j = ey;
2168: col[8].k = ez;
2169: col[8].loc = ELEMENT;
2170: col[8].c = 0;
2171: valA[8] = -1.0 / hx;
2172: }
2173: PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed));
2174: PetscCall(check_vals(ex, ey, ez, nEntries, valA, computed));
2175: }
2177: /* Equation on bottom face of this element */
2178: if (ey == 0) {
2179: /* Bottom boundary velocity Dirichlet */
2180: DMStagStencil row;
2181: const PetscScalar valA = 1.0;
2182: row.i = ex;
2183: row.j = ey;
2184: row.k = ez;
2185: row.loc = DOWN;
2186: row.c = 0;
2187: PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
2188: PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
2189: } else {
2190: /* Y-momentum equation, (v_xx + v_yy + v_zz) - p_y = f^y */
2191: DMStagStencil row, col[9];
2192: PetscScalar valA[9];
2193: PetscInt nEntries;
2195: row.i = ex;
2196: row.j = ey;
2197: row.k = ez;
2198: row.loc = DOWN;
2199: row.c = 0;
2200: if (ex == 0) {
2201: if (ez == 0) {
2202: nEntries = 7;
2203: col[0].i = ex;
2204: col[0].j = ey;
2205: col[0].k = ez;
2206: col[0].loc = DOWN;
2207: col[0].c = 0;
2208: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2209: col[1].i = ex;
2210: col[1].j = ey - 1;
2211: col[1].k = ez;
2212: col[1].loc = DOWN;
2213: col[1].c = 0;
2214: valA[1] = 1.0 / (hy * hy);
2215: col[2].i = ex;
2216: col[2].j = ey + 1;
2217: col[2].k = ez;
2218: col[2].loc = DOWN;
2219: col[2].c = 0;
2220: valA[2] = 1.0 / (hy * hy);
2221: /* Left term missing */
2222: col[3].i = ex + 1;
2223: col[3].j = ey;
2224: col[3].k = ez;
2225: col[3].loc = DOWN;
2226: col[3].c = 0;
2227: valA[3] = 1.0 / (hx * hx);
2228: /* Back term missing */
2229: col[4].i = ex;
2230: col[4].j = ey;
2231: col[4].k = ez + 1;
2232: col[4].loc = DOWN;
2233: col[4].c = 0;
2234: valA[4] = 1.0 / (hz * hz);
2235: col[5].i = ex;
2236: col[5].j = ey - 1;
2237: col[5].k = ez;
2238: col[5].loc = ELEMENT;
2239: col[5].c = 0;
2240: valA[5] = 1.0 / hy;
2241: col[6].i = ex;
2242: col[6].j = ey;
2243: col[6].k = ez;
2244: col[6].loc = ELEMENT;
2245: col[6].c = 0;
2246: valA[6] = -1.0 / hy;
2247: } else if (ez == N[2] - 1) {
2248: nEntries = 7;
2249: col[0].i = ex;
2250: col[0].j = ey;
2251: col[0].k = ez;
2252: col[0].loc = DOWN;
2253: col[0].c = 0;
2254: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2255: col[1].i = ex;
2256: col[1].j = ey - 1;
2257: col[1].k = ez;
2258: col[1].loc = DOWN;
2259: col[1].c = 0;
2260: valA[1] = 1.0 / (hy * hy);
2261: col[2].i = ex;
2262: col[2].j = ey + 1;
2263: col[2].k = ez;
2264: col[2].loc = DOWN;
2265: col[2].c = 0;
2266: valA[2] = 1.0 / (hy * hy);
2267: /* Left term missing */
2268: col[3].i = ex + 1;
2269: col[3].j = ey;
2270: col[3].k = ez;
2271: col[3].loc = DOWN;
2272: col[3].c = 0;
2273: valA[3] = 1.0 / (hx * hx);
2274: col[4].i = ex;
2275: col[4].j = ey;
2276: col[4].k = ez - 1;
2277: col[4].loc = DOWN;
2278: col[4].c = 0;
2279: valA[4] = 1.0 / (hz * hz);
2280: /* Front term missing */
2281: col[5].i = ex;
2282: col[5].j = ey - 1;
2283: col[5].k = ez;
2284: col[5].loc = ELEMENT;
2285: col[5].c = 0;
2286: valA[5] = 1.0 / hy;
2287: col[6].i = ex;
2288: col[6].j = ey;
2289: col[6].k = ez;
2290: col[6].loc = ELEMENT;
2291: col[6].c = 0;
2292: valA[6] = -1.0 / hy;
2293: } else {
2294: nEntries = 8;
2295: col[0].i = ex;
2296: col[0].j = ey;
2297: col[0].k = ez;
2298: col[0].loc = DOWN;
2299: col[0].c = 0;
2300: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2301: col[1].i = ex;
2302: col[1].j = ey - 1;
2303: col[1].k = ez;
2304: col[1].loc = DOWN;
2305: col[1].c = 0;
2306: valA[1] = 1.0 / (hy * hy);
2307: col[2].i = ex;
2308: col[2].j = ey + 1;
2309: col[2].k = ez;
2310: col[2].loc = DOWN;
2311: col[2].c = 0;
2312: valA[2] = 1.0 / (hy * hy);
2313: /* Left term missing */
2314: col[3].i = ex + 1;
2315: col[3].j = ey;
2316: col[3].k = ez;
2317: col[3].loc = DOWN;
2318: col[3].c = 0;
2319: valA[3] = 1.0 / (hx * hx);
2320: col[4].i = ex;
2321: col[4].j = ey;
2322: col[4].k = ez - 1;
2323: col[4].loc = DOWN;
2324: col[4].c = 0;
2325: valA[4] = 1.0 / (hz * hz);
2326: col[5].i = ex;
2327: col[5].j = ey;
2328: col[5].k = ez + 1;
2329: col[5].loc = DOWN;
2330: col[5].c = 0;
2331: valA[5] = 1.0 / (hz * hz);
2332: col[6].i = ex;
2333: col[6].j = ey - 1;
2334: col[6].k = ez;
2335: col[6].loc = ELEMENT;
2336: col[6].c = 0;
2337: valA[6] = 1.0 / hy;
2338: col[7].i = ex;
2339: col[7].j = ey;
2340: col[7].k = ez;
2341: col[7].loc = ELEMENT;
2342: col[7].c = 0;
2343: valA[7] = -1.0 / hy;
2344: }
2345: } else if (ex == N[0] - 1) {
2346: if (ez == 0) {
2347: nEntries = 7;
2348: col[0].i = ex;
2349: col[0].j = ey;
2350: col[0].k = ez;
2351: col[0].loc = DOWN;
2352: col[0].c = 0;
2353: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2354: col[1].i = ex;
2355: col[1].j = ey - 1;
2356: col[1].k = ez;
2357: col[1].loc = DOWN;
2358: col[1].c = 0;
2359: valA[1] = 1.0 / (hy * hy);
2360: col[2].i = ex;
2361: col[2].j = ey + 1;
2362: col[2].k = ez;
2363: col[2].loc = DOWN;
2364: col[2].c = 0;
2365: valA[2] = 1.0 / (hy * hy);
2366: col[3].i = ex - 1;
2367: col[3].j = ey;
2368: col[3].k = ez;
2369: col[3].loc = DOWN;
2370: col[3].c = 0;
2371: valA[3] = 1.0 / (hx * hx);
2372: /* Right term missing */
2373: /* Back term missing */
2374: col[4].i = ex;
2375: col[4].j = ey;
2376: col[4].k = ez + 1;
2377: col[4].loc = DOWN;
2378: col[4].c = 0;
2379: valA[4] = 1.0 / (hz * hz);
2380: col[5].i = ex;
2381: col[5].j = ey - 1;
2382: col[5].k = ez;
2383: col[5].loc = ELEMENT;
2384: col[5].c = 0;
2385: valA[5] = 1.0 / hy;
2386: col[6].i = ex;
2387: col[6].j = ey;
2388: col[6].k = ez;
2389: col[6].loc = ELEMENT;
2390: col[6].c = 0;
2391: valA[6] = -1.0 / hy;
2392: } else if (ez == N[2] - 1) {
2393: nEntries = 7;
2394: col[0].i = ex;
2395: col[0].j = ey;
2396: col[0].k = ez;
2397: col[0].loc = DOWN;
2398: col[0].c = 0;
2399: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2400: col[1].i = ex;
2401: col[1].j = ey - 1;
2402: col[1].k = ez;
2403: col[1].loc = DOWN;
2404: col[1].c = 0;
2405: valA[1] = 1.0 / (hy * hy);
2406: col[2].i = ex;
2407: col[2].j = ey + 1;
2408: col[2].k = ez;
2409: col[2].loc = DOWN;
2410: col[2].c = 0;
2411: valA[2] = 1.0 / (hy * hy);
2412: col[3].i = ex - 1;
2413: col[3].j = ey;
2414: col[3].k = ez;
2415: col[3].loc = DOWN;
2416: col[3].c = 0;
2417: valA[3] = 1.0 / (hx * hx);
2418: /* Right term missing */
2419: col[4].i = ex;
2420: col[4].j = ey;
2421: col[4].k = ez - 1;
2422: col[4].loc = DOWN;
2423: col[4].c = 0;
2424: valA[4] = 1.0 / (hz * hz);
2425: /* Front term missing */
2426: col[5].i = ex;
2427: col[5].j = ey - 1;
2428: col[5].k = ez;
2429: col[5].loc = ELEMENT;
2430: col[5].c = 0;
2431: valA[5] = 1.0 / hy;
2432: col[6].i = ex;
2433: col[6].j = ey;
2434: col[6].k = ez;
2435: col[6].loc = ELEMENT;
2436: col[6].c = 0;
2437: valA[6] = -1.0 / hy;
2438: } else {
2439: nEntries = 8;
2440: col[0].i = ex;
2441: col[0].j = ey;
2442: col[0].k = ez;
2443: col[0].loc = DOWN;
2444: col[0].c = 0;
2445: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2446: col[1].i = ex;
2447: col[1].j = ey - 1;
2448: col[1].k = ez;
2449: col[1].loc = DOWN;
2450: col[1].c = 0;
2451: valA[1] = 1.0 / (hy * hy);
2452: col[2].i = ex;
2453: col[2].j = ey + 1;
2454: col[2].k = ez;
2455: col[2].loc = DOWN;
2456: col[2].c = 0;
2457: valA[2] = 1.0 / (hy * hy);
2458: col[3].i = ex - 1;
2459: col[3].j = ey;
2460: col[3].k = ez;
2461: col[3].loc = DOWN;
2462: col[3].c = 0;
2463: valA[3] = 1.0 / (hx * hx);
2464: /* Right term missing */
2465: col[4].i = ex;
2466: col[4].j = ey;
2467: col[4].k = ez - 1;
2468: col[4].loc = DOWN;
2469: col[4].c = 0;
2470: valA[4] = 1.0 / (hz * hz);
2471: col[5].i = ex;
2472: col[5].j = ey;
2473: col[5].k = ez + 1;
2474: col[5].loc = DOWN;
2475: col[5].c = 0;
2476: valA[5] = 1.0 / (hz * hz);
2477: col[6].i = ex;
2478: col[6].j = ey - 1;
2479: col[6].k = ez;
2480: col[6].loc = ELEMENT;
2481: col[6].c = 0;
2482: valA[6] = 1.0 / hy;
2483: col[7].i = ex;
2484: col[7].j = ey;
2485: col[7].k = ez;
2486: col[7].loc = ELEMENT;
2487: col[7].c = 0;
2488: valA[7] = -1.0 / hy;
2489: }
2490: } else if (ez == 0) {
2491: nEntries = 8;
2492: col[0].i = ex;
2493: col[0].j = ey;
2494: col[0].k = ez;
2495: col[0].loc = DOWN;
2496: col[0].c = 0;
2497: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2498: col[1].i = ex;
2499: col[1].j = ey - 1;
2500: col[1].k = ez;
2501: col[1].loc = DOWN;
2502: col[1].c = 0;
2503: valA[1] = 1.0 / (hy * hy);
2504: col[2].i = ex;
2505: col[2].j = ey + 1;
2506: col[2].k = ez;
2507: col[2].loc = DOWN;
2508: col[2].c = 0;
2509: valA[2] = 1.0 / (hy * hy);
2510: col[3].i = ex - 1;
2511: col[3].j = ey;
2512: col[3].k = ez;
2513: col[3].loc = DOWN;
2514: col[3].c = 0;
2515: valA[3] = 1.0 / (hx * hx);
2516: col[4].i = ex + 1;
2517: col[4].j = ey;
2518: col[4].k = ez;
2519: col[4].loc = DOWN;
2520: col[4].c = 0;
2521: valA[4] = 1.0 / (hx * hx);
2522: /* Back term missing */
2523: col[5].i = ex;
2524: col[5].j = ey;
2525: col[5].k = ez + 1;
2526: col[5].loc = DOWN;
2527: col[5].c = 0;
2528: valA[5] = 1.0 / (hz * hz);
2529: col[6].i = ex;
2530: col[6].j = ey - 1;
2531: col[6].k = ez;
2532: col[6].loc = ELEMENT;
2533: col[6].c = 0;
2534: valA[6] = 1.0 / hy;
2535: col[7].i = ex;
2536: col[7].j = ey;
2537: col[7].k = ez;
2538: col[7].loc = ELEMENT;
2539: col[7].c = 0;
2540: valA[7] = -1.0 / hy;
2541: } else if (ez == N[2] - 1) {
2542: nEntries = 8;
2543: col[0].i = ex;
2544: col[0].j = ey;
2545: col[0].k = ez;
2546: col[0].loc = DOWN;
2547: col[0].c = 0;
2548: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2549: col[1].i = ex;
2550: col[1].j = ey - 1;
2551: col[1].k = ez;
2552: col[1].loc = DOWN;
2553: col[1].c = 0;
2554: valA[1] = 1.0 / (hy * hy);
2555: col[2].i = ex;
2556: col[2].j = ey + 1;
2557: col[2].k = ez;
2558: col[2].loc = DOWN;
2559: col[2].c = 0;
2560: valA[2] = 1.0 / (hy * hy);
2561: col[3].i = ex - 1;
2562: col[3].j = ey;
2563: col[3].k = ez;
2564: col[3].loc = DOWN;
2565: col[3].c = 0;
2566: valA[3] = 1.0 / (hx * hx);
2567: col[4].i = ex + 1;
2568: col[4].j = ey;
2569: col[4].k = ez;
2570: col[4].loc = DOWN;
2571: col[4].c = 0;
2572: valA[4] = 1.0 / (hx * hx);
2573: col[5].i = ex;
2574: col[5].j = ey;
2575: col[5].k = ez - 1;
2576: col[5].loc = DOWN;
2577: col[5].c = 0;
2578: valA[5] = 1.0 / (hz * hz);
2579: /* Front term missing */
2580: col[6].i = ex;
2581: col[6].j = ey - 1;
2582: col[6].k = ez;
2583: col[6].loc = ELEMENT;
2584: col[6].c = 0;
2585: valA[6] = 1.0 / hy;
2586: col[7].i = ex;
2587: col[7].j = ey;
2588: col[7].k = ez;
2589: col[7].loc = ELEMENT;
2590: col[7].c = 0;
2591: valA[7] = -1.0 / hy;
2592: } else {
2593: nEntries = 9;
2594: col[0].i = ex;
2595: col[0].j = ey;
2596: col[0].k = ez;
2597: col[0].loc = DOWN;
2598: col[0].c = 0;
2599: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2600: col[1].i = ex;
2601: col[1].j = ey - 1;
2602: col[1].k = ez;
2603: col[1].loc = DOWN;
2604: col[1].c = 0;
2605: valA[1] = 1.0 / (hy * hy);
2606: col[2].i = ex;
2607: col[2].j = ey + 1;
2608: col[2].k = ez;
2609: col[2].loc = DOWN;
2610: col[2].c = 0;
2611: valA[2] = 1.0 / (hy * hy);
2612: col[3].i = ex - 1;
2613: col[3].j = ey;
2614: col[3].k = ez;
2615: col[3].loc = DOWN;
2616: col[3].c = 0;
2617: valA[3] = 1.0 / (hx * hx);
2618: col[4].i = ex + 1;
2619: col[4].j = ey;
2620: col[4].k = ez;
2621: col[4].loc = DOWN;
2622: col[4].c = 0;
2623: valA[4] = 1.0 / (hx * hx);
2624: col[5].i = ex;
2625: col[5].j = ey;
2626: col[5].k = ez - 1;
2627: col[5].loc = DOWN;
2628: col[5].c = 0;
2629: valA[5] = 1.0 / (hz * hz);
2630: col[6].i = ex;
2631: col[6].j = ey;
2632: col[6].k = ez + 1;
2633: col[6].loc = DOWN;
2634: col[6].c = 0;
2635: valA[6] = 1.0 / (hz * hz);
2636: col[7].i = ex;
2637: col[7].j = ey - 1;
2638: col[7].k = ez;
2639: col[7].loc = ELEMENT;
2640: col[7].c = 0;
2641: valA[7] = 1.0 / hy;
2642: col[8].i = ex;
2643: col[8].j = ey;
2644: col[8].k = ez;
2645: col[8].loc = ELEMENT;
2646: col[8].c = 0;
2647: valA[8] = -1.0 / hy;
2648: }
2649: PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed));
2650: PetscCall(check_vals(ex, ey, ez, nEntries, valA, computed));
2651: }
2653: /* Equation on back face of this element */
2654: if (ez == 0) {
2655: /* Back boundary velocity Dirichlet */
2656: DMStagStencil row;
2657: const PetscScalar valA = 1.0;
2658: row.i = ex;
2659: row.j = ey;
2660: row.k = ez;
2661: row.loc = BACK;
2662: row.c = 0;
2663: PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
2664: PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
2665: } else {
2666: /* Z-momentum equation, (w_xx + w_yy + w_zz) - p_z = f^z */
2667: DMStagStencil row, col[9];
2668: PetscScalar valA[9];
2669: PetscInt nEntries;
2671: row.i = ex;
2672: row.j = ey;
2673: row.k = ez;
2674: row.loc = BACK;
2675: row.c = 0;
2676: if (ex == 0) {
2677: if (ey == 0) {
2678: nEntries = 7;
2679: col[0].i = ex;
2680: col[0].j = ey;
2681: col[0].k = ez;
2682: col[0].loc = BACK;
2683: col[0].c = 0;
2684: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2685: /* Down term missing */
2686: col[1].i = ex;
2687: col[1].j = ey + 1;
2688: col[1].k = ez;
2689: col[1].loc = BACK;
2690: col[1].c = 0;
2691: valA[1] = 1.0 / (hy * hy);
2692: /* Left term missing */
2693: col[2].i = ex + 1;
2694: col[2].j = ey;
2695: col[2].k = ez;
2696: col[2].loc = BACK;
2697: col[2].c = 0;
2698: valA[2] = 1.0 / (hx * hx);
2699: col[3].i = ex;
2700: col[3].j = ey;
2701: col[3].k = ez - 1;
2702: col[3].loc = BACK;
2703: col[3].c = 0;
2704: valA[3] = 1.0 / (hz * hz);
2705: col[4].i = ex;
2706: col[4].j = ey;
2707: col[4].k = ez + 1;
2708: col[4].loc = BACK;
2709: col[4].c = 0;
2710: valA[4] = 1.0 / (hz * hz);
2711: col[5].i = ex;
2712: col[5].j = ey;
2713: col[5].k = ez - 1;
2714: col[5].loc = ELEMENT;
2715: col[5].c = 0;
2716: valA[5] = 1.0 / hz;
2717: col[6].i = ex;
2718: col[6].j = ey;
2719: col[6].k = ez;
2720: col[6].loc = ELEMENT;
2721: col[6].c = 0;
2722: valA[6] = -1.0 / hz;
2723: } else if (ey == N[1] - 1) {
2724: nEntries = 7;
2725: col[0].i = ex;
2726: col[0].j = ey;
2727: col[0].k = ez;
2728: col[0].loc = BACK;
2729: col[0].c = 0;
2730: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2731: col[1].i = ex;
2732: col[1].j = ey - 1;
2733: col[1].k = ez;
2734: col[1].loc = BACK;
2735: col[1].c = 0;
2736: valA[1] = 1.0 / (hy * hy);
2737: /* Up term missing */
2738: /* Left term missing */
2739: col[2].i = ex + 1;
2740: col[2].j = ey;
2741: col[2].k = ez;
2742: col[2].loc = BACK;
2743: col[2].c = 0;
2744: valA[2] = 1.0 / (hx * hx);
2745: col[3].i = ex;
2746: col[3].j = ey;
2747: col[3].k = ez - 1;
2748: col[3].loc = BACK;
2749: col[3].c = 0;
2750: valA[3] = 1.0 / (hz * hz);
2751: col[4].i = ex;
2752: col[4].j = ey;
2753: col[4].k = ez + 1;
2754: col[4].loc = BACK;
2755: col[4].c = 0;
2756: valA[4] = 1.0 / (hz * hz);
2757: col[5].i = ex;
2758: col[5].j = ey;
2759: col[5].k = ez - 1;
2760: col[5].loc = ELEMENT;
2761: col[5].c = 0;
2762: valA[5] = 1.0 / hz;
2763: col[6].i = ex;
2764: col[6].j = ey;
2765: col[6].k = ez;
2766: col[6].loc = ELEMENT;
2767: col[6].c = 0;
2768: valA[6] = -1.0 / hz;
2769: } else {
2770: nEntries = 8;
2771: col[0].i = ex;
2772: col[0].j = ey;
2773: col[0].k = ez;
2774: col[0].loc = BACK;
2775: col[0].c = 0;
2776: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2777: col[1].i = ex;
2778: col[1].j = ey - 1;
2779: col[1].k = ez;
2780: col[1].loc = BACK;
2781: col[1].c = 0;
2782: valA[1] = 1.0 / (hy * hy);
2783: col[2].i = ex;
2784: col[2].j = ey + 1;
2785: col[2].k = ez;
2786: col[2].loc = BACK;
2787: col[2].c = 0;
2788: valA[2] = 1.0 / (hy * hy);
2789: /* Left term missing */
2790: col[3].i = ex + 1;
2791: col[3].j = ey;
2792: col[3].k = ez;
2793: col[3].loc = BACK;
2794: col[3].c = 0;
2795: valA[3] = 1.0 / (hx * hx);
2796: col[4].i = ex;
2797: col[4].j = ey;
2798: col[4].k = ez - 1;
2799: col[4].loc = BACK;
2800: col[4].c = 0;
2801: valA[4] = 1.0 / (hz * hz);
2802: col[5].i = ex;
2803: col[5].j = ey;
2804: col[5].k = ez + 1;
2805: col[5].loc = BACK;
2806: col[5].c = 0;
2807: valA[5] = 1.0 / (hz * hz);
2808: col[6].i = ex;
2809: col[6].j = ey;
2810: col[6].k = ez - 1;
2811: col[6].loc = ELEMENT;
2812: col[6].c = 0;
2813: valA[6] = 1.0 / hz;
2814: col[7].i = ex;
2815: col[7].j = ey;
2816: col[7].k = ez;
2817: col[7].loc = ELEMENT;
2818: col[7].c = 0;
2819: valA[7] = -1.0 / hz;
2820: }
2821: } else if (ex == N[0] - 1) {
2822: if (ey == 0) {
2823: nEntries = 7;
2824: col[0].i = ex;
2825: col[0].j = ey;
2826: col[0].k = ez;
2827: col[0].loc = BACK;
2828: col[0].c = 0;
2829: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2830: /* Down term missing */
2831: col[1].i = ex;
2832: col[1].j = ey + 1;
2833: col[1].k = ez;
2834: col[1].loc = BACK;
2835: col[1].c = 0;
2836: valA[1] = 1.0 / (hy * hy);
2837: col[2].i = ex - 1;
2838: col[2].j = ey;
2839: col[2].k = ez;
2840: col[2].loc = BACK;
2841: col[2].c = 0;
2842: valA[2] = 1.0 / (hx * hx);
2843: /* Right term missing */
2844: col[3].i = ex;
2845: col[3].j = ey;
2846: col[3].k = ez - 1;
2847: col[3].loc = BACK;
2848: col[3].c = 0;
2849: valA[3] = 1.0 / (hz * hz);
2850: col[4].i = ex;
2851: col[4].j = ey;
2852: col[4].k = ez + 1;
2853: col[4].loc = BACK;
2854: col[4].c = 0;
2855: valA[4] = 1.0 / (hz * hz);
2856: col[5].i = ex;
2857: col[5].j = ey;
2858: col[5].k = ez - 1;
2859: col[5].loc = ELEMENT;
2860: col[5].c = 0;
2861: valA[5] = 1.0 / hz;
2862: col[6].i = ex;
2863: col[6].j = ey;
2864: col[6].k = ez;
2865: col[6].loc = ELEMENT;
2866: col[6].c = 0;
2867: valA[6] = -1.0 / hz;
2868: } else if (ey == N[1] - 1) {
2869: nEntries = 7;
2870: col[0].i = ex;
2871: col[0].j = ey;
2872: col[0].k = ez;
2873: col[0].loc = BACK;
2874: col[0].c = 0;
2875: valA[0] = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2876: col[1].i = ex;
2877: col[1].j = ey - 1;
2878: col[1].k = ez;
2879: col[1].loc = BACK;
2880: col[1].c = 0;
2881: valA[1] = 1.0 / (hy * hy);
2882: /* Up term missing */
2883: col[2].i = ex - 1;
2884: col[2].j = ey;
2885: col[2].k = ez;
2886: col[2].loc = BACK;
2887: col[2].c = 0;
2888: valA[2] = 1.0 / (hx * hx);
2889: /* Right term missing */
2890: col[3].i = ex;
2891: col[3].j = ey;
2892: col[3].k = ez - 1;
2893: col[3].loc = BACK;
2894: col[3].c = 0;
2895: valA[3] = 1.0 / (hz * hz);
2896: col[4].i = ex;
2897: col[4].j = ey;
2898: col[4].k = ez + 1;
2899: col[4].loc = BACK;
2900: col[4].c = 0;
2901: valA[4] = 1.0 / (hz * hz);
2902: col[5].i = ex;
2903: col[5].j = ey;
2904: col[5].k = ez - 1;
2905: col[5].loc = ELEMENT;
2906: col[5].c = 0;
2907: valA[5] = 1.0 / hz;
2908: col[6].i = ex;
2909: col[6].j = ey;
2910: col[6].k = ez;
2911: col[6].loc = ELEMENT;
2912: col[6].c = 0;
2913: valA[6] = -1.0 / hz;
2914: } else {
2915: nEntries = 8;
2916: col[0].i = ex;
2917: col[0].j = ey;
2918: col[0].k = ez;
2919: col[0].loc = BACK;
2920: col[0].c = 0;
2921: valA[0] = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2922: col[1].i = ex;
2923: col[1].j = ey - 1;
2924: col[1].k = ez;
2925: col[1].loc = BACK;
2926: col[1].c = 0;
2927: valA[1] = 1.0 / (hy * hy);
2928: col[2].i = ex;
2929: col[2].j = ey + 1;
2930: col[2].k = ez;
2931: col[2].loc = BACK;
2932: col[2].c = 0;
2933: valA[2] = 1.0 / (hy * hy);
2934: col[3].i = ex - 1;
2935: col[3].j = ey;
2936: col[3].k = ez;
2937: col[3].loc = BACK;
2938: col[3].c = 0;
2939: valA[3] = 1.0 / (hx * hx);
2940: /* Right term missing */
2941: col[4].i = ex;
2942: col[4].j = ey;
2943: col[4].k = ez - 1;
2944: col[4].loc = BACK;
2945: col[4].c = 0;
2946: valA[4] = 1.0 / (hz * hz);
2947: col[5].i = ex;
2948: col[5].j = ey;
2949: col[5].k = ez + 1;
2950: col[5].loc = BACK;
2951: col[5].c = 0;
2952: valA[5] = 1.0 / (hz * hz);
2953: col[6].i = ex;
2954: col[6].j = ey;
2955: col[6].k = ez - 1;
2956: col[6].loc = ELEMENT;
2957: col[6].c = 0;
2958: valA[6] = 1.0 / hz;
2959: col[7].i = ex;
2960: col[7].j = ey;
2961: col[7].k = ez;
2962: col[7].loc = ELEMENT;
2963: col[7].c = 0;
2964: valA[7] = -1.0 / hz;
2965: }
2966: } else if (ey == 0) {
2967: nEntries = 8;
2968: col[0].i = ex;
2969: col[0].j = ey;
2970: col[0].k = ez;
2971: col[0].loc = BACK;
2972: col[0].c = 0;
2973: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2974: /* Down term missing */
2975: col[1].i = ex;
2976: col[1].j = ey + 1;
2977: col[1].k = ez;
2978: col[1].loc = BACK;
2979: col[1].c = 0;
2980: valA[1] = 1.0 / (hy * hy);
2981: col[2].i = ex - 1;
2982: col[2].j = ey;
2983: col[2].k = ez;
2984: col[2].loc = BACK;
2985: col[2].c = 0;
2986: valA[2] = 1.0 / (hx * hx);
2987: col[3].i = ex + 1;
2988: col[3].j = ey;
2989: col[3].k = ez;
2990: col[3].loc = BACK;
2991: col[3].c = 0;
2992: valA[3] = 1.0 / (hx * hx);
2993: col[4].i = ex;
2994: col[4].j = ey;
2995: col[4].k = ez - 1;
2996: col[4].loc = BACK;
2997: col[4].c = 0;
2998: valA[4] = 1.0 / (hz * hz);
2999: col[5].i = ex;
3000: col[5].j = ey;
3001: col[5].k = ez + 1;
3002: col[5].loc = BACK;
3003: col[5].c = 0;
3004: valA[5] = 1.0 / (hz * hz);
3005: col[6].i = ex;
3006: col[6].j = ey;
3007: col[6].k = ez - 1;
3008: col[6].loc = ELEMENT;
3009: col[6].c = 0;
3010: valA[6] = 1.0 / hz;
3011: col[7].i = ex;
3012: col[7].j = ey;
3013: col[7].k = ez;
3014: col[7].loc = ELEMENT;
3015: col[7].c = 0;
3016: valA[7] = -1.0 / hz;
3017: } else if (ey == N[1] - 1) {
3018: nEntries = 8;
3019: col[0].i = ex;
3020: col[0].j = ey;
3021: col[0].k = ez;
3022: col[0].loc = BACK;
3023: col[0].c = 0;
3024: valA[0] = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
3025: col[1].i = ex;
3026: col[1].j = ey - 1;
3027: col[1].k = ez;
3028: col[1].loc = BACK;
3029: col[1].c = 0;
3030: valA[1] = 1.0 / (hy * hy);
3031: /* Up term missing */
3032: col[2].i = ex - 1;
3033: col[2].j = ey;
3034: col[2].k = ez;
3035: col[2].loc = BACK;
3036: col[2].c = 0;
3037: valA[2] = 1.0 / (hx * hx);
3038: col[3].i = ex + 1;
3039: col[3].j = ey;
3040: col[3].k = ez;
3041: col[3].loc = BACK;
3042: col[3].c = 0;
3043: valA[3] = 1.0 / (hx * hx);
3044: col[4].i = ex;
3045: col[4].j = ey;
3046: col[4].k = ez - 1;
3047: col[4].loc = BACK;
3048: col[4].c = 0;
3049: valA[4] = 1.0 / (hz * hz);
3050: col[5].i = ex;
3051: col[5].j = ey;
3052: col[5].k = ez + 1;
3053: col[5].loc = BACK;
3054: col[5].c = 0;
3055: valA[5] = 1.0 / (hz * hz);
3056: col[6].i = ex;
3057: col[6].j = ey;
3058: col[6].k = ez - 1;
3059: col[6].loc = ELEMENT;
3060: col[6].c = 0;
3061: valA[6] = 1.0 / hz;
3062: col[7].i = ex;
3063: col[7].j = ey;
3064: col[7].k = ez;
3065: col[7].loc = ELEMENT;
3066: col[7].c = 0;
3067: valA[7] = -1.0 / hz;
3068: } else {
3069: nEntries = 9;
3070: col[0].i = ex;
3071: col[0].j = ey;
3072: col[0].k = ez;
3073: col[0].loc = BACK;
3074: col[0].c = 0;
3075: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
3076: col[1].i = ex;
3077: col[1].j = ey - 1;
3078: col[1].k = ez;
3079: col[1].loc = BACK;
3080: col[1].c = 0;
3081: valA[1] = 1.0 / (hy * hy);
3082: col[2].i = ex;
3083: col[2].j = ey + 1;
3084: col[2].k = ez;
3085: col[2].loc = BACK;
3086: col[2].c = 0;
3087: valA[2] = 1.0 / (hy * hy);
3088: col[3].i = ex - 1;
3089: col[3].j = ey;
3090: col[3].k = ez;
3091: col[3].loc = BACK;
3092: col[3].c = 0;
3093: valA[3] = 1.0 / (hx * hx);
3094: col[4].i = ex + 1;
3095: col[4].j = ey;
3096: col[4].k = ez;
3097: col[4].loc = BACK;
3098: col[4].c = 0;
3099: valA[4] = 1.0 / (hx * hx);
3100: col[5].i = ex;
3101: col[5].j = ey;
3102: col[5].k = ez - 1;
3103: col[5].loc = BACK;
3104: col[5].c = 0;
3105: valA[5] = 1.0 / (hz * hz);
3106: col[6].i = ex;
3107: col[6].j = ey;
3108: col[6].k = ez + 1;
3109: col[6].loc = BACK;
3110: col[6].c = 0;
3111: valA[6] = 1.0 / (hz * hz);
3112: col[7].i = ex;
3113: col[7].j = ey;
3114: col[7].k = ez - 1;
3115: col[7].loc = ELEMENT;
3116: col[7].c = 0;
3117: valA[7] = 1.0 / hz;
3118: col[8].i = ex;
3119: col[8].j = ey;
3120: col[8].k = ez;
3121: col[8].loc = ELEMENT;
3122: col[8].c = 0;
3123: valA[8] = -1.0 / hz;
3124: }
3125: PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed));
3126: PetscCall(check_vals(ex, ey, ez, nEntries, valA, computed));
3127: }
3129: /* P equation : u_x + v_y + w_z = g
3130: Note that this includes an explicit zero on the diagonal. This is only needed for
3131: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
3132: {
3133: DMStagStencil row, col[7];
3134: PetscScalar valA[7];
3136: row.i = ex;
3137: row.j = ey;
3138: row.k = ez;
3139: row.loc = ELEMENT;
3140: row.c = 0;
3141: col[0].i = ex;
3142: col[0].j = ey;
3143: col[0].k = ez;
3144: col[0].loc = LEFT;
3145: col[0].c = 0;
3146: valA[0] = -1.0 / hx;
3147: col[1].i = ex;
3148: col[1].j = ey;
3149: col[1].k = ez;
3150: col[1].loc = RIGHT;
3151: col[1].c = 0;
3152: valA[1] = 1.0 / hx;
3153: col[2].i = ex;
3154: col[2].j = ey;
3155: col[2].k = ez;
3156: col[2].loc = DOWN;
3157: col[2].c = 0;
3158: valA[2] = -1.0 / hy;
3159: col[3].i = ex;
3160: col[3].j = ey;
3161: col[3].k = ez;
3162: col[3].loc = UP;
3163: col[3].c = 0;
3164: valA[3] = 1.0 / hy;
3165: col[4].i = ex;
3166: col[4].j = ey;
3167: col[4].k = ez;
3168: col[4].loc = BACK;
3169: col[4].c = 0;
3170: valA[4] = -1.0 / hz;
3171: col[5].i = ex;
3172: col[5].j = ey;
3173: col[5].k = ez;
3174: col[5].loc = FRONT;
3175: col[5].c = 0;
3176: valA[5] = 1.0 / hz;
3177: col[6] = row;
3178: valA[6] = 0.0;
3179: PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 7, col, computed));
3180: PetscCall(check_vals(ex, ey, ez, 7, valA, computed));
3181: }
3182: }
3183: }
3184: }
3185: PetscCall(DMStagVecRestoreArrayRead(dmCoord, coordLocal, &arrCoord));
3186: PetscFunctionReturn(PETSC_SUCCESS);
3187: }
3189: /*TEST
3191: test:
3192: suffix: 1
3193: nsize: 1
3194: output_file: output/empty.out
3196: test:
3197: suffix: 2
3198: nsize: 4
3199: output_file: output/empty.out
3201: TEST*/