1 static char help[] = "Test DMStagMatGetValuesStencil() in 3D\n\n";
2
3 #include <petscdm.h>
4 #include <petscksp.h>
5 #include <petscdmstag.h>
6
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
35
36 static PetscErrorCode CreateMat(DM, Mat *);
37 static PetscErrorCode CheckMat(DM, Mat);
38
main(int argc,char ** argv)39 int main(int argc, char **argv)
40 {
41 DM dmSol;
42 Mat A;
43
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 }
61
CreateMat(DM dmSol,Mat * pA)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;
72
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 }
96
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 }
133
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;
150
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 }
607
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;
624
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 }
1081
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;
1098
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 }
1555
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];
1562
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 }
1616
1617 /* A helper function to check values */
check_vals(PetscInt ex,PetscInt ey,PetscInt ez,PetscInt n,const PetscScalar * ref,const PetscScalar * computed)1618 static PetscErrorCode check_vals(PetscInt ex, PetscInt ey, PetscInt ez, PetscInt n, const PetscScalar *ref, const PetscScalar *computed)
1619 {
1620 PetscInt i;
1621
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 }
1628
1629 /* The same function as above, but getting and checking values, instead of setting them */
CheckMat(DM dmSol,Mat A)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];
1640
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 }
1662
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 }
1702
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;
1720
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 }
2178
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;
2196
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 }
2654
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;
2672
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 }
3130
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];
3137
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 }
3190
3191 /*TEST
3192
3193 test:
3194 suffix: 1
3195 nsize: 1
3196 output_file: output/empty.out
3197
3198 test:
3199 suffix: 2
3200 nsize: 4
3201 output_file: output/empty.out
3202
3203 TEST*/
3204