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, (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:     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));

1615:   PetscFunctionReturn(PETSC_SUCCESS);
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;

1623:   PetscFunctionBeginUser;
1624:   for (i = 0; i < n; ++i) {
1625:     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);
1626:   }
1627:   PetscFunctionReturn(PETSC_SUCCESS);
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];

1642:   PetscFunctionBeginUser;
1643:   PetscCall(DMStagGetCorners(dmSol, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
1644:   PetscCall(DMStagGetGlobalSizes(dmSol, &N[0], &N[1], &N[2]));
1645:   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");
1646:   PetscCall(DMStagGetIsLastRank(dmSol, &isLastRankx, &isLastRanky, &isLastRankz));
1647:   PetscCall(DMStagGetIsFirstRank(dmSol, &isFirstRankx, &isFirstRanky, &isFirstRankz));
1648:   hx = 1.0 / N[0];
1649:   hy = 1.0 / N[1];
1650:   hz = 1.0 / N[2];
1651:   PetscCall(DMGetCoordinateDM(dmSol, &dmCoord));
1652:   PetscCall(DMGetCoordinatesLocal(dmSol, &coordLocal));
1653:   PetscCall(DMStagVecGetArrayRead(dmCoord, coordLocal, &arrCoord));
1654:   for (d = 0; d < 3; ++d) {
1655:     PetscCall(DMStagGetLocationSlot(dmCoord, ELEMENT, d, &icp[d]));
1656:     PetscCall(DMStagGetLocationSlot(dmCoord, LEFT, d, &icux[d]));
1657:     PetscCall(DMStagGetLocationSlot(dmCoord, DOWN, d, &icuy[d]));
1658:     PetscCall(DMStagGetLocationSlot(dmCoord, BACK, d, &icuz[d]));
1659:     PetscCall(DMStagGetLocationSlot(dmCoord, RIGHT, d, &icux_right[d]));
1660:     PetscCall(DMStagGetLocationSlot(dmCoord, UP, d, &icuy_up[d]));
1661:     PetscCall(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:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1677:           PetscCall(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:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1689:           PetscCall(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:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1701:           PetscCall(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:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1715:           PetscCall(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:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed));
2177:           PetscCall(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:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
2191:           PetscCall(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:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed));
2653:           PetscCall(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:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
2667:           PetscCall(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:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed));
3129:           PetscCall(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:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 7, col, computed));
3183:           PetscCall(check_vals(ex, ey, ez, 7, valA, computed));
3184:         }
3185:       }
3186:     }
3187:   }
3188:   PetscCall(DMStagVecRestoreArrayRead(dmCoord, coordLocal, &arrCoord));
3189:   PetscFunctionReturn(PETSC_SUCCESS);
3190: }

3192: /*TEST

3194:    test:
3195:       suffix: 1
3196:       nsize: 1

3198:    test:
3199:       suffix: 2
3200:       nsize: 4

3202: TEST*/