Actual source code: ex30.c

  1: static char help[] = "Test DMStagMatGetValuesStencil() in 3D\n\n";

  3: #include <petscdm.h>
  4: #include <petscksp.h>
  5: #include <petscdmstag.h>

  7: /* Shorter, more convenient names for DMStagStencilLocation entries */
  8: #define BACK_DOWN_LEFT   DMSTAG_BACK_DOWN_LEFT
  9: #define BACK_DOWN        DMSTAG_BACK_DOWN
 10: #define BACK_DOWN_RIGHT  DMSTAG_BACK_DOWN_RIGHT
 11: #define BACK_LEFT        DMSTAG_BACK_LEFT
 12: #define BACK             DMSTAG_BACK
 13: #define BACK_RIGHT       DMSTAG_BACK_RIGHT
 14: #define BACK_UP_LEFT     DMSTAG_BACK_UP_LEFT
 15: #define BACK_UP          DMSTAG_BACK_UP
 16: #define BACK_UP_RIGHT    DMSTAG_BACK_UP_RIGHT
 17: #define DOWN_LEFT        DMSTAG_DOWN_LEFT
 18: #define DOWN             DMSTAG_DOWN
 19: #define DOWN_RIGHT       DMSTAG_DOWN_RIGHT
 20: #define LEFT             DMSTAG_LEFT
 21: #define ELEMENT          DMSTAG_ELEMENT
 22: #define RIGHT            DMSTAG_RIGHT
 23: #define UP_LEFT          DMSTAG_UP_LEFT
 24: #define UP               DMSTAG_UP
 25: #define UP_RIGHT         DMSTAG_UP_RIGHT
 26: #define FRONT_DOWN_LEFT  DMSTAG_FRONT_DOWN_LEFT
 27: #define FRONT_DOWN       DMSTAG_FRONT_DOWN
 28: #define FRONT_DOWN_RIGHT DMSTAG_FRONT_DOWN_RIGHT
 29: #define FRONT_LEFT       DMSTAG_FRONT_LEFT
 30: #define FRONT            DMSTAG_FRONT
 31: #define FRONT_RIGHT      DMSTAG_FRONT_RIGHT
 32: #define FRONT_UP_LEFT    DMSTAG_FRONT_UP_LEFT
 33: #define FRONT_UP         DMSTAG_FRONT_UP
 34: #define FRONT_UP_RIGHT   DMSTAG_FRONT_UP_RIGHT

 36: static PetscErrorCode CreateMat(DM, Mat *);
 37: static PetscErrorCode CheckMat(DM, Mat);

 39: int main(int argc, char **argv)
 40: {
 41:   DM  dmSol;
 42:   Mat A;

 44:   PetscFunctionBeginUser;
 45:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 46:   {
 47:     const PetscInt dof0 = 0, dof1 = 0, dof2 = 1, dof3 = 1; /* 1 dof on each face and element center */
 48:     const PetscInt stencilWidth = 1;
 49:     PetscCall(DMStagCreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 4, 5, 6, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, dof3, DMSTAG_STENCIL_BOX, stencilWidth, NULL, NULL, NULL, &dmSol));
 50:     PetscCall(DMSetFromOptions(dmSol));
 51:     PetscCall(DMSetUp(dmSol));
 52:     PetscCall(DMStagSetUniformCoordinatesExplicit(dmSol, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
 53:   }
 54:   PetscCall(CreateMat(dmSol, &A));
 55:   PetscCall(CheckMat(dmSol, A));
 56:   PetscCall(MatDestroy(&A));
 57:   PetscCall(DMDestroy(&dmSol));
 58:   PetscCall(PetscFinalize());
 59:   return 0;
 60: }

 62: static PetscErrorCode CreateMat(DM dmSol, Mat *pA)
 63: {
 64:   Vec             coordLocal;
 65:   Mat             A;
 66:   PetscInt        startx, starty, startz, N[3], nx, ny, nz, ex, ey, ez, d;
 67:   PetscInt        icp[3], icux[3], icuy[3], icuz[3], icux_right[3], icuy_up[3], icuz_front[3];
 68:   PetscBool       isLastRankx, isLastRanky, isLastRankz, isFirstRankx, isFirstRanky, isFirstRankz;
 69:   PetscReal       hx, hy, hz;
 70:   DM              dmCoord;
 71:   PetscScalar ****arrCoord;

 73:   PetscFunctionBeginUser;
 74:   PetscCall(DMCreateMatrix(dmSol, pA));
 75:   A = *pA;
 76:   PetscCall(DMStagGetCorners(dmSol, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
 77:   PetscCall(DMStagGetGlobalSizes(dmSol, &N[0], &N[1], &N[2]));
 78:   PetscCheck(N[0] >= 2 && N[1] >= 2 && N[2] >= 2, PetscObjectComm((PetscObject)dmSol), PETSC_ERR_ARG_SIZ, "This example requires at least two elements in each dimensions");
 79:   PetscCall(DMStagGetIsLastRank(dmSol, &isLastRankx, &isLastRanky, &isLastRankz));
 80:   PetscCall(DMStagGetIsFirstRank(dmSol, &isFirstRankx, &isFirstRanky, &isFirstRankz));
 81:   hx = 1.0 / N[0];
 82:   hy = 1.0 / N[1];
 83:   hz = 1.0 / N[2];
 84:   PetscCall(DMGetCoordinateDM(dmSol, &dmCoord));
 85:   PetscCall(DMGetCoordinatesLocal(dmSol, &coordLocal));
 86:   PetscCall(DMStagVecGetArrayRead(dmCoord, coordLocal, &arrCoord));
 87:   for (d = 0; d < 3; ++d) {
 88:     PetscCall(DMStagGetLocationSlot(dmCoord, ELEMENT, d, &icp[d]));
 89:     PetscCall(DMStagGetLocationSlot(dmCoord, LEFT, d, &icux[d]));
 90:     PetscCall(DMStagGetLocationSlot(dmCoord, DOWN, d, &icuy[d]));
 91:     PetscCall(DMStagGetLocationSlot(dmCoord, BACK, d, &icuz[d]));
 92:     PetscCall(DMStagGetLocationSlot(dmCoord, RIGHT, d, &icux_right[d]));
 93:     PetscCall(DMStagGetLocationSlot(dmCoord, UP, d, &icuy_up[d]));
 94:     PetscCall(DMStagGetLocationSlot(dmCoord, FRONT, d, &icuz_front[d]));
 95:   }

 97:   for (ez = startz; ez < startz + nz; ++ez) {
 98:     for (ey = starty; ey < starty + ny; ++ey) {
 99:       for (ex = startx; ex < startx + nx; ++ex) {
100:         if (ex == N[0] - 1) {
101:           /* Right Boundary velocity Dirichlet */
102:           DMStagStencil     row;
103:           const PetscScalar valA = 1.0;
104:           row.i                  = ex;
105:           row.j                  = ey;
106:           row.k                  = ez;
107:           row.loc                = RIGHT;
108:           row.c                  = 0;
109:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
110:         }
111:         if (ey == N[1] - 1) {
112:           /* Top boundary velocity Dirichlet */
113:           DMStagStencil     row;
114:           const PetscScalar valA = 1.0;
115:           row.i                  = ex;
116:           row.j                  = ey;
117:           row.k                  = ez;
118:           row.loc                = UP;
119:           row.c                  = 0;
120:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
121:         }
122:         if (ez == N[2] - 1) {
123:           /* Top boundary velocity Dirichlet */
124:           DMStagStencil     row;
125:           const PetscScalar valA = 1.0;
126:           row.i                  = ex;
127:           row.j                  = ey;
128:           row.k                  = ez;
129:           row.loc                = FRONT;
130:           row.c                  = 0;
131:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
132:         }

134:         /* Equation on left face of this element */
135:         if (ex == 0) {
136:           /* Left velocity Dirichlet */
137:           DMStagStencil     row;
138:           const PetscScalar valA = 1.0;
139:           row.i                  = ex;
140:           row.j                  = ey;
141:           row.k                  = ez;
142:           row.loc                = LEFT;
143:           row.c                  = 0;
144:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
145:         } else {
146:           /* X-momentum interior equation : (u_xx + u_yy + u_zz) - p_x = f^x */
147:           DMStagStencil row, col[9];
148:           PetscScalar   valA[9];
149:           PetscInt      nEntries;

151:           row.i   = ex;
152:           row.j   = ey;
153:           row.k   = ez;
154:           row.loc = LEFT;
155:           row.c   = 0;
156:           if (ey == 0) {
157:             if (ez == 0) {
158:               nEntries   = 7;
159:               col[0].i   = ex;
160:               col[0].j   = ey;
161:               col[0].k   = ez;
162:               col[0].loc = LEFT;
163:               col[0].c   = 0;
164:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
165:               /* Missing down term */
166:               col[1].i   = ex;
167:               col[1].j   = ey + 1;
168:               col[1].k   = ez;
169:               col[1].loc = LEFT;
170:               col[1].c   = 0;
171:               valA[1]    = 1.0 / (hy * hy);
172:               col[2].i   = ex - 1;
173:               col[2].j   = ey;
174:               col[2].k   = ez;
175:               col[2].loc = LEFT;
176:               col[2].c   = 0;
177:               valA[2]    = 1.0 / (hx * hx);
178:               col[3].i   = ex + 1;
179:               col[3].j   = ey;
180:               col[3].k   = ez;
181:               col[3].loc = LEFT;
182:               col[3].c   = 0;
183:               valA[3]    = 1.0 / (hx * hx);
184:               /* Missing back term */
185:               col[4].i   = ex;
186:               col[4].j   = ey;
187:               col[4].k   = ez + 1;
188:               col[4].loc = LEFT;
189:               col[4].c   = 0;
190:               valA[4]    = 1.0 / (hz * hz);
191:               col[5].i   = ex - 1;
192:               col[5].j   = ey;
193:               col[5].k   = ez;
194:               col[5].loc = ELEMENT;
195:               col[5].c   = 0;
196:               valA[5]    = 1.0 / hx;
197:               col[6].i   = ex;
198:               col[6].j   = ey;
199:               col[6].k   = ez;
200:               col[6].loc = ELEMENT;
201:               col[6].c   = 0;
202:               valA[6]    = -1.0 / hx;
203:             } else if (ez == N[2] - 1) {
204:               nEntries   = 7;
205:               col[0].i   = ex;
206:               col[0].j   = ey;
207:               col[0].k   = ez;
208:               col[0].loc = LEFT;
209:               col[0].c   = 0;
210:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
211:               /* Missing down term */
212:               col[1].i   = ex;
213:               col[1].j   = ey + 1;
214:               col[1].k   = ez;
215:               col[1].loc = LEFT;
216:               col[1].c   = 0;
217:               valA[1]    = 1.0 / (hy * hy);
218:               col[2].i   = ex - 1;
219:               col[2].j   = ey;
220:               col[2].k   = ez;
221:               col[2].loc = LEFT;
222:               col[2].c   = 0;
223:               valA[2]    = 1.0 / (hx * hx);
224:               col[3].i   = ex + 1;
225:               col[3].j   = ey;
226:               col[3].k   = ez;
227:               col[3].loc = LEFT;
228:               col[3].c   = 0;
229:               valA[3]    = 1.0 / (hx * hx);
230:               col[4].i   = ex;
231:               col[4].j   = ey;
232:               col[4].k   = ez - 1;
233:               col[4].loc = LEFT;
234:               col[4].c   = 0;
235:               valA[4]    = 1.0 / (hz * hz);
236:               /* Missing front term */
237:               col[5].i   = ex - 1;
238:               col[5].j   = ey;
239:               col[5].k   = ez;
240:               col[5].loc = ELEMENT;
241:               col[5].c   = 0;
242:               valA[5]    = 1.0 / hx;
243:               col[6].i   = ex;
244:               col[6].j   = ey;
245:               col[6].k   = ez;
246:               col[6].loc = ELEMENT;
247:               col[6].c   = 0;
248:               valA[6]    = -1.0 / hx;
249:             } else {
250:               nEntries   = 8;
251:               col[0].i   = ex;
252:               col[0].j   = ey;
253:               col[0].k   = ez;
254:               col[0].loc = LEFT;
255:               col[0].c   = 0;
256:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
257:               /* Missing down term */
258:               col[1].i   = ex;
259:               col[1].j   = ey + 1;
260:               col[1].k   = ez;
261:               col[1].loc = LEFT;
262:               col[1].c   = 0;
263:               valA[1]    = 1.0 / (hy * hy);
264:               col[2].i   = ex - 1;
265:               col[2].j   = ey;
266:               col[2].k   = ez;
267:               col[2].loc = LEFT;
268:               col[2].c   = 0;
269:               valA[2]    = 1.0 / (hx * hx);
270:               col[3].i   = ex + 1;
271:               col[3].j   = ey;
272:               col[3].k   = ez;
273:               col[3].loc = LEFT;
274:               col[3].c   = 0;
275:               valA[3]    = 1.0 / (hx * hx);
276:               col[4].i   = ex;
277:               col[4].j   = ey;
278:               col[4].k   = ez - 1;
279:               col[4].loc = LEFT;
280:               col[4].c   = 0;
281:               valA[4]    = 1.0 / (hz * hz);
282:               col[5].i   = ex;
283:               col[5].j   = ey;
284:               col[5].k   = ez + 1;
285:               col[5].loc = LEFT;
286:               col[5].c   = 0;
287:               valA[5]    = 1.0 / (hz * hz);
288:               col[6].i   = ex - 1;
289:               col[6].j   = ey;
290:               col[6].k   = ez;
291:               col[6].loc = ELEMENT;
292:               col[6].c   = 0;
293:               valA[6]    = 1.0 / hx;
294:               col[7].i   = ex;
295:               col[7].j   = ey;
296:               col[7].k   = ez;
297:               col[7].loc = ELEMENT;
298:               col[7].c   = 0;
299:               valA[7]    = -1.0 / hx;
300:             }
301:           } else if (ey == N[1] - 1) {
302:             if (ez == 0) {
303:               nEntries   = 7;
304:               col[0].i   = ex;
305:               col[0].j   = ey;
306:               col[0].k   = ez;
307:               col[0].loc = LEFT;
308:               col[0].c   = 0;
309:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
310:               col[1].i   = ex;
311:               col[1].j   = ey - 1;
312:               col[1].k   = ez;
313:               col[1].loc = LEFT;
314:               col[1].c   = 0;
315:               valA[1]    = 1.0 / (hy * hy);
316:               /* Missing up term */
317:               col[2].i   = ex - 1;
318:               col[2].j   = ey;
319:               col[2].k   = ez;
320:               col[2].loc = LEFT;
321:               col[2].c   = 0;
322:               valA[2]    = 1.0 / (hx * hx);
323:               col[3].i   = ex + 1;
324:               col[3].j   = ey;
325:               col[3].k   = ez;
326:               col[3].loc = LEFT;
327:               col[3].c   = 0;
328:               valA[3]    = 1.0 / (hx * hx);
329:               /* Missing back entry */
330:               col[4].i   = ex;
331:               col[4].j   = ey;
332:               col[4].k   = ez + 1;
333:               col[4].loc = LEFT;
334:               col[4].c   = 0;
335:               valA[4]    = 1.0 / (hz * hz);
336:               col[5].i   = ex - 1;
337:               col[5].j   = ey;
338:               col[5].k   = ez;
339:               col[5].loc = ELEMENT;
340:               col[5].c   = 0;
341:               valA[5]    = 1.0 / hx;
342:               col[6].i   = ex;
343:               col[6].j   = ey;
344:               col[6].k   = ez;
345:               col[6].loc = ELEMENT;
346:               col[6].c   = 0;
347:               valA[6]    = -1.0 / hx;
348:             } else if (ez == N[2] - 1) {
349:               nEntries   = 7;
350:               col[0].i   = ex;
351:               col[0].j   = ey;
352:               col[0].k   = ez;
353:               col[0].loc = LEFT;
354:               col[0].c   = 0;
355:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
356:               col[1].i   = ex;
357:               col[1].j   = ey - 1;
358:               col[1].k   = ez;
359:               col[1].loc = LEFT;
360:               col[1].c   = 0;
361:               valA[1]    = 1.0 / (hy * hy);
362:               /* Missing up term */
363:               col[2].i   = ex - 1;
364:               col[2].j   = ey;
365:               col[2].k   = ez;
366:               col[2].loc = LEFT;
367:               col[2].c   = 0;
368:               valA[2]    = 1.0 / (hx * hx);
369:               col[3].i   = ex + 1;
370:               col[3].j   = ey;
371:               col[3].k   = ez;
372:               col[3].loc = LEFT;
373:               col[3].c   = 0;
374:               valA[3]    = 1.0 / (hx * hx);
375:               col[4].i   = ex;
376:               col[4].j   = ey;
377:               col[4].k   = ez - 1;
378:               col[4].loc = LEFT;
379:               col[4].c   = 0;
380:               valA[4]    = 1.0 / (hz * hz);
381:               /* Missing front term */
382:               col[5].i   = ex - 1;
383:               col[5].j   = ey;
384:               col[5].k   = ez;
385:               col[5].loc = ELEMENT;
386:               col[5].c   = 0;
387:               valA[5]    = 1.0 / hx;
388:               col[6].i   = ex;
389:               col[6].j   = ey;
390:               col[6].k   = ez;
391:               col[6].loc = ELEMENT;
392:               col[6].c   = 0;
393:               valA[6]    = -1.0 / hx;
394:             } else {
395:               nEntries   = 8;
396:               col[0].i   = ex;
397:               col[0].j   = ey;
398:               col[0].k   = ez;
399:               col[0].loc = LEFT;
400:               col[0].c   = 0;
401:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
402:               col[1].i   = ex;
403:               col[1].j   = ey - 1;
404:               col[1].k   = ez;
405:               col[1].loc = LEFT;
406:               col[1].c   = 0;
407:               valA[1]    = 1.0 / (hy * hy);
408:               /* Missing up term */
409:               col[2].i   = ex - 1;
410:               col[2].j   = ey;
411:               col[2].k   = ez;
412:               col[2].loc = LEFT;
413:               col[2].c   = 0;
414:               valA[2]    = 1.0 / (hx * hx);
415:               col[3].i   = ex + 1;
416:               col[3].j   = ey;
417:               col[3].k   = ez;
418:               col[3].loc = LEFT;
419:               col[3].c   = 0;
420:               valA[3]    = 1.0 / (hx * hx);
421:               col[4].i   = ex;
422:               col[4].j   = ey;
423:               col[4].k   = ez - 1;
424:               col[4].loc = LEFT;
425:               col[4].c   = 0;
426:               valA[4]    = 1.0 / (hz * hz);
427:               col[5].i   = ex;
428:               col[5].j   = ey;
429:               col[5].k   = ez + 1;
430:               col[5].loc = LEFT;
431:               col[5].c   = 0;
432:               valA[5]    = 1.0 / (hz * hz);
433:               col[6].i   = ex - 1;
434:               col[6].j   = ey;
435:               col[6].k   = ez;
436:               col[6].loc = ELEMENT;
437:               col[6].c   = 0;
438:               valA[6]    = 1.0 / hx;
439:               col[7].i   = ex;
440:               col[7].j   = ey;
441:               col[7].k   = ez;
442:               col[7].loc = ELEMENT;
443:               col[7].c   = 0;
444:               valA[7]    = -1.0 / hx;
445:             }
446:           } else if (ez == 0) {
447:             nEntries   = 8;
448:             col[0].i   = ex;
449:             col[0].j   = ey;
450:             col[0].k   = ez;
451:             col[0].loc = LEFT;
452:             col[0].c   = 0;
453:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
454:             col[1].i   = ex;
455:             col[1].j   = ey - 1;
456:             col[1].k   = ez;
457:             col[1].loc = LEFT;
458:             col[1].c   = 0;
459:             valA[1]    = 1.0 / (hy * hy);
460:             col[2].i   = ex;
461:             col[2].j   = ey + 1;
462:             col[2].k   = ez;
463:             col[2].loc = LEFT;
464:             col[2].c   = 0;
465:             valA[2]    = 1.0 / (hy * hy);
466:             col[3].i   = ex - 1;
467:             col[3].j   = ey;
468:             col[3].k   = ez;
469:             col[3].loc = LEFT;
470:             col[3].c   = 0;
471:             valA[3]    = 1.0 / (hx * hx);
472:             col[4].i   = ex + 1;
473:             col[4].j   = ey;
474:             col[4].k   = ez;
475:             col[4].loc = LEFT;
476:             col[4].c   = 0;
477:             valA[4]    = 1.0 / (hx * hx);
478:             /* Missing back term */
479:             col[5].i   = ex;
480:             col[5].j   = ey;
481:             col[5].k   = ez + 1;
482:             col[5].loc = LEFT;
483:             col[5].c   = 0;
484:             valA[5]    = 1.0 / (hz * hz);
485:             col[6].i   = ex - 1;
486:             col[6].j   = ey;
487:             col[6].k   = ez;
488:             col[6].loc = ELEMENT;
489:             col[6].c   = 0;
490:             valA[6]    = 1.0 / hx;
491:             col[7].i   = ex;
492:             col[7].j   = ey;
493:             col[7].k   = ez;
494:             col[7].loc = ELEMENT;
495:             col[7].c   = 0;
496:             valA[7]    = -1.0 / hx;
497:           } else if (ez == N[2] - 1) {
498:             nEntries   = 8;
499:             col[0].i   = ex;
500:             col[0].j   = ey;
501:             col[0].k   = ez;
502:             col[0].loc = LEFT;
503:             col[0].c   = 0;
504:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
505:             col[1].i   = ex;
506:             col[1].j   = ey - 1;
507:             col[1].k   = ez;
508:             col[1].loc = LEFT;
509:             col[1].c   = 0;
510:             valA[1]    = 1.0 / (hy * hy);
511:             col[2].i   = ex;
512:             col[2].j   = ey + 1;
513:             col[2].k   = ez;
514:             col[2].loc = LEFT;
515:             col[2].c   = 0;
516:             valA[2]    = 1.0 / (hy * hy);
517:             col[3].i   = ex - 1;
518:             col[3].j   = ey;
519:             col[3].k   = ez;
520:             col[3].loc = LEFT;
521:             col[3].c   = 0;
522:             valA[3]    = 1.0 / (hx * hx);
523:             col[4].i   = ex + 1;
524:             col[4].j   = ey;
525:             col[4].k   = ez;
526:             col[4].loc = LEFT;
527:             col[4].c   = 0;
528:             valA[4]    = 1.0 / (hx * hx);
529:             col[5].i   = ex;
530:             col[5].j   = ey;
531:             col[5].k   = ez - 1;
532:             col[5].loc = LEFT;
533:             col[5].c   = 0;
534:             valA[5]    = 1.0 / (hz * hz);
535:             /* Missing front term */
536:             col[6].i   = ex - 1;
537:             col[6].j   = ey;
538:             col[6].k   = ez;
539:             col[6].loc = ELEMENT;
540:             col[6].c   = 0;
541:             valA[6]    = 1.0 / hx;
542:             col[7].i   = ex;
543:             col[7].j   = ey;
544:             col[7].k   = ez;
545:             col[7].loc = ELEMENT;
546:             col[7].c   = 0;
547:             valA[7]    = -1.0 / hx;
548:           } else {
549:             nEntries   = 9;
550:             col[0].i   = ex;
551:             col[0].j   = ey;
552:             col[0].k   = ez;
553:             col[0].loc = LEFT;
554:             col[0].c   = 0;
555:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
556:             col[1].i   = ex;
557:             col[1].j   = ey - 1;
558:             col[1].k   = ez;
559:             col[1].loc = LEFT;
560:             col[1].c   = 0;
561:             valA[1]    = 1.0 / (hy * hy);
562:             col[2].i   = ex;
563:             col[2].j   = ey + 1;
564:             col[2].k   = ez;
565:             col[2].loc = LEFT;
566:             col[2].c   = 0;
567:             valA[2]    = 1.0 / (hy * hy);
568:             col[3].i   = ex - 1;
569:             col[3].j   = ey;
570:             col[3].k   = ez;
571:             col[3].loc = LEFT;
572:             col[3].c   = 0;
573:             valA[3]    = 1.0 / (hx * hx);
574:             col[4].i   = ex + 1;
575:             col[4].j   = ey;
576:             col[4].k   = ez;
577:             col[4].loc = LEFT;
578:             col[4].c   = 0;
579:             valA[4]    = 1.0 / (hx * hx);
580:             col[5].i   = ex;
581:             col[5].j   = ey;
582:             col[5].k   = ez - 1;
583:             col[5].loc = LEFT;
584:             col[5].c   = 0;
585:             valA[5]    = 1.0 / (hz * hz);
586:             col[6].i   = ex;
587:             col[6].j   = ey;
588:             col[6].k   = ez + 1;
589:             col[6].loc = LEFT;
590:             col[6].c   = 0;
591:             valA[6]    = 1.0 / (hz * hz);
592:             col[7].i   = ex - 1;
593:             col[7].j   = ey;
594:             col[7].k   = ez;
595:             col[7].loc = ELEMENT;
596:             col[7].c   = 0;
597:             valA[7]    = 1.0 / hx;
598:             col[8].i   = ex;
599:             col[8].j   = ey;
600:             col[8].k   = ez;
601:             col[8].loc = ELEMENT;
602:             col[8].c   = 0;
603:             valA[8]    = -1.0 / hx;
604:           }
605:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
606:         }

608:         /* Equation on bottom face of this element */
609:         if (ey == 0) {
610:           /* Bottom boundary velocity Dirichlet */
611:           DMStagStencil     row;
612:           const PetscScalar valA = 1.0;
613:           row.i                  = ex;
614:           row.j                  = ey;
615:           row.k                  = ez;
616:           row.loc                = DOWN;
617:           row.c                  = 0;
618:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
619:         } else {
620:           /* Y-momentum equation, (v_xx + v_yy + v_zz) - p_y = f^y */
621:           DMStagStencil row, col[9];
622:           PetscScalar   valA[9];
623:           PetscInt      nEntries;

625:           row.i   = ex;
626:           row.j   = ey;
627:           row.k   = ez;
628:           row.loc = DOWN;
629:           row.c   = 0;
630:           if (ex == 0) {
631:             if (ez == 0) {
632:               nEntries   = 7;
633:               col[0].i   = ex;
634:               col[0].j   = ey;
635:               col[0].k   = ez;
636:               col[0].loc = DOWN;
637:               col[0].c   = 0;
638:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
639:               col[1].i   = ex;
640:               col[1].j   = ey - 1;
641:               col[1].k   = ez;
642:               col[1].loc = DOWN;
643:               col[1].c   = 0;
644:               valA[1]    = 1.0 / (hy * hy);
645:               col[2].i   = ex;
646:               col[2].j   = ey + 1;
647:               col[2].k   = ez;
648:               col[2].loc = DOWN;
649:               col[2].c   = 0;
650:               valA[2]    = 1.0 / (hy * hy);
651:               /* Left term missing */
652:               col[3].i   = ex + 1;
653:               col[3].j   = ey;
654:               col[3].k   = ez;
655:               col[3].loc = DOWN;
656:               col[3].c   = 0;
657:               valA[3]    = 1.0 / (hx * hx);
658:               /* Back term missing */
659:               col[4].i   = ex;
660:               col[4].j   = ey;
661:               col[4].k   = ez + 1;
662:               col[4].loc = DOWN;
663:               col[4].c   = 0;
664:               valA[4]    = 1.0 / (hz * hz);
665:               col[5].i   = ex;
666:               col[5].j   = ey - 1;
667:               col[5].k   = ez;
668:               col[5].loc = ELEMENT;
669:               col[5].c   = 0;
670:               valA[5]    = 1.0 / hy;
671:               col[6].i   = ex;
672:               col[6].j   = ey;
673:               col[6].k   = ez;
674:               col[6].loc = ELEMENT;
675:               col[6].c   = 0;
676:               valA[6]    = -1.0 / hy;
677:             } else if (ez == N[2] - 1) {
678:               nEntries   = 7;
679:               col[0].i   = ex;
680:               col[0].j   = ey;
681:               col[0].k   = ez;
682:               col[0].loc = DOWN;
683:               col[0].c   = 0;
684:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
685:               col[1].i   = ex;
686:               col[1].j   = ey - 1;
687:               col[1].k   = ez;
688:               col[1].loc = DOWN;
689:               col[1].c   = 0;
690:               valA[1]    = 1.0 / (hy * hy);
691:               col[2].i   = ex;
692:               col[2].j   = ey + 1;
693:               col[2].k   = ez;
694:               col[2].loc = DOWN;
695:               col[2].c   = 0;
696:               valA[2]    = 1.0 / (hy * hy);
697:               /* Left term missing */
698:               col[3].i   = ex + 1;
699:               col[3].j   = ey;
700:               col[3].k   = ez;
701:               col[3].loc = DOWN;
702:               col[3].c   = 0;
703:               valA[3]    = 1.0 / (hx * hx);
704:               col[4].i   = ex;
705:               col[4].j   = ey;
706:               col[4].k   = ez - 1;
707:               col[4].loc = DOWN;
708:               col[4].c   = 0;
709:               valA[4]    = 1.0 / (hz * hz);
710:               /* Front term missing */
711:               col[5].i   = ex;
712:               col[5].j   = ey - 1;
713:               col[5].k   = ez;
714:               col[5].loc = ELEMENT;
715:               col[5].c   = 0;
716:               valA[5]    = 1.0 / hy;
717:               col[6].i   = ex;
718:               col[6].j   = ey;
719:               col[6].k   = ez;
720:               col[6].loc = ELEMENT;
721:               col[6].c   = 0;
722:               valA[6]    = -1.0 / hy;
723:             } else {
724:               nEntries   = 8;
725:               col[0].i   = ex;
726:               col[0].j   = ey;
727:               col[0].k   = ez;
728:               col[0].loc = DOWN;
729:               col[0].c   = 0;
730:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
731:               col[1].i   = ex;
732:               col[1].j   = ey - 1;
733:               col[1].k   = ez;
734:               col[1].loc = DOWN;
735:               col[1].c   = 0;
736:               valA[1]    = 1.0 / (hy * hy);
737:               col[2].i   = ex;
738:               col[2].j   = ey + 1;
739:               col[2].k   = ez;
740:               col[2].loc = DOWN;
741:               col[2].c   = 0;
742:               valA[2]    = 1.0 / (hy * hy);
743:               /* Left term missing */
744:               col[3].i   = ex + 1;
745:               col[3].j   = ey;
746:               col[3].k   = ez;
747:               col[3].loc = DOWN;
748:               col[3].c   = 0;
749:               valA[3]    = 1.0 / (hx * hx);
750:               col[4].i   = ex;
751:               col[4].j   = ey;
752:               col[4].k   = ez - 1;
753:               col[4].loc = DOWN;
754:               col[4].c   = 0;
755:               valA[4]    = 1.0 / (hz * hz);
756:               col[5].i   = ex;
757:               col[5].j   = ey;
758:               col[5].k   = ez + 1;
759:               col[5].loc = DOWN;
760:               col[5].c   = 0;
761:               valA[5]    = 1.0 / (hz * hz);
762:               col[6].i   = ex;
763:               col[6].j   = ey - 1;
764:               col[6].k   = ez;
765:               col[6].loc = ELEMENT;
766:               col[6].c   = 0;
767:               valA[6]    = 1.0 / hy;
768:               col[7].i   = ex;
769:               col[7].j   = ey;
770:               col[7].k   = ez;
771:               col[7].loc = ELEMENT;
772:               col[7].c   = 0;
773:               valA[7]    = -1.0 / hy;
774:             }
775:           } else if (ex == N[0] - 1) {
776:             if (ez == 0) {
777:               nEntries   = 7;
778:               col[0].i   = ex;
779:               col[0].j   = ey;
780:               col[0].k   = ez;
781:               col[0].loc = DOWN;
782:               col[0].c   = 0;
783:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
784:               col[1].i   = ex;
785:               col[1].j   = ey - 1;
786:               col[1].k   = ez;
787:               col[1].loc = DOWN;
788:               col[1].c   = 0;
789:               valA[1]    = 1.0 / (hy * hy);
790:               col[2].i   = ex;
791:               col[2].j   = ey + 1;
792:               col[2].k   = ez;
793:               col[2].loc = DOWN;
794:               col[2].c   = 0;
795:               valA[2]    = 1.0 / (hy * hy);
796:               col[3].i   = ex - 1;
797:               col[3].j   = ey;
798:               col[3].k   = ez;
799:               col[3].loc = DOWN;
800:               col[3].c   = 0;
801:               valA[3]    = 1.0 / (hx * hx);
802:               /* Right term missing */
803:               /* Back term missing */
804:               col[4].i   = ex;
805:               col[4].j   = ey;
806:               col[4].k   = ez + 1;
807:               col[4].loc = DOWN;
808:               col[4].c   = 0;
809:               valA[4]    = 1.0 / (hz * hz);
810:               col[5].i   = ex;
811:               col[5].j   = ey - 1;
812:               col[5].k   = ez;
813:               col[5].loc = ELEMENT;
814:               col[5].c   = 0;
815:               valA[5]    = 1.0 / hy;
816:               col[6].i   = ex;
817:               col[6].j   = ey;
818:               col[6].k   = ez;
819:               col[6].loc = ELEMENT;
820:               col[6].c   = 0;
821:               valA[6]    = -1.0 / hy;
822:             } else if (ez == N[2] - 1) {
823:               nEntries   = 7;
824:               col[0].i   = ex;
825:               col[0].j   = ey;
826:               col[0].k   = ez;
827:               col[0].loc = DOWN;
828:               col[0].c   = 0;
829:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
830:               col[1].i   = ex;
831:               col[1].j   = ey - 1;
832:               col[1].k   = ez;
833:               col[1].loc = DOWN;
834:               col[1].c   = 0;
835:               valA[1]    = 1.0 / (hy * hy);
836:               col[2].i   = ex;
837:               col[2].j   = ey + 1;
838:               col[2].k   = ez;
839:               col[2].loc = DOWN;
840:               col[2].c   = 0;
841:               valA[2]    = 1.0 / (hy * hy);
842:               col[3].i   = ex - 1;
843:               col[3].j   = ey;
844:               col[3].k   = ez;
845:               col[3].loc = DOWN;
846:               col[3].c   = 0;
847:               valA[3]    = 1.0 / (hx * hx);
848:               /* Right term missing */
849:               col[4].i   = ex;
850:               col[4].j   = ey;
851:               col[4].k   = ez - 1;
852:               col[4].loc = DOWN;
853:               col[4].c   = 0;
854:               valA[4]    = 1.0 / (hz * hz);
855:               /* Front term missing */
856:               col[5].i   = ex;
857:               col[5].j   = ey - 1;
858:               col[5].k   = ez;
859:               col[5].loc = ELEMENT;
860:               col[5].c   = 0;
861:               valA[5]    = 1.0 / hy;
862:               col[6].i   = ex;
863:               col[6].j   = ey;
864:               col[6].k   = ez;
865:               col[6].loc = ELEMENT;
866:               col[6].c   = 0;
867:               valA[6]    = -1.0 / hy;
868:             } else {
869:               nEntries   = 8;
870:               col[0].i   = ex;
871:               col[0].j   = ey;
872:               col[0].k   = ez;
873:               col[0].loc = DOWN;
874:               col[0].c   = 0;
875:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
876:               col[1].i   = ex;
877:               col[1].j   = ey - 1;
878:               col[1].k   = ez;
879:               col[1].loc = DOWN;
880:               col[1].c   = 0;
881:               valA[1]    = 1.0 / (hy * hy);
882:               col[2].i   = ex;
883:               col[2].j   = ey + 1;
884:               col[2].k   = ez;
885:               col[2].loc = DOWN;
886:               col[2].c   = 0;
887:               valA[2]    = 1.0 / (hy * hy);
888:               col[3].i   = ex - 1;
889:               col[3].j   = ey;
890:               col[3].k   = ez;
891:               col[3].loc = DOWN;
892:               col[3].c   = 0;
893:               valA[3]    = 1.0 / (hx * hx);
894:               /* Right term missing */
895:               col[4].i   = ex;
896:               col[4].j   = ey;
897:               col[4].k   = ez - 1;
898:               col[4].loc = DOWN;
899:               col[4].c   = 0;
900:               valA[4]    = 1.0 / (hz * hz);
901:               col[5].i   = ex;
902:               col[5].j   = ey;
903:               col[5].k   = ez + 1;
904:               col[5].loc = DOWN;
905:               col[5].c   = 0;
906:               valA[5]    = 1.0 / (hz * hz);
907:               col[6].i   = ex;
908:               col[6].j   = ey - 1;
909:               col[6].k   = ez;
910:               col[6].loc = ELEMENT;
911:               col[6].c   = 0;
912:               valA[6]    = 1.0 / hy;
913:               col[7].i   = ex;
914:               col[7].j   = ey;
915:               col[7].k   = ez;
916:               col[7].loc = ELEMENT;
917:               col[7].c   = 0;
918:               valA[7]    = -1.0 / hy;
919:             }
920:           } else if (ez == 0) {
921:             nEntries   = 8;
922:             col[0].i   = ex;
923:             col[0].j   = ey;
924:             col[0].k   = ez;
925:             col[0].loc = DOWN;
926:             col[0].c   = 0;
927:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
928:             col[1].i   = ex;
929:             col[1].j   = ey - 1;
930:             col[1].k   = ez;
931:             col[1].loc = DOWN;
932:             col[1].c   = 0;
933:             valA[1]    = 1.0 / (hy * hy);
934:             col[2].i   = ex;
935:             col[2].j   = ey + 1;
936:             col[2].k   = ez;
937:             col[2].loc = DOWN;
938:             col[2].c   = 0;
939:             valA[2]    = 1.0 / (hy * hy);
940:             col[3].i   = ex - 1;
941:             col[3].j   = ey;
942:             col[3].k   = ez;
943:             col[3].loc = DOWN;
944:             col[3].c   = 0;
945:             valA[3]    = 1.0 / (hx * hx);
946:             col[4].i   = ex + 1;
947:             col[4].j   = ey;
948:             col[4].k   = ez;
949:             col[4].loc = DOWN;
950:             col[4].c   = 0;
951:             valA[4]    = 1.0 / (hx * hx);
952:             /* Back term missing */
953:             col[5].i   = ex;
954:             col[5].j   = ey;
955:             col[5].k   = ez + 1;
956:             col[5].loc = DOWN;
957:             col[5].c   = 0;
958:             valA[5]    = 1.0 / (hz * hz);
959:             col[6].i   = ex;
960:             col[6].j   = ey - 1;
961:             col[6].k   = ez;
962:             col[6].loc = ELEMENT;
963:             col[6].c   = 0;
964:             valA[6]    = 1.0 / hy;
965:             col[7].i   = ex;
966:             col[7].j   = ey;
967:             col[7].k   = ez;
968:             col[7].loc = ELEMENT;
969:             col[7].c   = 0;
970:             valA[7]    = -1.0 / hy;
971:           } else if (ez == N[2] - 1) {
972:             nEntries   = 8;
973:             col[0].i   = ex;
974:             col[0].j   = ey;
975:             col[0].k   = ez;
976:             col[0].loc = DOWN;
977:             col[0].c   = 0;
978:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
979:             col[1].i   = ex;
980:             col[1].j   = ey - 1;
981:             col[1].k   = ez;
982:             col[1].loc = DOWN;
983:             col[1].c   = 0;
984:             valA[1]    = 1.0 / (hy * hy);
985:             col[2].i   = ex;
986:             col[2].j   = ey + 1;
987:             col[2].k   = ez;
988:             col[2].loc = DOWN;
989:             col[2].c   = 0;
990:             valA[2]    = 1.0 / (hy * hy);
991:             col[3].i   = ex - 1;
992:             col[3].j   = ey;
993:             col[3].k   = ez;
994:             col[3].loc = DOWN;
995:             col[3].c   = 0;
996:             valA[3]    = 1.0 / (hx * hx);
997:             col[4].i   = ex + 1;
998:             col[4].j   = ey;
999:             col[4].k   = ez;
1000:             col[4].loc = DOWN;
1001:             col[4].c   = 0;
1002:             valA[4]    = 1.0 / (hx * hx);
1003:             col[5].i   = ex;
1004:             col[5].j   = ey;
1005:             col[5].k   = ez - 1;
1006:             col[5].loc = DOWN;
1007:             col[5].c   = 0;
1008:             valA[5]    = 1.0 / (hz * hz);
1009:             /* Front term missing */
1010:             col[6].i   = ex;
1011:             col[6].j   = ey - 1;
1012:             col[6].k   = ez;
1013:             col[6].loc = ELEMENT;
1014:             col[6].c   = 0;
1015:             valA[6]    = 1.0 / hy;
1016:             col[7].i   = ex;
1017:             col[7].j   = ey;
1018:             col[7].k   = ez;
1019:             col[7].loc = ELEMENT;
1020:             col[7].c   = 0;
1021:             valA[7]    = -1.0 / hy;
1022:           } else {
1023:             nEntries   = 9;
1024:             col[0].i   = ex;
1025:             col[0].j   = ey;
1026:             col[0].k   = ez;
1027:             col[0].loc = DOWN;
1028:             col[0].c   = 0;
1029:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1030:             col[1].i   = ex;
1031:             col[1].j   = ey - 1;
1032:             col[1].k   = ez;
1033:             col[1].loc = DOWN;
1034:             col[1].c   = 0;
1035:             valA[1]    = 1.0 / (hy * hy);
1036:             col[2].i   = ex;
1037:             col[2].j   = ey + 1;
1038:             col[2].k   = ez;
1039:             col[2].loc = DOWN;
1040:             col[2].c   = 0;
1041:             valA[2]    = 1.0 / (hy * hy);
1042:             col[3].i   = ex - 1;
1043:             col[3].j   = ey;
1044:             col[3].k   = ez;
1045:             col[3].loc = DOWN;
1046:             col[3].c   = 0;
1047:             valA[3]    = 1.0 / (hx * hx);
1048:             col[4].i   = ex + 1;
1049:             col[4].j   = ey;
1050:             col[4].k   = ez;
1051:             col[4].loc = DOWN;
1052:             col[4].c   = 0;
1053:             valA[4]    = 1.0 / (hx * hx);
1054:             col[5].i   = ex;
1055:             col[5].j   = ey;
1056:             col[5].k   = ez - 1;
1057:             col[5].loc = DOWN;
1058:             col[5].c   = 0;
1059:             valA[5]    = 1.0 / (hz * hz);
1060:             col[6].i   = ex;
1061:             col[6].j   = ey;
1062:             col[6].k   = ez + 1;
1063:             col[6].loc = DOWN;
1064:             col[6].c   = 0;
1065:             valA[6]    = 1.0 / (hz * hz);
1066:             col[7].i   = ex;
1067:             col[7].j   = ey - 1;
1068:             col[7].k   = ez;
1069:             col[7].loc = ELEMENT;
1070:             col[7].c   = 0;
1071:             valA[7]    = 1.0 / hy;
1072:             col[8].i   = ex;
1073:             col[8].j   = ey;
1074:             col[8].k   = ez;
1075:             col[8].loc = ELEMENT;
1076:             col[8].c   = 0;
1077:             valA[8]    = -1.0 / hy;
1078:           }
1079:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
1080:         }

1082:         /* Equation on back face of this element */
1083:         if (ez == 0) {
1084:           /* Back boundary velocity Dirichlet */
1085:           DMStagStencil     row;
1086:           const PetscScalar valA = 1.0;
1087:           row.i                  = ex;
1088:           row.j                  = ey;
1089:           row.k                  = ez;
1090:           row.loc                = BACK;
1091:           row.c                  = 0;
1092:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
1093:         } else {
1094:           /* Z-momentum equation, (w_xx + w_yy + w_zz) - p_z = f^z */
1095:           DMStagStencil row, col[9];
1096:           PetscScalar   valA[9];
1097:           PetscInt      nEntries;

1099:           row.i   = ex;
1100:           row.j   = ey;
1101:           row.k   = ez;
1102:           row.loc = BACK;
1103:           row.c   = 0;
1104:           if (ex == 0) {
1105:             if (ey == 0) {
1106:               nEntries   = 7;
1107:               col[0].i   = ex;
1108:               col[0].j   = ey;
1109:               col[0].k   = ez;
1110:               col[0].loc = BACK;
1111:               col[0].c   = 0;
1112:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1113:               /* Down term missing */
1114:               col[1].i   = ex;
1115:               col[1].j   = ey + 1;
1116:               col[1].k   = ez;
1117:               col[1].loc = BACK;
1118:               col[1].c   = 0;
1119:               valA[1]    = 1.0 / (hy * hy);
1120:               /* Left term missing */
1121:               col[2].i   = ex + 1;
1122:               col[2].j   = ey;
1123:               col[2].k   = ez;
1124:               col[2].loc = BACK;
1125:               col[2].c   = 0;
1126:               valA[2]    = 1.0 / (hx * hx);
1127:               col[3].i   = ex;
1128:               col[3].j   = ey;
1129:               col[3].k   = ez - 1;
1130:               col[3].loc = BACK;
1131:               col[3].c   = 0;
1132:               valA[3]    = 1.0 / (hz * hz);
1133:               col[4].i   = ex;
1134:               col[4].j   = ey;
1135:               col[4].k   = ez + 1;
1136:               col[4].loc = BACK;
1137:               col[4].c   = 0;
1138:               valA[4]    = 1.0 / (hz * hz);
1139:               col[5].i   = ex;
1140:               col[5].j   = ey;
1141:               col[5].k   = ez - 1;
1142:               col[5].loc = ELEMENT;
1143:               col[5].c   = 0;
1144:               valA[5]    = 1.0 / hz;
1145:               col[6].i   = ex;
1146:               col[6].j   = ey;
1147:               col[6].k   = ez;
1148:               col[6].loc = ELEMENT;
1149:               col[6].c   = 0;
1150:               valA[6]    = -1.0 / hz;
1151:             } else if (ey == N[1] - 1) {
1152:               nEntries   = 7;
1153:               col[0].i   = ex;
1154:               col[0].j   = ey;
1155:               col[0].k   = ez;
1156:               col[0].loc = BACK;
1157:               col[0].c   = 0;
1158:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1159:               col[1].i   = ex;
1160:               col[1].j   = ey - 1;
1161:               col[1].k   = ez;
1162:               col[1].loc = BACK;
1163:               col[1].c   = 0;
1164:               valA[1]    = 1.0 / (hy * hy);
1165:               /* Up term missing */
1166:               /* Left term missing */
1167:               col[2].i   = ex + 1;
1168:               col[2].j   = ey;
1169:               col[2].k   = ez;
1170:               col[2].loc = BACK;
1171:               col[2].c   = 0;
1172:               valA[2]    = 1.0 / (hx * hx);
1173:               col[3].i   = ex;
1174:               col[3].j   = ey;
1175:               col[3].k   = ez - 1;
1176:               col[3].loc = BACK;
1177:               col[3].c   = 0;
1178:               valA[3]    = 1.0 / (hz * hz);
1179:               col[4].i   = ex;
1180:               col[4].j   = ey;
1181:               col[4].k   = ez + 1;
1182:               col[4].loc = BACK;
1183:               col[4].c   = 0;
1184:               valA[4]    = 1.0 / (hz * hz);
1185:               col[5].i   = ex;
1186:               col[5].j   = ey;
1187:               col[5].k   = ez - 1;
1188:               col[5].loc = ELEMENT;
1189:               col[5].c   = 0;
1190:               valA[5]    = 1.0 / hz;
1191:               col[6].i   = ex;
1192:               col[6].j   = ey;
1193:               col[6].k   = ez;
1194:               col[6].loc = ELEMENT;
1195:               col[6].c   = 0;
1196:               valA[6]    = -1.0 / hz;
1197:             } else {
1198:               nEntries   = 8;
1199:               col[0].i   = ex;
1200:               col[0].j   = ey;
1201:               col[0].k   = ez;
1202:               col[0].loc = BACK;
1203:               col[0].c   = 0;
1204:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1205:               col[1].i   = ex;
1206:               col[1].j   = ey - 1;
1207:               col[1].k   = ez;
1208:               col[1].loc = BACK;
1209:               col[1].c   = 0;
1210:               valA[1]    = 1.0 / (hy * hy);
1211:               col[2].i   = ex;
1212:               col[2].j   = ey + 1;
1213:               col[2].k   = ez;
1214:               col[2].loc = BACK;
1215:               col[2].c   = 0;
1216:               valA[2]    = 1.0 / (hy * hy);
1217:               /* Left term missing */
1218:               col[3].i   = ex + 1;
1219:               col[3].j   = ey;
1220:               col[3].k   = ez;
1221:               col[3].loc = BACK;
1222:               col[3].c   = 0;
1223:               valA[3]    = 1.0 / (hx * hx);
1224:               col[4].i   = ex;
1225:               col[4].j   = ey;
1226:               col[4].k   = ez - 1;
1227:               col[4].loc = BACK;
1228:               col[4].c   = 0;
1229:               valA[4]    = 1.0 / (hz * hz);
1230:               col[5].i   = ex;
1231:               col[5].j   = ey;
1232:               col[5].k   = ez + 1;
1233:               col[5].loc = BACK;
1234:               col[5].c   = 0;
1235:               valA[5]    = 1.0 / (hz * hz);
1236:               col[6].i   = ex;
1237:               col[6].j   = ey;
1238:               col[6].k   = ez - 1;
1239:               col[6].loc = ELEMENT;
1240:               col[6].c   = 0;
1241:               valA[6]    = 1.0 / hz;
1242:               col[7].i   = ex;
1243:               col[7].j   = ey;
1244:               col[7].k   = ez;
1245:               col[7].loc = ELEMENT;
1246:               col[7].c   = 0;
1247:               valA[7]    = -1.0 / hz;
1248:             }
1249:           } else if (ex == N[0] - 1) {
1250:             if (ey == 0) {
1251:               nEntries   = 7;
1252:               col[0].i   = ex;
1253:               col[0].j   = ey;
1254:               col[0].k   = ez;
1255:               col[0].loc = BACK;
1256:               col[0].c   = 0;
1257:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1258:               /* Down term missing */
1259:               col[1].i   = ex;
1260:               col[1].j   = ey + 1;
1261:               col[1].k   = ez;
1262:               col[1].loc = BACK;
1263:               col[1].c   = 0;
1264:               valA[1]    = 1.0 / (hy * hy);
1265:               col[2].i   = ex - 1;
1266:               col[2].j   = ey;
1267:               col[2].k   = ez;
1268:               col[2].loc = BACK;
1269:               col[2].c   = 0;
1270:               valA[2]    = 1.0 / (hx * hx);
1271:               /* Right term missing */
1272:               col[3].i   = ex;
1273:               col[3].j   = ey;
1274:               col[3].k   = ez - 1;
1275:               col[3].loc = BACK;
1276:               col[3].c   = 0;
1277:               valA[3]    = 1.0 / (hz * hz);
1278:               col[4].i   = ex;
1279:               col[4].j   = ey;
1280:               col[4].k   = ez + 1;
1281:               col[4].loc = BACK;
1282:               col[4].c   = 0;
1283:               valA[4]    = 1.0 / (hz * hz);
1284:               col[5].i   = ex;
1285:               col[5].j   = ey;
1286:               col[5].k   = ez - 1;
1287:               col[5].loc = ELEMENT;
1288:               col[5].c   = 0;
1289:               valA[5]    = 1.0 / hz;
1290:               col[6].i   = ex;
1291:               col[6].j   = ey;
1292:               col[6].k   = ez;
1293:               col[6].loc = ELEMENT;
1294:               col[6].c   = 0;
1295:               valA[6]    = -1.0 / hz;
1296:             } else if (ey == N[1] - 1) {
1297:               nEntries   = 7;
1298:               col[0].i   = ex;
1299:               col[0].j   = ey;
1300:               col[0].k   = ez;
1301:               col[0].loc = BACK;
1302:               col[0].c   = 0;
1303:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1304:               col[1].i   = ex;
1305:               col[1].j   = ey - 1;
1306:               col[1].k   = ez;
1307:               col[1].loc = BACK;
1308:               col[1].c   = 0;
1309:               valA[1]    = 1.0 / (hy * hy);
1310:               /* Up term missing */
1311:               col[2].i   = ex - 1;
1312:               col[2].j   = ey;
1313:               col[2].k   = ez;
1314:               col[2].loc = BACK;
1315:               col[2].c   = 0;
1316:               valA[2]    = 1.0 / (hx * hx);
1317:               /* Right term missing */
1318:               col[3].i   = ex;
1319:               col[3].j   = ey;
1320:               col[3].k   = ez - 1;
1321:               col[3].loc = BACK;
1322:               col[3].c   = 0;
1323:               valA[3]    = 1.0 / (hz * hz);
1324:               col[4].i   = ex;
1325:               col[4].j   = ey;
1326:               col[4].k   = ez + 1;
1327:               col[4].loc = BACK;
1328:               col[4].c   = 0;
1329:               valA[4]    = 1.0 / (hz * hz);
1330:               col[5].i   = ex;
1331:               col[5].j   = ey;
1332:               col[5].k   = ez - 1;
1333:               col[5].loc = ELEMENT;
1334:               col[5].c   = 0;
1335:               valA[5]    = 1.0 / hz;
1336:               col[6].i   = ex;
1337:               col[6].j   = ey;
1338:               col[6].k   = ez;
1339:               col[6].loc = ELEMENT;
1340:               col[6].c   = 0;
1341:               valA[6]    = -1.0 / hz;
1342:             } else {
1343:               nEntries   = 8;
1344:               col[0].i   = ex;
1345:               col[0].j   = ey;
1346:               col[0].k   = ez;
1347:               col[0].loc = BACK;
1348:               col[0].c   = 0;
1349:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1350:               col[1].i   = ex;
1351:               col[1].j   = ey - 1;
1352:               col[1].k   = ez;
1353:               col[1].loc = BACK;
1354:               col[1].c   = 0;
1355:               valA[1]    = 1.0 / (hy * hy);
1356:               col[2].i   = ex;
1357:               col[2].j   = ey + 1;
1358:               col[2].k   = ez;
1359:               col[2].loc = BACK;
1360:               col[2].c   = 0;
1361:               valA[2]    = 1.0 / (hy * hy);
1362:               col[3].i   = ex - 1;
1363:               col[3].j   = ey;
1364:               col[3].k   = ez;
1365:               col[3].loc = BACK;
1366:               col[3].c   = 0;
1367:               valA[3]    = 1.0 / (hx * hx);
1368:               /* Right term missing */
1369:               col[4].i   = ex;
1370:               col[4].j   = ey;
1371:               col[4].k   = ez - 1;
1372:               col[4].loc = BACK;
1373:               col[4].c   = 0;
1374:               valA[4]    = 1.0 / (hz * hz);
1375:               col[5].i   = ex;
1376:               col[5].j   = ey;
1377:               col[5].k   = ez + 1;
1378:               col[5].loc = BACK;
1379:               col[5].c   = 0;
1380:               valA[5]    = 1.0 / (hz * hz);
1381:               col[6].i   = ex;
1382:               col[6].j   = ey;
1383:               col[6].k   = ez - 1;
1384:               col[6].loc = ELEMENT;
1385:               col[6].c   = 0;
1386:               valA[6]    = 1.0 / hz;
1387:               col[7].i   = ex;
1388:               col[7].j   = ey;
1389:               col[7].k   = ez;
1390:               col[7].loc = ELEMENT;
1391:               col[7].c   = 0;
1392:               valA[7]    = -1.0 / hz;
1393:             }
1394:           } else if (ey == 0) {
1395:             nEntries   = 8;
1396:             col[0].i   = ex;
1397:             col[0].j   = ey;
1398:             col[0].k   = ez;
1399:             col[0].loc = BACK;
1400:             col[0].c   = 0;
1401:             valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1402:             /* Down term missing */
1403:             col[1].i   = ex;
1404:             col[1].j   = ey + 1;
1405:             col[1].k   = ez;
1406:             col[1].loc = BACK;
1407:             col[1].c   = 0;
1408:             valA[1]    = 1.0 / (hy * hy);
1409:             col[2].i   = ex - 1;
1410:             col[2].j   = ey;
1411:             col[2].k   = ez;
1412:             col[2].loc = BACK;
1413:             col[2].c   = 0;
1414:             valA[2]    = 1.0 / (hx * hx);
1415:             col[3].i   = ex + 1;
1416:             col[3].j   = ey;
1417:             col[3].k   = ez;
1418:             col[3].loc = BACK;
1419:             col[3].c   = 0;
1420:             valA[3]    = 1.0 / (hx * hx);
1421:             col[4].i   = ex;
1422:             col[4].j   = ey;
1423:             col[4].k   = ez - 1;
1424:             col[4].loc = BACK;
1425:             col[4].c   = 0;
1426:             valA[4]    = 1.0 / (hz * hz);
1427:             col[5].i   = ex;
1428:             col[5].j   = ey;
1429:             col[5].k   = ez + 1;
1430:             col[5].loc = BACK;
1431:             col[5].c   = 0;
1432:             valA[5]    = 1.0 / (hz * hz);
1433:             col[6].i   = ex;
1434:             col[6].j   = ey;
1435:             col[6].k   = ez - 1;
1436:             col[6].loc = ELEMENT;
1437:             col[6].c   = 0;
1438:             valA[6]    = 1.0 / hz;
1439:             col[7].i   = ex;
1440:             col[7].j   = ey;
1441:             col[7].k   = ez;
1442:             col[7].loc = ELEMENT;
1443:             col[7].c   = 0;
1444:             valA[7]    = -1.0 / hz;
1445:           } else if (ey == N[1] - 1) {
1446:             nEntries   = 8;
1447:             col[0].i   = ex;
1448:             col[0].j   = ey;
1449:             col[0].k   = ez;
1450:             col[0].loc = BACK;
1451:             col[0].c   = 0;
1452:             valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1453:             col[1].i   = ex;
1454:             col[1].j   = ey - 1;
1455:             col[1].k   = ez;
1456:             col[1].loc = BACK;
1457:             col[1].c   = 0;
1458:             valA[1]    = 1.0 / (hy * hy);
1459:             /* Up term missing */
1460:             col[2].i   = ex - 1;
1461:             col[2].j   = ey;
1462:             col[2].k   = ez;
1463:             col[2].loc = BACK;
1464:             col[2].c   = 0;
1465:             valA[2]    = 1.0 / (hx * hx);
1466:             col[3].i   = ex + 1;
1467:             col[3].j   = ey;
1468:             col[3].k   = ez;
1469:             col[3].loc = BACK;
1470:             col[3].c   = 0;
1471:             valA[3]    = 1.0 / (hx * hx);
1472:             col[4].i   = ex;
1473:             col[4].j   = ey;
1474:             col[4].k   = ez - 1;
1475:             col[4].loc = BACK;
1476:             col[4].c   = 0;
1477:             valA[4]    = 1.0 / (hz * hz);
1478:             col[5].i   = ex;
1479:             col[5].j   = ey;
1480:             col[5].k   = ez + 1;
1481:             col[5].loc = BACK;
1482:             col[5].c   = 0;
1483:             valA[5]    = 1.0 / (hz * hz);
1484:             col[6].i   = ex;
1485:             col[6].j   = ey;
1486:             col[6].k   = ez - 1;
1487:             col[6].loc = ELEMENT;
1488:             col[6].c   = 0;
1489:             valA[6]    = 1.0 / hz;
1490:             col[7].i   = ex;
1491:             col[7].j   = ey;
1492:             col[7].k   = ez;
1493:             col[7].loc = ELEMENT;
1494:             col[7].c   = 0;
1495:             valA[7]    = -1.0 / hz;
1496:           } else {
1497:             nEntries   = 9;
1498:             col[0].i   = ex;
1499:             col[0].j   = ey;
1500:             col[0].k   = ez;
1501:             col[0].loc = BACK;
1502:             col[0].c   = 0;
1503:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1504:             col[1].i   = ex;
1505:             col[1].j   = ey - 1;
1506:             col[1].k   = ez;
1507:             col[1].loc = BACK;
1508:             col[1].c   = 0;
1509:             valA[1]    = 1.0 / (hy * hy);
1510:             col[2].i   = ex;
1511:             col[2].j   = ey + 1;
1512:             col[2].k   = ez;
1513:             col[2].loc = BACK;
1514:             col[2].c   = 0;
1515:             valA[2]    = 1.0 / (hy * hy);
1516:             col[3].i   = ex - 1;
1517:             col[3].j   = ey;
1518:             col[3].k   = ez;
1519:             col[3].loc = BACK;
1520:             col[3].c   = 0;
1521:             valA[3]    = 1.0 / (hx * hx);
1522:             col[4].i   = ex + 1;
1523:             col[4].j   = ey;
1524:             col[4].k   = ez;
1525:             col[4].loc = BACK;
1526:             col[4].c   = 0;
1527:             valA[4]    = 1.0 / (hx * hx);
1528:             col[5].i   = ex;
1529:             col[5].j   = ey;
1530:             col[5].k   = ez - 1;
1531:             col[5].loc = BACK;
1532:             col[5].c   = 0;
1533:             valA[5]    = 1.0 / (hz * hz);
1534:             col[6].i   = ex;
1535:             col[6].j   = ey;
1536:             col[6].k   = ez + 1;
1537:             col[6].loc = BACK;
1538:             col[6].c   = 0;
1539:             valA[6]    = 1.0 / (hz * hz);
1540:             col[7].i   = ex;
1541:             col[7].j   = ey;
1542:             col[7].k   = ez - 1;
1543:             col[7].loc = ELEMENT;
1544:             col[7].c   = 0;
1545:             valA[7]    = 1.0 / hz;
1546:             col[8].i   = ex;
1547:             col[8].j   = ey;
1548:             col[8].k   = ez;
1549:             col[8].loc = ELEMENT;
1550:             col[8].c   = 0;
1551:             valA[8]    = -1.0 / hz;
1552:           }
1553:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
1554:         }

1556:         /* P equation : u_x + v_y + w_z = g
1557:            Note that this includes an explicit zero on the diagonal. This is only needed for
1558:            direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
1559:         {
1560:           DMStagStencil row, col[7];
1561:           PetscScalar   valA[7];

1563:           row.i      = ex;
1564:           row.j      = ey;
1565:           row.k      = ez;
1566:           row.loc    = ELEMENT;
1567:           row.c      = 0;
1568:           col[0].i   = ex;
1569:           col[0].j   = ey;
1570:           col[0].k   = ez;
1571:           col[0].loc = LEFT;
1572:           col[0].c   = 0;
1573:           valA[0]    = -1.0 / hx;
1574:           col[1].i   = ex;
1575:           col[1].j   = ey;
1576:           col[1].k   = ez;
1577:           col[1].loc = RIGHT;
1578:           col[1].c   = 0;
1579:           valA[1]    = 1.0 / hx;
1580:           col[2].i   = ex;
1581:           col[2].j   = ey;
1582:           col[2].k   = ez;
1583:           col[2].loc = DOWN;
1584:           col[2].c   = 0;
1585:           valA[2]    = -1.0 / hy;
1586:           col[3].i   = ex;
1587:           col[3].j   = ey;
1588:           col[3].k   = ez;
1589:           col[3].loc = UP;
1590:           col[3].c   = 0;
1591:           valA[3]    = 1.0 / hy;
1592:           col[4].i   = ex;
1593:           col[4].j   = ey;
1594:           col[4].k   = ez;
1595:           col[4].loc = BACK;
1596:           col[4].c   = 0;
1597:           valA[4]    = -1.0 / hz;
1598:           col[5].i   = ex;
1599:           col[5].j   = ey;
1600:           col[5].k   = ez;
1601:           col[5].loc = FRONT;
1602:           col[5].c   = 0;
1603:           valA[5]    = 1.0 / hz;
1604:           col[6]     = row;
1605:           valA[6]    = 0.0;
1606:           PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 7, col, valA, INSERT_VALUES));
1607:         }
1608:       }
1609:     }
1610:   }
1611:   PetscCall(DMStagVecRestoreArrayRead(dmCoord, coordLocal, &arrCoord));
1612:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1613:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1614:   PetscFunctionReturn(PETSC_SUCCESS);
1615: }

1617: /* A helper function to check values */
1618: static PetscErrorCode check_vals(PetscInt ex, PetscInt ey, PetscInt ez, PetscInt n, const PetscScalar *ref, const PetscScalar *computed)
1619: {
1620:   PetscInt i;

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

1629: /* The same function as above, but getting and checking values, instead of setting them */
1630: static PetscErrorCode CheckMat(DM dmSol, Mat A)
1631: {
1632:   Vec             coordLocal;
1633:   PetscInt        startx, starty, startz, N[3], nx, ny, nz, ex, ey, ez, d;
1634:   PetscInt        icp[3], icux[3], icuy[3], icuz[3], icux_right[3], icuy_up[3], icuz_front[3];
1635:   PetscBool       isLastRankx, isLastRanky, isLastRankz, isFirstRankx, isFirstRanky, isFirstRankz;
1636:   PetscReal       hx, hy, hz;
1637:   DM              dmCoord;
1638:   PetscScalar ****arrCoord;
1639:   PetscScalar     computed[1024];

1641:   PetscFunctionBeginUser;
1642:   PetscCall(DMStagGetCorners(dmSol, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
1643:   PetscCall(DMStagGetGlobalSizes(dmSol, &N[0], &N[1], &N[2]));
1644:   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");
1645:   PetscCall(DMStagGetIsLastRank(dmSol, &isLastRankx, &isLastRanky, &isLastRankz));
1646:   PetscCall(DMStagGetIsFirstRank(dmSol, &isFirstRankx, &isFirstRanky, &isFirstRankz));
1647:   hx = 1.0 / N[0];
1648:   hy = 1.0 / N[1];
1649:   hz = 1.0 / N[2];
1650:   PetscCall(DMGetCoordinateDM(dmSol, &dmCoord));
1651:   PetscCall(DMGetCoordinatesLocal(dmSol, &coordLocal));
1652:   PetscCall(DMStagVecGetArrayRead(dmCoord, coordLocal, &arrCoord));
1653:   for (d = 0; d < 3; ++d) {
1654:     PetscCall(DMStagGetLocationSlot(dmCoord, ELEMENT, d, &icp[d]));
1655:     PetscCall(DMStagGetLocationSlot(dmCoord, LEFT, d, &icux[d]));
1656:     PetscCall(DMStagGetLocationSlot(dmCoord, DOWN, d, &icuy[d]));
1657:     PetscCall(DMStagGetLocationSlot(dmCoord, BACK, d, &icuz[d]));
1658:     PetscCall(DMStagGetLocationSlot(dmCoord, RIGHT, d, &icux_right[d]));
1659:     PetscCall(DMStagGetLocationSlot(dmCoord, UP, d, &icuy_up[d]));
1660:     PetscCall(DMStagGetLocationSlot(dmCoord, FRONT, d, &icuz_front[d]));
1661:   }

1663:   for (ez = startz; ez < startz + nz; ++ez) {
1664:     for (ey = starty; ey < starty + ny; ++ey) {
1665:       for (ex = startx; ex < startx + nx; ++ex) {
1666:         if (ex == N[0] - 1) {
1667:           /* Right Boundary velocity Dirichlet */
1668:           DMStagStencil     row;
1669:           const PetscScalar valA = 1.0;
1670:           row.i                  = ex;
1671:           row.j                  = ey;
1672:           row.k                  = ez;
1673:           row.loc                = RIGHT;
1674:           row.c                  = 0;
1675:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1676:           PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
1677:         }
1678:         if (ey == N[1] - 1) {
1679:           /* Top boundary velocity Dirichlet */
1680:           DMStagStencil     row;
1681:           const PetscScalar valA = 1.0;
1682:           row.i                  = ex;
1683:           row.j                  = ey;
1684:           row.k                  = ez;
1685:           row.loc                = UP;
1686:           row.c                  = 0;
1687:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1688:           PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
1689:         }
1690:         if (ez == N[2] - 1) {
1691:           /* Top boundary velocity Dirichlet */
1692:           DMStagStencil     row;
1693:           const PetscScalar valA = 1.0;
1694:           row.i                  = ex;
1695:           row.j                  = ey;
1696:           row.k                  = ez;
1697:           row.loc                = FRONT;
1698:           row.c                  = 0;
1699:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1700:           PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
1701:         }

1703:         /* Equation on left face of this element */
1704:         if (ex == 0) {
1705:           /* Left velocity Dirichlet */
1706:           DMStagStencil     row;
1707:           const PetscScalar valA = 1.0;
1708:           row.i                  = ex;
1709:           row.j                  = ey;
1710:           row.k                  = ez;
1711:           row.loc                = LEFT;
1712:           row.c                  = 0;
1713:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
1714:           PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
1715:         } else {
1716:           /* X-momentum interior equation : (u_xx + u_yy + u_zz) - p_x = f^x */
1717:           DMStagStencil row, col[9];
1718:           PetscScalar   valA[9];
1719:           PetscInt      nEntries;

1721:           row.i   = ex;
1722:           row.j   = ey;
1723:           row.k   = ez;
1724:           row.loc = LEFT;
1725:           row.c   = 0;
1726:           if (ey == 0) {
1727:             if (ez == 0) {
1728:               nEntries   = 7;
1729:               col[0].i   = ex;
1730:               col[0].j   = ey;
1731:               col[0].k   = ez;
1732:               col[0].loc = LEFT;
1733:               col[0].c   = 0;
1734:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
1735:               /* Missing down term */
1736:               col[1].i   = ex;
1737:               col[1].j   = ey + 1;
1738:               col[1].k   = ez;
1739:               col[1].loc = LEFT;
1740:               col[1].c   = 0;
1741:               valA[1]    = 1.0 / (hy * hy);
1742:               col[2].i   = ex - 1;
1743:               col[2].j   = ey;
1744:               col[2].k   = ez;
1745:               col[2].loc = LEFT;
1746:               col[2].c   = 0;
1747:               valA[2]    = 1.0 / (hx * hx);
1748:               col[3].i   = ex + 1;
1749:               col[3].j   = ey;
1750:               col[3].k   = ez;
1751:               col[3].loc = LEFT;
1752:               col[3].c   = 0;
1753:               valA[3]    = 1.0 / (hx * hx);
1754:               /* Missing back term */
1755:               col[4].i   = ex;
1756:               col[4].j   = ey;
1757:               col[4].k   = ez + 1;
1758:               col[4].loc = LEFT;
1759:               col[4].c   = 0;
1760:               valA[4]    = 1.0 / (hz * hz);
1761:               col[5].i   = ex - 1;
1762:               col[5].j   = ey;
1763:               col[5].k   = ez;
1764:               col[5].loc = ELEMENT;
1765:               col[5].c   = 0;
1766:               valA[5]    = 1.0 / hx;
1767:               col[6].i   = ex;
1768:               col[6].j   = ey;
1769:               col[6].k   = ez;
1770:               col[6].loc = ELEMENT;
1771:               col[6].c   = 0;
1772:               valA[6]    = -1.0 / hx;
1773:             } else if (ez == N[2] - 1) {
1774:               nEntries   = 7;
1775:               col[0].i   = ex;
1776:               col[0].j   = ey;
1777:               col[0].k   = ez;
1778:               col[0].loc = LEFT;
1779:               col[0].c   = 0;
1780:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
1781:               /* Missing down term */
1782:               col[1].i   = ex;
1783:               col[1].j   = ey + 1;
1784:               col[1].k   = ez;
1785:               col[1].loc = LEFT;
1786:               col[1].c   = 0;
1787:               valA[1]    = 1.0 / (hy * hy);
1788:               col[2].i   = ex - 1;
1789:               col[2].j   = ey;
1790:               col[2].k   = ez;
1791:               col[2].loc = LEFT;
1792:               col[2].c   = 0;
1793:               valA[2]    = 1.0 / (hx * hx);
1794:               col[3].i   = ex + 1;
1795:               col[3].j   = ey;
1796:               col[3].k   = ez;
1797:               col[3].loc = LEFT;
1798:               col[3].c   = 0;
1799:               valA[3]    = 1.0 / (hx * hx);
1800:               col[4].i   = ex;
1801:               col[4].j   = ey;
1802:               col[4].k   = ez - 1;
1803:               col[4].loc = LEFT;
1804:               col[4].c   = 0;
1805:               valA[4]    = 1.0 / (hz * hz);
1806:               /* Missing front term */
1807:               col[5].i   = ex - 1;
1808:               col[5].j   = ey;
1809:               col[5].k   = ez;
1810:               col[5].loc = ELEMENT;
1811:               col[5].c   = 0;
1812:               valA[5]    = 1.0 / hx;
1813:               col[6].i   = ex;
1814:               col[6].j   = ey;
1815:               col[6].k   = ez;
1816:               col[6].loc = ELEMENT;
1817:               col[6].c   = 0;
1818:               valA[6]    = -1.0 / hx;
1819:             } else {
1820:               nEntries   = 8;
1821:               col[0].i   = ex;
1822:               col[0].j   = ey;
1823:               col[0].k   = ez;
1824:               col[0].loc = LEFT;
1825:               col[0].c   = 0;
1826:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1827:               /* Missing down term */
1828:               col[1].i   = ex;
1829:               col[1].j   = ey + 1;
1830:               col[1].k   = ez;
1831:               col[1].loc = LEFT;
1832:               col[1].c   = 0;
1833:               valA[1]    = 1.0 / (hy * hy);
1834:               col[2].i   = ex - 1;
1835:               col[2].j   = ey;
1836:               col[2].k   = ez;
1837:               col[2].loc = LEFT;
1838:               col[2].c   = 0;
1839:               valA[2]    = 1.0 / (hx * hx);
1840:               col[3].i   = ex + 1;
1841:               col[3].j   = ey;
1842:               col[3].k   = ez;
1843:               col[3].loc = LEFT;
1844:               col[3].c   = 0;
1845:               valA[3]    = 1.0 / (hx * hx);
1846:               col[4].i   = ex;
1847:               col[4].j   = ey;
1848:               col[4].k   = ez - 1;
1849:               col[4].loc = LEFT;
1850:               col[4].c   = 0;
1851:               valA[4]    = 1.0 / (hz * hz);
1852:               col[5].i   = ex;
1853:               col[5].j   = ey;
1854:               col[5].k   = ez + 1;
1855:               col[5].loc = LEFT;
1856:               col[5].c   = 0;
1857:               valA[5]    = 1.0 / (hz * hz);
1858:               col[6].i   = ex - 1;
1859:               col[6].j   = ey;
1860:               col[6].k   = ez;
1861:               col[6].loc = ELEMENT;
1862:               col[6].c   = 0;
1863:               valA[6]    = 1.0 / hx;
1864:               col[7].i   = ex;
1865:               col[7].j   = ey;
1866:               col[7].k   = ez;
1867:               col[7].loc = ELEMENT;
1868:               col[7].c   = 0;
1869:               valA[7]    = -1.0 / hx;
1870:             }
1871:           } else if (ey == N[1] - 1) {
1872:             if (ez == 0) {
1873:               nEntries   = 7;
1874:               col[0].i   = ex;
1875:               col[0].j   = ey;
1876:               col[0].k   = ez;
1877:               col[0].loc = LEFT;
1878:               col[0].c   = 0;
1879:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
1880:               col[1].i   = ex;
1881:               col[1].j   = ey - 1;
1882:               col[1].k   = ez;
1883:               col[1].loc = LEFT;
1884:               col[1].c   = 0;
1885:               valA[1]    = 1.0 / (hy * hy);
1886:               /* Missing up term */
1887:               col[2].i   = ex - 1;
1888:               col[2].j   = ey;
1889:               col[2].k   = ez;
1890:               col[2].loc = LEFT;
1891:               col[2].c   = 0;
1892:               valA[2]    = 1.0 / (hx * hx);
1893:               col[3].i   = ex + 1;
1894:               col[3].j   = ey;
1895:               col[3].k   = ez;
1896:               col[3].loc = LEFT;
1897:               col[3].c   = 0;
1898:               valA[3]    = 1.0 / (hx * hx);
1899:               /* Missing back entry */
1900:               col[4].i   = ex;
1901:               col[4].j   = ey;
1902:               col[4].k   = ez + 1;
1903:               col[4].loc = LEFT;
1904:               col[4].c   = 0;
1905:               valA[4]    = 1.0 / (hz * hz);
1906:               col[5].i   = ex - 1;
1907:               col[5].j   = ey;
1908:               col[5].k   = ez;
1909:               col[5].loc = ELEMENT;
1910:               col[5].c   = 0;
1911:               valA[5]    = 1.0 / hx;
1912:               col[6].i   = ex;
1913:               col[6].j   = ey;
1914:               col[6].k   = ez;
1915:               col[6].loc = ELEMENT;
1916:               col[6].c   = 0;
1917:               valA[6]    = -1.0 / hx;
1918:             } else if (ez == N[2] - 1) {
1919:               nEntries   = 7;
1920:               col[0].i   = ex;
1921:               col[0].j   = ey;
1922:               col[0].k   = ez;
1923:               col[0].loc = LEFT;
1924:               col[0].c   = 0;
1925:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
1926:               col[1].i   = ex;
1927:               col[1].j   = ey - 1;
1928:               col[1].k   = ez;
1929:               col[1].loc = LEFT;
1930:               col[1].c   = 0;
1931:               valA[1]    = 1.0 / (hy * hy);
1932:               /* Missing up term */
1933:               col[2].i   = ex - 1;
1934:               col[2].j   = ey;
1935:               col[2].k   = ez;
1936:               col[2].loc = LEFT;
1937:               col[2].c   = 0;
1938:               valA[2]    = 1.0 / (hx * hx);
1939:               col[3].i   = ex + 1;
1940:               col[3].j   = ey;
1941:               col[3].k   = ez;
1942:               col[3].loc = LEFT;
1943:               col[3].c   = 0;
1944:               valA[3]    = 1.0 / (hx * hx);
1945:               col[4].i   = ex;
1946:               col[4].j   = ey;
1947:               col[4].k   = ez - 1;
1948:               col[4].loc = LEFT;
1949:               col[4].c   = 0;
1950:               valA[4]    = 1.0 / (hz * hz);
1951:               /* Missing front term */
1952:               col[5].i   = ex - 1;
1953:               col[5].j   = ey;
1954:               col[5].k   = ez;
1955:               col[5].loc = ELEMENT;
1956:               col[5].c   = 0;
1957:               valA[5]    = 1.0 / hx;
1958:               col[6].i   = ex;
1959:               col[6].j   = ey;
1960:               col[6].k   = ez;
1961:               col[6].loc = ELEMENT;
1962:               col[6].c   = 0;
1963:               valA[6]    = -1.0 / hx;
1964:             } else {
1965:               nEntries   = 8;
1966:               col[0].i   = ex;
1967:               col[0].j   = ey;
1968:               col[0].k   = ez;
1969:               col[0].loc = LEFT;
1970:               col[0].c   = 0;
1971:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1972:               col[1].i   = ex;
1973:               col[1].j   = ey - 1;
1974:               col[1].k   = ez;
1975:               col[1].loc = LEFT;
1976:               col[1].c   = 0;
1977:               valA[1]    = 1.0 / (hy * hy);
1978:               /* Missing up term */
1979:               col[2].i   = ex - 1;
1980:               col[2].j   = ey;
1981:               col[2].k   = ez;
1982:               col[2].loc = LEFT;
1983:               col[2].c   = 0;
1984:               valA[2]    = 1.0 / (hx * hx);
1985:               col[3].i   = ex + 1;
1986:               col[3].j   = ey;
1987:               col[3].k   = ez;
1988:               col[3].loc = LEFT;
1989:               col[3].c   = 0;
1990:               valA[3]    = 1.0 / (hx * hx);
1991:               col[4].i   = ex;
1992:               col[4].j   = ey;
1993:               col[4].k   = ez - 1;
1994:               col[4].loc = LEFT;
1995:               col[4].c   = 0;
1996:               valA[4]    = 1.0 / (hz * hz);
1997:               col[5].i   = ex;
1998:               col[5].j   = ey;
1999:               col[5].k   = ez + 1;
2000:               col[5].loc = LEFT;
2001:               col[5].c   = 0;
2002:               valA[5]    = 1.0 / (hz * hz);
2003:               col[6].i   = ex - 1;
2004:               col[6].j   = ey;
2005:               col[6].k   = ez;
2006:               col[6].loc = ELEMENT;
2007:               col[6].c   = 0;
2008:               valA[6]    = 1.0 / hx;
2009:               col[7].i   = ex;
2010:               col[7].j   = ey;
2011:               col[7].k   = ez;
2012:               col[7].loc = ELEMENT;
2013:               col[7].c   = 0;
2014:               valA[7]    = -1.0 / hx;
2015:             }
2016:           } else if (ez == 0) {
2017:             nEntries   = 8;
2018:             col[0].i   = ex;
2019:             col[0].j   = ey;
2020:             col[0].k   = ez;
2021:             col[0].loc = LEFT;
2022:             col[0].c   = 0;
2023:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2024:             col[1].i   = ex;
2025:             col[1].j   = ey - 1;
2026:             col[1].k   = ez;
2027:             col[1].loc = LEFT;
2028:             col[1].c   = 0;
2029:             valA[1]    = 1.0 / (hy * hy);
2030:             col[2].i   = ex;
2031:             col[2].j   = ey + 1;
2032:             col[2].k   = ez;
2033:             col[2].loc = LEFT;
2034:             col[2].c   = 0;
2035:             valA[2]    = 1.0 / (hy * hy);
2036:             col[3].i   = ex - 1;
2037:             col[3].j   = ey;
2038:             col[3].k   = ez;
2039:             col[3].loc = LEFT;
2040:             col[3].c   = 0;
2041:             valA[3]    = 1.0 / (hx * hx);
2042:             col[4].i   = ex + 1;
2043:             col[4].j   = ey;
2044:             col[4].k   = ez;
2045:             col[4].loc = LEFT;
2046:             col[4].c   = 0;
2047:             valA[4]    = 1.0 / (hx * hx);
2048:             /* Missing back term */
2049:             col[5].i   = ex;
2050:             col[5].j   = ey;
2051:             col[5].k   = ez + 1;
2052:             col[5].loc = LEFT;
2053:             col[5].c   = 0;
2054:             valA[5]    = 1.0 / (hz * hz);
2055:             col[6].i   = ex - 1;
2056:             col[6].j   = ey;
2057:             col[6].k   = ez;
2058:             col[6].loc = ELEMENT;
2059:             col[6].c   = 0;
2060:             valA[6]    = 1.0 / hx;
2061:             col[7].i   = ex;
2062:             col[7].j   = ey;
2063:             col[7].k   = ez;
2064:             col[7].loc = ELEMENT;
2065:             col[7].c   = 0;
2066:             valA[7]    = -1.0 / hx;
2067:           } else if (ez == N[2] - 1) {
2068:             nEntries   = 8;
2069:             col[0].i   = ex;
2070:             col[0].j   = ey;
2071:             col[0].k   = ez;
2072:             col[0].loc = LEFT;
2073:             col[0].c   = 0;
2074:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2075:             col[1].i   = ex;
2076:             col[1].j   = ey - 1;
2077:             col[1].k   = ez;
2078:             col[1].loc = LEFT;
2079:             col[1].c   = 0;
2080:             valA[1]    = 1.0 / (hy * hy);
2081:             col[2].i   = ex;
2082:             col[2].j   = ey + 1;
2083:             col[2].k   = ez;
2084:             col[2].loc = LEFT;
2085:             col[2].c   = 0;
2086:             valA[2]    = 1.0 / (hy * hy);
2087:             col[3].i   = ex - 1;
2088:             col[3].j   = ey;
2089:             col[3].k   = ez;
2090:             col[3].loc = LEFT;
2091:             col[3].c   = 0;
2092:             valA[3]    = 1.0 / (hx * hx);
2093:             col[4].i   = ex + 1;
2094:             col[4].j   = ey;
2095:             col[4].k   = ez;
2096:             col[4].loc = LEFT;
2097:             col[4].c   = 0;
2098:             valA[4]    = 1.0 / (hx * hx);
2099:             col[5].i   = ex;
2100:             col[5].j   = ey;
2101:             col[5].k   = ez - 1;
2102:             col[5].loc = LEFT;
2103:             col[5].c   = 0;
2104:             valA[5]    = 1.0 / (hz * hz);
2105:             /* Missing front term */
2106:             col[6].i   = ex - 1;
2107:             col[6].j   = ey;
2108:             col[6].k   = ez;
2109:             col[6].loc = ELEMENT;
2110:             col[6].c   = 0;
2111:             valA[6]    = 1.0 / hx;
2112:             col[7].i   = ex;
2113:             col[7].j   = ey;
2114:             col[7].k   = ez;
2115:             col[7].loc = ELEMENT;
2116:             col[7].c   = 0;
2117:             valA[7]    = -1.0 / hx;
2118:           } else {
2119:             nEntries   = 9;
2120:             col[0].i   = ex;
2121:             col[0].j   = ey;
2122:             col[0].k   = ez;
2123:             col[0].loc = LEFT;
2124:             col[0].c   = 0;
2125:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2126:             col[1].i   = ex;
2127:             col[1].j   = ey - 1;
2128:             col[1].k   = ez;
2129:             col[1].loc = LEFT;
2130:             col[1].c   = 0;
2131:             valA[1]    = 1.0 / (hy * hy);
2132:             col[2].i   = ex;
2133:             col[2].j   = ey + 1;
2134:             col[2].k   = ez;
2135:             col[2].loc = LEFT;
2136:             col[2].c   = 0;
2137:             valA[2]    = 1.0 / (hy * hy);
2138:             col[3].i   = ex - 1;
2139:             col[3].j   = ey;
2140:             col[3].k   = ez;
2141:             col[3].loc = LEFT;
2142:             col[3].c   = 0;
2143:             valA[3]    = 1.0 / (hx * hx);
2144:             col[4].i   = ex + 1;
2145:             col[4].j   = ey;
2146:             col[4].k   = ez;
2147:             col[4].loc = LEFT;
2148:             col[4].c   = 0;
2149:             valA[4]    = 1.0 / (hx * hx);
2150:             col[5].i   = ex;
2151:             col[5].j   = ey;
2152:             col[5].k   = ez - 1;
2153:             col[5].loc = LEFT;
2154:             col[5].c   = 0;
2155:             valA[5]    = 1.0 / (hz * hz);
2156:             col[6].i   = ex;
2157:             col[6].j   = ey;
2158:             col[6].k   = ez + 1;
2159:             col[6].loc = LEFT;
2160:             col[6].c   = 0;
2161:             valA[6]    = 1.0 / (hz * hz);
2162:             col[7].i   = ex - 1;
2163:             col[7].j   = ey;
2164:             col[7].k   = ez;
2165:             col[7].loc = ELEMENT;
2166:             col[7].c   = 0;
2167:             valA[7]    = 1.0 / hx;
2168:             col[8].i   = ex;
2169:             col[8].j   = ey;
2170:             col[8].k   = ez;
2171:             col[8].loc = ELEMENT;
2172:             col[8].c   = 0;
2173:             valA[8]    = -1.0 / hx;
2174:           }
2175:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed));
2176:           PetscCall(check_vals(ex, ey, ez, nEntries, valA, computed));
2177:         }

2179:         /* Equation on bottom face of this element */
2180:         if (ey == 0) {
2181:           /* Bottom boundary velocity Dirichlet */
2182:           DMStagStencil     row;
2183:           const PetscScalar valA = 1.0;
2184:           row.i                  = ex;
2185:           row.j                  = ey;
2186:           row.k                  = ez;
2187:           row.loc                = DOWN;
2188:           row.c                  = 0;
2189:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
2190:           PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
2191:         } else {
2192:           /* Y-momentum equation, (v_xx + v_yy + v_zz) - p_y = f^y */
2193:           DMStagStencil row, col[9];
2194:           PetscScalar   valA[9];
2195:           PetscInt      nEntries;

2197:           row.i   = ex;
2198:           row.j   = ey;
2199:           row.k   = ez;
2200:           row.loc = DOWN;
2201:           row.c   = 0;
2202:           if (ex == 0) {
2203:             if (ez == 0) {
2204:               nEntries   = 7;
2205:               col[0].i   = ex;
2206:               col[0].j   = ey;
2207:               col[0].k   = ez;
2208:               col[0].loc = DOWN;
2209:               col[0].c   = 0;
2210:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2211:               col[1].i   = ex;
2212:               col[1].j   = ey - 1;
2213:               col[1].k   = ez;
2214:               col[1].loc = DOWN;
2215:               col[1].c   = 0;
2216:               valA[1]    = 1.0 / (hy * hy);
2217:               col[2].i   = ex;
2218:               col[2].j   = ey + 1;
2219:               col[2].k   = ez;
2220:               col[2].loc = DOWN;
2221:               col[2].c   = 0;
2222:               valA[2]    = 1.0 / (hy * hy);
2223:               /* Left term missing */
2224:               col[3].i   = ex + 1;
2225:               col[3].j   = ey;
2226:               col[3].k   = ez;
2227:               col[3].loc = DOWN;
2228:               col[3].c   = 0;
2229:               valA[3]    = 1.0 / (hx * hx);
2230:               /* Back term missing */
2231:               col[4].i   = ex;
2232:               col[4].j   = ey;
2233:               col[4].k   = ez + 1;
2234:               col[4].loc = DOWN;
2235:               col[4].c   = 0;
2236:               valA[4]    = 1.0 / (hz * hz);
2237:               col[5].i   = ex;
2238:               col[5].j   = ey - 1;
2239:               col[5].k   = ez;
2240:               col[5].loc = ELEMENT;
2241:               col[5].c   = 0;
2242:               valA[5]    = 1.0 / hy;
2243:               col[6].i   = ex;
2244:               col[6].j   = ey;
2245:               col[6].k   = ez;
2246:               col[6].loc = ELEMENT;
2247:               col[6].c   = 0;
2248:               valA[6]    = -1.0 / hy;
2249:             } else if (ez == N[2] - 1) {
2250:               nEntries   = 7;
2251:               col[0].i   = ex;
2252:               col[0].j   = ey;
2253:               col[0].k   = ez;
2254:               col[0].loc = DOWN;
2255:               col[0].c   = 0;
2256:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2257:               col[1].i   = ex;
2258:               col[1].j   = ey - 1;
2259:               col[1].k   = ez;
2260:               col[1].loc = DOWN;
2261:               col[1].c   = 0;
2262:               valA[1]    = 1.0 / (hy * hy);
2263:               col[2].i   = ex;
2264:               col[2].j   = ey + 1;
2265:               col[2].k   = ez;
2266:               col[2].loc = DOWN;
2267:               col[2].c   = 0;
2268:               valA[2]    = 1.0 / (hy * hy);
2269:               /* Left term missing */
2270:               col[3].i   = ex + 1;
2271:               col[3].j   = ey;
2272:               col[3].k   = ez;
2273:               col[3].loc = DOWN;
2274:               col[3].c   = 0;
2275:               valA[3]    = 1.0 / (hx * hx);
2276:               col[4].i   = ex;
2277:               col[4].j   = ey;
2278:               col[4].k   = ez - 1;
2279:               col[4].loc = DOWN;
2280:               col[4].c   = 0;
2281:               valA[4]    = 1.0 / (hz * hz);
2282:               /* Front term missing */
2283:               col[5].i   = ex;
2284:               col[5].j   = ey - 1;
2285:               col[5].k   = ez;
2286:               col[5].loc = ELEMENT;
2287:               col[5].c   = 0;
2288:               valA[5]    = 1.0 / hy;
2289:               col[6].i   = ex;
2290:               col[6].j   = ey;
2291:               col[6].k   = ez;
2292:               col[6].loc = ELEMENT;
2293:               col[6].c   = 0;
2294:               valA[6]    = -1.0 / hy;
2295:             } else {
2296:               nEntries   = 8;
2297:               col[0].i   = ex;
2298:               col[0].j   = ey;
2299:               col[0].k   = ez;
2300:               col[0].loc = DOWN;
2301:               col[0].c   = 0;
2302:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2303:               col[1].i   = ex;
2304:               col[1].j   = ey - 1;
2305:               col[1].k   = ez;
2306:               col[1].loc = DOWN;
2307:               col[1].c   = 0;
2308:               valA[1]    = 1.0 / (hy * hy);
2309:               col[2].i   = ex;
2310:               col[2].j   = ey + 1;
2311:               col[2].k   = ez;
2312:               col[2].loc = DOWN;
2313:               col[2].c   = 0;
2314:               valA[2]    = 1.0 / (hy * hy);
2315:               /* Left term missing */
2316:               col[3].i   = ex + 1;
2317:               col[3].j   = ey;
2318:               col[3].k   = ez;
2319:               col[3].loc = DOWN;
2320:               col[3].c   = 0;
2321:               valA[3]    = 1.0 / (hx * hx);
2322:               col[4].i   = ex;
2323:               col[4].j   = ey;
2324:               col[4].k   = ez - 1;
2325:               col[4].loc = DOWN;
2326:               col[4].c   = 0;
2327:               valA[4]    = 1.0 / (hz * hz);
2328:               col[5].i   = ex;
2329:               col[5].j   = ey;
2330:               col[5].k   = ez + 1;
2331:               col[5].loc = DOWN;
2332:               col[5].c   = 0;
2333:               valA[5]    = 1.0 / (hz * hz);
2334:               col[6].i   = ex;
2335:               col[6].j   = ey - 1;
2336:               col[6].k   = ez;
2337:               col[6].loc = ELEMENT;
2338:               col[6].c   = 0;
2339:               valA[6]    = 1.0 / hy;
2340:               col[7].i   = ex;
2341:               col[7].j   = ey;
2342:               col[7].k   = ez;
2343:               col[7].loc = ELEMENT;
2344:               col[7].c   = 0;
2345:               valA[7]    = -1.0 / hy;
2346:             }
2347:           } else if (ex == N[0] - 1) {
2348:             if (ez == 0) {
2349:               nEntries   = 7;
2350:               col[0].i   = ex;
2351:               col[0].j   = ey;
2352:               col[0].k   = ez;
2353:               col[0].loc = DOWN;
2354:               col[0].c   = 0;
2355:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2356:               col[1].i   = ex;
2357:               col[1].j   = ey - 1;
2358:               col[1].k   = ez;
2359:               col[1].loc = DOWN;
2360:               col[1].c   = 0;
2361:               valA[1]    = 1.0 / (hy * hy);
2362:               col[2].i   = ex;
2363:               col[2].j   = ey + 1;
2364:               col[2].k   = ez;
2365:               col[2].loc = DOWN;
2366:               col[2].c   = 0;
2367:               valA[2]    = 1.0 / (hy * hy);
2368:               col[3].i   = ex - 1;
2369:               col[3].j   = ey;
2370:               col[3].k   = ez;
2371:               col[3].loc = DOWN;
2372:               col[3].c   = 0;
2373:               valA[3]    = 1.0 / (hx * hx);
2374:               /* Right term missing */
2375:               /* Back term missing */
2376:               col[4].i   = ex;
2377:               col[4].j   = ey;
2378:               col[4].k   = ez + 1;
2379:               col[4].loc = DOWN;
2380:               col[4].c   = 0;
2381:               valA[4]    = 1.0 / (hz * hz);
2382:               col[5].i   = ex;
2383:               col[5].j   = ey - 1;
2384:               col[5].k   = ez;
2385:               col[5].loc = ELEMENT;
2386:               col[5].c   = 0;
2387:               valA[5]    = 1.0 / hy;
2388:               col[6].i   = ex;
2389:               col[6].j   = ey;
2390:               col[6].k   = ez;
2391:               col[6].loc = ELEMENT;
2392:               col[6].c   = 0;
2393:               valA[6]    = -1.0 / hy;
2394:             } else if (ez == N[2] - 1) {
2395:               nEntries   = 7;
2396:               col[0].i   = ex;
2397:               col[0].j   = ey;
2398:               col[0].k   = ez;
2399:               col[0].loc = DOWN;
2400:               col[0].c   = 0;
2401:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2402:               col[1].i   = ex;
2403:               col[1].j   = ey - 1;
2404:               col[1].k   = ez;
2405:               col[1].loc = DOWN;
2406:               col[1].c   = 0;
2407:               valA[1]    = 1.0 / (hy * hy);
2408:               col[2].i   = ex;
2409:               col[2].j   = ey + 1;
2410:               col[2].k   = ez;
2411:               col[2].loc = DOWN;
2412:               col[2].c   = 0;
2413:               valA[2]    = 1.0 / (hy * hy);
2414:               col[3].i   = ex - 1;
2415:               col[3].j   = ey;
2416:               col[3].k   = ez;
2417:               col[3].loc = DOWN;
2418:               col[3].c   = 0;
2419:               valA[3]    = 1.0 / (hx * hx);
2420:               /* Right term missing */
2421:               col[4].i   = ex;
2422:               col[4].j   = ey;
2423:               col[4].k   = ez - 1;
2424:               col[4].loc = DOWN;
2425:               col[4].c   = 0;
2426:               valA[4]    = 1.0 / (hz * hz);
2427:               /* Front term missing */
2428:               col[5].i   = ex;
2429:               col[5].j   = ey - 1;
2430:               col[5].k   = ez;
2431:               col[5].loc = ELEMENT;
2432:               col[5].c   = 0;
2433:               valA[5]    = 1.0 / hy;
2434:               col[6].i   = ex;
2435:               col[6].j   = ey;
2436:               col[6].k   = ez;
2437:               col[6].loc = ELEMENT;
2438:               col[6].c   = 0;
2439:               valA[6]    = -1.0 / hy;
2440:             } else {
2441:               nEntries   = 8;
2442:               col[0].i   = ex;
2443:               col[0].j   = ey;
2444:               col[0].k   = ez;
2445:               col[0].loc = DOWN;
2446:               col[0].c   = 0;
2447:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2448:               col[1].i   = ex;
2449:               col[1].j   = ey - 1;
2450:               col[1].k   = ez;
2451:               col[1].loc = DOWN;
2452:               col[1].c   = 0;
2453:               valA[1]    = 1.0 / (hy * hy);
2454:               col[2].i   = ex;
2455:               col[2].j   = ey + 1;
2456:               col[2].k   = ez;
2457:               col[2].loc = DOWN;
2458:               col[2].c   = 0;
2459:               valA[2]    = 1.0 / (hy * hy);
2460:               col[3].i   = ex - 1;
2461:               col[3].j   = ey;
2462:               col[3].k   = ez;
2463:               col[3].loc = DOWN;
2464:               col[3].c   = 0;
2465:               valA[3]    = 1.0 / (hx * hx);
2466:               /* Right term missing */
2467:               col[4].i   = ex;
2468:               col[4].j   = ey;
2469:               col[4].k   = ez - 1;
2470:               col[4].loc = DOWN;
2471:               col[4].c   = 0;
2472:               valA[4]    = 1.0 / (hz * hz);
2473:               col[5].i   = ex;
2474:               col[5].j   = ey;
2475:               col[5].k   = ez + 1;
2476:               col[5].loc = DOWN;
2477:               col[5].c   = 0;
2478:               valA[5]    = 1.0 / (hz * hz);
2479:               col[6].i   = ex;
2480:               col[6].j   = ey - 1;
2481:               col[6].k   = ez;
2482:               col[6].loc = ELEMENT;
2483:               col[6].c   = 0;
2484:               valA[6]    = 1.0 / hy;
2485:               col[7].i   = ex;
2486:               col[7].j   = ey;
2487:               col[7].k   = ez;
2488:               col[7].loc = ELEMENT;
2489:               col[7].c   = 0;
2490:               valA[7]    = -1.0 / hy;
2491:             }
2492:           } else if (ez == 0) {
2493:             nEntries   = 8;
2494:             col[0].i   = ex;
2495:             col[0].j   = ey;
2496:             col[0].k   = ez;
2497:             col[0].loc = DOWN;
2498:             col[0].c   = 0;
2499:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2500:             col[1].i   = ex;
2501:             col[1].j   = ey - 1;
2502:             col[1].k   = ez;
2503:             col[1].loc = DOWN;
2504:             col[1].c   = 0;
2505:             valA[1]    = 1.0 / (hy * hy);
2506:             col[2].i   = ex;
2507:             col[2].j   = ey + 1;
2508:             col[2].k   = ez;
2509:             col[2].loc = DOWN;
2510:             col[2].c   = 0;
2511:             valA[2]    = 1.0 / (hy * hy);
2512:             col[3].i   = ex - 1;
2513:             col[3].j   = ey;
2514:             col[3].k   = ez;
2515:             col[3].loc = DOWN;
2516:             col[3].c   = 0;
2517:             valA[3]    = 1.0 / (hx * hx);
2518:             col[4].i   = ex + 1;
2519:             col[4].j   = ey;
2520:             col[4].k   = ez;
2521:             col[4].loc = DOWN;
2522:             col[4].c   = 0;
2523:             valA[4]    = 1.0 / (hx * hx);
2524:             /* Back term missing */
2525:             col[5].i   = ex;
2526:             col[5].j   = ey;
2527:             col[5].k   = ez + 1;
2528:             col[5].loc = DOWN;
2529:             col[5].c   = 0;
2530:             valA[5]    = 1.0 / (hz * hz);
2531:             col[6].i   = ex;
2532:             col[6].j   = ey - 1;
2533:             col[6].k   = ez;
2534:             col[6].loc = ELEMENT;
2535:             col[6].c   = 0;
2536:             valA[6]    = 1.0 / hy;
2537:             col[7].i   = ex;
2538:             col[7].j   = ey;
2539:             col[7].k   = ez;
2540:             col[7].loc = ELEMENT;
2541:             col[7].c   = 0;
2542:             valA[7]    = -1.0 / hy;
2543:           } else if (ez == N[2] - 1) {
2544:             nEntries   = 8;
2545:             col[0].i   = ex;
2546:             col[0].j   = ey;
2547:             col[0].k   = ez;
2548:             col[0].loc = DOWN;
2549:             col[0].c   = 0;
2550:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
2551:             col[1].i   = ex;
2552:             col[1].j   = ey - 1;
2553:             col[1].k   = ez;
2554:             col[1].loc = DOWN;
2555:             col[1].c   = 0;
2556:             valA[1]    = 1.0 / (hy * hy);
2557:             col[2].i   = ex;
2558:             col[2].j   = ey + 1;
2559:             col[2].k   = ez;
2560:             col[2].loc = DOWN;
2561:             col[2].c   = 0;
2562:             valA[2]    = 1.0 / (hy * hy);
2563:             col[3].i   = ex - 1;
2564:             col[3].j   = ey;
2565:             col[3].k   = ez;
2566:             col[3].loc = DOWN;
2567:             col[3].c   = 0;
2568:             valA[3]    = 1.0 / (hx * hx);
2569:             col[4].i   = ex + 1;
2570:             col[4].j   = ey;
2571:             col[4].k   = ez;
2572:             col[4].loc = DOWN;
2573:             col[4].c   = 0;
2574:             valA[4]    = 1.0 / (hx * hx);
2575:             col[5].i   = ex;
2576:             col[5].j   = ey;
2577:             col[5].k   = ez - 1;
2578:             col[5].loc = DOWN;
2579:             col[5].c   = 0;
2580:             valA[5]    = 1.0 / (hz * hz);
2581:             /* Front term missing */
2582:             col[6].i   = ex;
2583:             col[6].j   = ey - 1;
2584:             col[6].k   = ez;
2585:             col[6].loc = ELEMENT;
2586:             col[6].c   = 0;
2587:             valA[6]    = 1.0 / hy;
2588:             col[7].i   = ex;
2589:             col[7].j   = ey;
2590:             col[7].k   = ez;
2591:             col[7].loc = ELEMENT;
2592:             col[7].c   = 0;
2593:             valA[7]    = -1.0 / hy;
2594:           } else {
2595:             nEntries   = 9;
2596:             col[0].i   = ex;
2597:             col[0].j   = ey;
2598:             col[0].k   = ez;
2599:             col[0].loc = DOWN;
2600:             col[0].c   = 0;
2601:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2602:             col[1].i   = ex;
2603:             col[1].j   = ey - 1;
2604:             col[1].k   = ez;
2605:             col[1].loc = DOWN;
2606:             col[1].c   = 0;
2607:             valA[1]    = 1.0 / (hy * hy);
2608:             col[2].i   = ex;
2609:             col[2].j   = ey + 1;
2610:             col[2].k   = ez;
2611:             col[2].loc = DOWN;
2612:             col[2].c   = 0;
2613:             valA[2]    = 1.0 / (hy * hy);
2614:             col[3].i   = ex - 1;
2615:             col[3].j   = ey;
2616:             col[3].k   = ez;
2617:             col[3].loc = DOWN;
2618:             col[3].c   = 0;
2619:             valA[3]    = 1.0 / (hx * hx);
2620:             col[4].i   = ex + 1;
2621:             col[4].j   = ey;
2622:             col[4].k   = ez;
2623:             col[4].loc = DOWN;
2624:             col[4].c   = 0;
2625:             valA[4]    = 1.0 / (hx * hx);
2626:             col[5].i   = ex;
2627:             col[5].j   = ey;
2628:             col[5].k   = ez - 1;
2629:             col[5].loc = DOWN;
2630:             col[5].c   = 0;
2631:             valA[5]    = 1.0 / (hz * hz);
2632:             col[6].i   = ex;
2633:             col[6].j   = ey;
2634:             col[6].k   = ez + 1;
2635:             col[6].loc = DOWN;
2636:             col[6].c   = 0;
2637:             valA[6]    = 1.0 / (hz * hz);
2638:             col[7].i   = ex;
2639:             col[7].j   = ey - 1;
2640:             col[7].k   = ez;
2641:             col[7].loc = ELEMENT;
2642:             col[7].c   = 0;
2643:             valA[7]    = 1.0 / hy;
2644:             col[8].i   = ex;
2645:             col[8].j   = ey;
2646:             col[8].k   = ez;
2647:             col[8].loc = ELEMENT;
2648:             col[8].c   = 0;
2649:             valA[8]    = -1.0 / hy;
2650:           }
2651:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed));
2652:           PetscCall(check_vals(ex, ey, ez, nEntries, valA, computed));
2653:         }

2655:         /* Equation on back face of this element */
2656:         if (ez == 0) {
2657:           /* Back boundary velocity Dirichlet */
2658:           DMStagStencil     row;
2659:           const PetscScalar valA = 1.0;
2660:           row.i                  = ex;
2661:           row.j                  = ey;
2662:           row.k                  = ez;
2663:           row.loc                = BACK;
2664:           row.c                  = 0;
2665:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 1, &row, computed));
2666:           PetscCall(check_vals(ex, ey, ez, 1, &valA, computed));
2667:         } else {
2668:           /* Z-momentum equation, (w_xx + w_yy + w_zz) - p_z = f^z */
2669:           DMStagStencil row, col[9];
2670:           PetscScalar   valA[9];
2671:           PetscInt      nEntries;

2673:           row.i   = ex;
2674:           row.j   = ey;
2675:           row.k   = ez;
2676:           row.loc = BACK;
2677:           row.c   = 0;
2678:           if (ex == 0) {
2679:             if (ey == 0) {
2680:               nEntries   = 7;
2681:               col[0].i   = ex;
2682:               col[0].j   = ey;
2683:               col[0].k   = ez;
2684:               col[0].loc = BACK;
2685:               col[0].c   = 0;
2686:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2687:               /* Down term missing */
2688:               col[1].i   = ex;
2689:               col[1].j   = ey + 1;
2690:               col[1].k   = ez;
2691:               col[1].loc = BACK;
2692:               col[1].c   = 0;
2693:               valA[1]    = 1.0 / (hy * hy);
2694:               /* Left term missing */
2695:               col[2].i   = ex + 1;
2696:               col[2].j   = ey;
2697:               col[2].k   = ez;
2698:               col[2].loc = BACK;
2699:               col[2].c   = 0;
2700:               valA[2]    = 1.0 / (hx * hx);
2701:               col[3].i   = ex;
2702:               col[3].j   = ey;
2703:               col[3].k   = ez - 1;
2704:               col[3].loc = BACK;
2705:               col[3].c   = 0;
2706:               valA[3]    = 1.0 / (hz * hz);
2707:               col[4].i   = ex;
2708:               col[4].j   = ey;
2709:               col[4].k   = ez + 1;
2710:               col[4].loc = BACK;
2711:               col[4].c   = 0;
2712:               valA[4]    = 1.0 / (hz * hz);
2713:               col[5].i   = ex;
2714:               col[5].j   = ey;
2715:               col[5].k   = ez - 1;
2716:               col[5].loc = ELEMENT;
2717:               col[5].c   = 0;
2718:               valA[5]    = 1.0 / hz;
2719:               col[6].i   = ex;
2720:               col[6].j   = ey;
2721:               col[6].k   = ez;
2722:               col[6].loc = ELEMENT;
2723:               col[6].c   = 0;
2724:               valA[6]    = -1.0 / hz;
2725:             } else if (ey == N[1] - 1) {
2726:               nEntries   = 7;
2727:               col[0].i   = ex;
2728:               col[0].j   = ey;
2729:               col[0].k   = ez;
2730:               col[0].loc = BACK;
2731:               col[0].c   = 0;
2732:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2733:               col[1].i   = ex;
2734:               col[1].j   = ey - 1;
2735:               col[1].k   = ez;
2736:               col[1].loc = BACK;
2737:               col[1].c   = 0;
2738:               valA[1]    = 1.0 / (hy * hy);
2739:               /* Up term missing */
2740:               /* Left term missing */
2741:               col[2].i   = ex + 1;
2742:               col[2].j   = ey;
2743:               col[2].k   = ez;
2744:               col[2].loc = BACK;
2745:               col[2].c   = 0;
2746:               valA[2]    = 1.0 / (hx * hx);
2747:               col[3].i   = ex;
2748:               col[3].j   = ey;
2749:               col[3].k   = ez - 1;
2750:               col[3].loc = BACK;
2751:               col[3].c   = 0;
2752:               valA[3]    = 1.0 / (hz * hz);
2753:               col[4].i   = ex;
2754:               col[4].j   = ey;
2755:               col[4].k   = ez + 1;
2756:               col[4].loc = BACK;
2757:               col[4].c   = 0;
2758:               valA[4]    = 1.0 / (hz * hz);
2759:               col[5].i   = ex;
2760:               col[5].j   = ey;
2761:               col[5].k   = ez - 1;
2762:               col[5].loc = ELEMENT;
2763:               col[5].c   = 0;
2764:               valA[5]    = 1.0 / hz;
2765:               col[6].i   = ex;
2766:               col[6].j   = ey;
2767:               col[6].k   = ez;
2768:               col[6].loc = ELEMENT;
2769:               col[6].c   = 0;
2770:               valA[6]    = -1.0 / hz;
2771:             } else {
2772:               nEntries   = 8;
2773:               col[0].i   = ex;
2774:               col[0].j   = ey;
2775:               col[0].k   = ez;
2776:               col[0].loc = BACK;
2777:               col[0].c   = 0;
2778:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2779:               col[1].i   = ex;
2780:               col[1].j   = ey - 1;
2781:               col[1].k   = ez;
2782:               col[1].loc = BACK;
2783:               col[1].c   = 0;
2784:               valA[1]    = 1.0 / (hy * hy);
2785:               col[2].i   = ex;
2786:               col[2].j   = ey + 1;
2787:               col[2].k   = ez;
2788:               col[2].loc = BACK;
2789:               col[2].c   = 0;
2790:               valA[2]    = 1.0 / (hy * hy);
2791:               /* Left term missing */
2792:               col[3].i   = ex + 1;
2793:               col[3].j   = ey;
2794:               col[3].k   = ez;
2795:               col[3].loc = BACK;
2796:               col[3].c   = 0;
2797:               valA[3]    = 1.0 / (hx * hx);
2798:               col[4].i   = ex;
2799:               col[4].j   = ey;
2800:               col[4].k   = ez - 1;
2801:               col[4].loc = BACK;
2802:               col[4].c   = 0;
2803:               valA[4]    = 1.0 / (hz * hz);
2804:               col[5].i   = ex;
2805:               col[5].j   = ey;
2806:               col[5].k   = ez + 1;
2807:               col[5].loc = BACK;
2808:               col[5].c   = 0;
2809:               valA[5]    = 1.0 / (hz * hz);
2810:               col[6].i   = ex;
2811:               col[6].j   = ey;
2812:               col[6].k   = ez - 1;
2813:               col[6].loc = ELEMENT;
2814:               col[6].c   = 0;
2815:               valA[6]    = 1.0 / hz;
2816:               col[7].i   = ex;
2817:               col[7].j   = ey;
2818:               col[7].k   = ez;
2819:               col[7].loc = ELEMENT;
2820:               col[7].c   = 0;
2821:               valA[7]    = -1.0 / hz;
2822:             }
2823:           } else if (ex == N[0] - 1) {
2824:             if (ey == 0) {
2825:               nEntries   = 7;
2826:               col[0].i   = ex;
2827:               col[0].j   = ey;
2828:               col[0].k   = ez;
2829:               col[0].loc = BACK;
2830:               col[0].c   = 0;
2831:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2832:               /* Down term missing */
2833:               col[1].i   = ex;
2834:               col[1].j   = ey + 1;
2835:               col[1].k   = ez;
2836:               col[1].loc = BACK;
2837:               col[1].c   = 0;
2838:               valA[1]    = 1.0 / (hy * hy);
2839:               col[2].i   = ex - 1;
2840:               col[2].j   = ey;
2841:               col[2].k   = ez;
2842:               col[2].loc = BACK;
2843:               col[2].c   = 0;
2844:               valA[2]    = 1.0 / (hx * hx);
2845:               /* Right term missing */
2846:               col[3].i   = ex;
2847:               col[3].j   = ey;
2848:               col[3].k   = ez - 1;
2849:               col[3].loc = BACK;
2850:               col[3].c   = 0;
2851:               valA[3]    = 1.0 / (hz * hz);
2852:               col[4].i   = ex;
2853:               col[4].j   = ey;
2854:               col[4].k   = ez + 1;
2855:               col[4].loc = BACK;
2856:               col[4].c   = 0;
2857:               valA[4]    = 1.0 / (hz * hz);
2858:               col[5].i   = ex;
2859:               col[5].j   = ey;
2860:               col[5].k   = ez - 1;
2861:               col[5].loc = ELEMENT;
2862:               col[5].c   = 0;
2863:               valA[5]    = 1.0 / hz;
2864:               col[6].i   = ex;
2865:               col[6].j   = ey;
2866:               col[6].k   = ez;
2867:               col[6].loc = ELEMENT;
2868:               col[6].c   = 0;
2869:               valA[6]    = -1.0 / hz;
2870:             } else if (ey == N[1] - 1) {
2871:               nEntries   = 7;
2872:               col[0].i   = ex;
2873:               col[0].j   = ey;
2874:               col[0].k   = ez;
2875:               col[0].loc = BACK;
2876:               col[0].c   = 0;
2877:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2878:               col[1].i   = ex;
2879:               col[1].j   = ey - 1;
2880:               col[1].k   = ez;
2881:               col[1].loc = BACK;
2882:               col[1].c   = 0;
2883:               valA[1]    = 1.0 / (hy * hy);
2884:               /* Up term missing */
2885:               col[2].i   = ex - 1;
2886:               col[2].j   = ey;
2887:               col[2].k   = ez;
2888:               col[2].loc = BACK;
2889:               col[2].c   = 0;
2890:               valA[2]    = 1.0 / (hx * hx);
2891:               /* Right term missing */
2892:               col[3].i   = ex;
2893:               col[3].j   = ey;
2894:               col[3].k   = ez - 1;
2895:               col[3].loc = BACK;
2896:               col[3].c   = 0;
2897:               valA[3]    = 1.0 / (hz * hz);
2898:               col[4].i   = ex;
2899:               col[4].j   = ey;
2900:               col[4].k   = ez + 1;
2901:               col[4].loc = BACK;
2902:               col[4].c   = 0;
2903:               valA[4]    = 1.0 / (hz * hz);
2904:               col[5].i   = ex;
2905:               col[5].j   = ey;
2906:               col[5].k   = ez - 1;
2907:               col[5].loc = ELEMENT;
2908:               col[5].c   = 0;
2909:               valA[5]    = 1.0 / hz;
2910:               col[6].i   = ex;
2911:               col[6].j   = ey;
2912:               col[6].k   = ez;
2913:               col[6].loc = ELEMENT;
2914:               col[6].c   = 0;
2915:               valA[6]    = -1.0 / hz;
2916:             } else {
2917:               nEntries   = 8;
2918:               col[0].i   = ex;
2919:               col[0].j   = ey;
2920:               col[0].k   = ez;
2921:               col[0].loc = BACK;
2922:               col[0].c   = 0;
2923:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
2924:               col[1].i   = ex;
2925:               col[1].j   = ey - 1;
2926:               col[1].k   = ez;
2927:               col[1].loc = BACK;
2928:               col[1].c   = 0;
2929:               valA[1]    = 1.0 / (hy * hy);
2930:               col[2].i   = ex;
2931:               col[2].j   = ey + 1;
2932:               col[2].k   = ez;
2933:               col[2].loc = BACK;
2934:               col[2].c   = 0;
2935:               valA[2]    = 1.0 / (hy * hy);
2936:               col[3].i   = ex - 1;
2937:               col[3].j   = ey;
2938:               col[3].k   = ez;
2939:               col[3].loc = BACK;
2940:               col[3].c   = 0;
2941:               valA[3]    = 1.0 / (hx * hx);
2942:               /* Right term missing */
2943:               col[4].i   = ex;
2944:               col[4].j   = ey;
2945:               col[4].k   = ez - 1;
2946:               col[4].loc = BACK;
2947:               col[4].c   = 0;
2948:               valA[4]    = 1.0 / (hz * hz);
2949:               col[5].i   = ex;
2950:               col[5].j   = ey;
2951:               col[5].k   = ez + 1;
2952:               col[5].loc = BACK;
2953:               col[5].c   = 0;
2954:               valA[5]    = 1.0 / (hz * hz);
2955:               col[6].i   = ex;
2956:               col[6].j   = ey;
2957:               col[6].k   = ez - 1;
2958:               col[6].loc = ELEMENT;
2959:               col[6].c   = 0;
2960:               valA[6]    = 1.0 / hz;
2961:               col[7].i   = ex;
2962:               col[7].j   = ey;
2963:               col[7].k   = ez;
2964:               col[7].loc = ELEMENT;
2965:               col[7].c   = 0;
2966:               valA[7]    = -1.0 / hz;
2967:             }
2968:           } else if (ey == 0) {
2969:             nEntries   = 8;
2970:             col[0].i   = ex;
2971:             col[0].j   = ey;
2972:             col[0].k   = ez;
2973:             col[0].loc = BACK;
2974:             col[0].c   = 0;
2975:             valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
2976:             /* Down term missing */
2977:             col[1].i   = ex;
2978:             col[1].j   = ey + 1;
2979:             col[1].k   = ez;
2980:             col[1].loc = BACK;
2981:             col[1].c   = 0;
2982:             valA[1]    = 1.0 / (hy * hy);
2983:             col[2].i   = ex - 1;
2984:             col[2].j   = ey;
2985:             col[2].k   = ez;
2986:             col[2].loc = BACK;
2987:             col[2].c   = 0;
2988:             valA[2]    = 1.0 / (hx * hx);
2989:             col[3].i   = ex + 1;
2990:             col[3].j   = ey;
2991:             col[3].k   = ez;
2992:             col[3].loc = BACK;
2993:             col[3].c   = 0;
2994:             valA[3]    = 1.0 / (hx * hx);
2995:             col[4].i   = ex;
2996:             col[4].j   = ey;
2997:             col[4].k   = ez - 1;
2998:             col[4].loc = BACK;
2999:             col[4].c   = 0;
3000:             valA[4]    = 1.0 / (hz * hz);
3001:             col[5].i   = ex;
3002:             col[5].j   = ey;
3003:             col[5].k   = ez + 1;
3004:             col[5].loc = BACK;
3005:             col[5].c   = 0;
3006:             valA[5]    = 1.0 / (hz * hz);
3007:             col[6].i   = ex;
3008:             col[6].j   = ey;
3009:             col[6].k   = ez - 1;
3010:             col[6].loc = ELEMENT;
3011:             col[6].c   = 0;
3012:             valA[6]    = 1.0 / hz;
3013:             col[7].i   = ex;
3014:             col[7].j   = ey;
3015:             col[7].k   = ez;
3016:             col[7].loc = ELEMENT;
3017:             col[7].c   = 0;
3018:             valA[7]    = -1.0 / hz;
3019:           } else if (ey == N[1] - 1) {
3020:             nEntries   = 8;
3021:             col[0].i   = ex;
3022:             col[0].j   = ey;
3023:             col[0].k   = ez;
3024:             col[0].loc = BACK;
3025:             col[0].c   = 0;
3026:             valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
3027:             col[1].i   = ex;
3028:             col[1].j   = ey - 1;
3029:             col[1].k   = ez;
3030:             col[1].loc = BACK;
3031:             col[1].c   = 0;
3032:             valA[1]    = 1.0 / (hy * hy);
3033:             /* Up term missing */
3034:             col[2].i   = ex - 1;
3035:             col[2].j   = ey;
3036:             col[2].k   = ez;
3037:             col[2].loc = BACK;
3038:             col[2].c   = 0;
3039:             valA[2]    = 1.0 / (hx * hx);
3040:             col[3].i   = ex + 1;
3041:             col[3].j   = ey;
3042:             col[3].k   = ez;
3043:             col[3].loc = BACK;
3044:             col[3].c   = 0;
3045:             valA[3]    = 1.0 / (hx * hx);
3046:             col[4].i   = ex;
3047:             col[4].j   = ey;
3048:             col[4].k   = ez - 1;
3049:             col[4].loc = BACK;
3050:             col[4].c   = 0;
3051:             valA[4]    = 1.0 / (hz * hz);
3052:             col[5].i   = ex;
3053:             col[5].j   = ey;
3054:             col[5].k   = ez + 1;
3055:             col[5].loc = BACK;
3056:             col[5].c   = 0;
3057:             valA[5]    = 1.0 / (hz * hz);
3058:             col[6].i   = ex;
3059:             col[6].j   = ey;
3060:             col[6].k   = ez - 1;
3061:             col[6].loc = ELEMENT;
3062:             col[6].c   = 0;
3063:             valA[6]    = 1.0 / hz;
3064:             col[7].i   = ex;
3065:             col[7].j   = ey;
3066:             col[7].k   = ez;
3067:             col[7].loc = ELEMENT;
3068:             col[7].c   = 0;
3069:             valA[7]    = -1.0 / hz;
3070:           } else {
3071:             nEntries   = 9;
3072:             col[0].i   = ex;
3073:             col[0].j   = ey;
3074:             col[0].k   = ez;
3075:             col[0].loc = BACK;
3076:             col[0].c   = 0;
3077:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
3078:             col[1].i   = ex;
3079:             col[1].j   = ey - 1;
3080:             col[1].k   = ez;
3081:             col[1].loc = BACK;
3082:             col[1].c   = 0;
3083:             valA[1]    = 1.0 / (hy * hy);
3084:             col[2].i   = ex;
3085:             col[2].j   = ey + 1;
3086:             col[2].k   = ez;
3087:             col[2].loc = BACK;
3088:             col[2].c   = 0;
3089:             valA[2]    = 1.0 / (hy * hy);
3090:             col[3].i   = ex - 1;
3091:             col[3].j   = ey;
3092:             col[3].k   = ez;
3093:             col[3].loc = BACK;
3094:             col[3].c   = 0;
3095:             valA[3]    = 1.0 / (hx * hx);
3096:             col[4].i   = ex + 1;
3097:             col[4].j   = ey;
3098:             col[4].k   = ez;
3099:             col[4].loc = BACK;
3100:             col[4].c   = 0;
3101:             valA[4]    = 1.0 / (hx * hx);
3102:             col[5].i   = ex;
3103:             col[5].j   = ey;
3104:             col[5].k   = ez - 1;
3105:             col[5].loc = BACK;
3106:             col[5].c   = 0;
3107:             valA[5]    = 1.0 / (hz * hz);
3108:             col[6].i   = ex;
3109:             col[6].j   = ey;
3110:             col[6].k   = ez + 1;
3111:             col[6].loc = BACK;
3112:             col[6].c   = 0;
3113:             valA[6]    = 1.0 / (hz * hz);
3114:             col[7].i   = ex;
3115:             col[7].j   = ey;
3116:             col[7].k   = ez - 1;
3117:             col[7].loc = ELEMENT;
3118:             col[7].c   = 0;
3119:             valA[7]    = 1.0 / hz;
3120:             col[8].i   = ex;
3121:             col[8].j   = ey;
3122:             col[8].k   = ez;
3123:             col[8].loc = ELEMENT;
3124:             col[8].c   = 0;
3125:             valA[8]    = -1.0 / hz;
3126:           }
3127:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, nEntries, col, computed));
3128:           PetscCall(check_vals(ex, ey, ez, nEntries, valA, computed));
3129:         }

3131:         /* P equation : u_x + v_y + w_z = g
3132:            Note that this includes an explicit zero on the diagonal. This is only needed for
3133:            direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
3134:         {
3135:           DMStagStencil row, col[7];
3136:           PetscScalar   valA[7];

3138:           row.i      = ex;
3139:           row.j      = ey;
3140:           row.k      = ez;
3141:           row.loc    = ELEMENT;
3142:           row.c      = 0;
3143:           col[0].i   = ex;
3144:           col[0].j   = ey;
3145:           col[0].k   = ez;
3146:           col[0].loc = LEFT;
3147:           col[0].c   = 0;
3148:           valA[0]    = -1.0 / hx;
3149:           col[1].i   = ex;
3150:           col[1].j   = ey;
3151:           col[1].k   = ez;
3152:           col[1].loc = RIGHT;
3153:           col[1].c   = 0;
3154:           valA[1]    = 1.0 / hx;
3155:           col[2].i   = ex;
3156:           col[2].j   = ey;
3157:           col[2].k   = ez;
3158:           col[2].loc = DOWN;
3159:           col[2].c   = 0;
3160:           valA[2]    = -1.0 / hy;
3161:           col[3].i   = ex;
3162:           col[3].j   = ey;
3163:           col[3].k   = ez;
3164:           col[3].loc = UP;
3165:           col[3].c   = 0;
3166:           valA[3]    = 1.0 / hy;
3167:           col[4].i   = ex;
3168:           col[4].j   = ey;
3169:           col[4].k   = ez;
3170:           col[4].loc = BACK;
3171:           col[4].c   = 0;
3172:           valA[4]    = -1.0 / hz;
3173:           col[5].i   = ex;
3174:           col[5].j   = ey;
3175:           col[5].k   = ez;
3176:           col[5].loc = FRONT;
3177:           col[5].c   = 0;
3178:           valA[5]    = 1.0 / hz;
3179:           col[6]     = row;
3180:           valA[6]    = 0.0;
3181:           PetscCall(DMStagMatGetValuesStencil(dmSol, A, 1, &row, 7, col, computed));
3182:           PetscCall(check_vals(ex, ey, ez, 7, valA, computed));
3183:         }
3184:       }
3185:     }
3186:   }
3187:   PetscCall(DMStagVecRestoreArrayRead(dmCoord, coordLocal, &arrCoord));
3188:   PetscFunctionReturn(PETSC_SUCCESS);
3189: }

3191: /*TEST

3193:    test:
3194:       suffix: 1
3195:       nsize: 1

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

3201: TEST*/