Actual source code: stagmulti.c

  1: /* Internal and DMStag-specific functions related to multigrid */
  2: #include <petsc/private/dmstagimpl.h>

  4: /*@C
  5:     DMStagRestrictSimple - restricts data from a fine to a coarse `DMSTAG`, in the simplest way

  7:     Values on coarse cells are averages of all fine cells that they cover.
  8:     Thus, values on vertices are injected, values on edges are averages
  9:     of the underlying two fine edges, and and values on elements in
 10:     d dimensions are averages of $2^d$ underlying elements.

 12:     Input Parameters:
 13: +   dmf - fine `DM`
 14: .   xf - data on fine `DM`
 15: -   dmc - coarse `DM`

 17:     Output Parameter:
 18: .   xc - data on coarse `DM`

 20:     Level: advanced

 22: .seealso: [](chapter_stag), `DMSTAG`, `DM`, `DMRestrict()`, `DMCoarsen()`, `DMSTAG`, `DMCreateInjection()`
 23: @*/
 24: PetscErrorCode DMStagRestrictSimple(DM dmf, Vec xf, DM dmc, Vec xc)
 25: {
 26:   PetscInt dim;

 28:   DMGetDimension(dmf, &dim);
 29:   switch (dim) {
 30:   case 1:
 31:     DMStagRestrictSimple_1d(dmf, xf, dmc, xc);
 32:     break;
 33:   case 2:
 34:     DMStagRestrictSimple_2d(dmf, xf, dmc, xc);
 35:     break;
 36:   case 3:
 37:     DMStagRestrictSimple_3d(dmf, xf, dmc, xc);
 38:     break;
 39:   default:
 40:     SETERRQ(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT "", dim);
 41:     break;
 42:   }
 43:   return 0;
 44: }

 46: /* Code duplication note: the next two functions are nearly identical, save the inclusion of the element terms */
 47: PETSC_INTERN PetscErrorCode DMStagPopulateInterpolation1d_a_b_Private(DM dmc, DM dmf, Mat A)
 48: {
 49:   PetscInt       exf, startexf, nexf, nextraxf, startexc;
 50:   PetscInt       dof[2];
 51:   const PetscInt dim = 1;

 53:   DMStagGetDOF(dmc, &dof[0], &dof[1], NULL, NULL);

 57:   /* In 1D, each fine point can receive data from at most 2 coarse points, at most one of which could be off-process */
 58:   MatSeqAIJSetPreallocation(A, 2, NULL);
 59:   MatMPIAIJSetPreallocation(A, 2, NULL, 1, NULL);

 61:   DMStagGetCorners(dmf, &startexf, NULL, NULL, &nexf, NULL, NULL, &nextraxf, NULL, NULL);
 62:   DMStagGetCorners(dmc, &startexc, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
 63:   for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
 64:     PetscInt exc, exf_local;
 65:     exf_local = exf - startexf;
 66:     exc       = startexc + exf_local / 2;
 67:     /* "even" vertices are just injected, odd vertices averaged */
 68:     if (exf_local % 2 == 0) {
 69:       DMStagStencil     rowf, colc;
 70:       PetscInt          ir, ic;
 71:       const PetscScalar one = 1.0;

 73:       rowf.i   = exf;
 74:       rowf.c   = 0;
 75:       rowf.loc = DMSTAG_LEFT;
 76:       colc.i   = exc;
 77:       colc.c   = 0;
 78:       colc.loc = DMSTAG_LEFT;
 79:       DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir);
 80:       DMStagStencilToIndexLocal(dmc, dim, 1, &colc, &ic);
 81:       MatSetValuesLocal(A, 1, &ir, 1, &ic, &one, INSERT_VALUES);
 82:     } else {
 83:       DMStagStencil     rowf, colc[2];
 84:       PetscInt          ir, ic[2];
 85:       const PetscScalar weight[2] = {0.5, 0.5};

 87:       rowf.i      = exf;
 88:       rowf.c      = 0;
 89:       rowf.loc    = DMSTAG_LEFT;
 90:       colc[0].i   = exc;
 91:       colc[0].c   = 0;
 92:       colc[0].loc = DMSTAG_LEFT;
 93:       colc[1].i   = exc;
 94:       colc[1].c   = 0;
 95:       colc[1].loc = DMSTAG_RIGHT;
 96:       DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir);
 97:       DMStagStencilToIndexLocal(dmc, dim, 2, colc, ic);
 98:       MatSetValuesLocal(A, 1, &ir, 2, ic, weight, INSERT_VALUES);
 99:     }
100:     /* Elements (excluding "extra" dummies) */
101:     if (dof[1] > 0 && exf < startexf + nexf) {
102:       DMStagStencil     rowf, colc;
103:       PetscInt          ir, ic;
104:       const PetscScalar weight = 1.0;

106:       rowf.i   = exf;
107:       rowf.c   = 0;
108:       rowf.loc = DMSTAG_ELEMENT; /* Note that this assumes only 1 dof */
109:       colc.i   = exc;
110:       colc.c   = 0;
111:       colc.loc = DMSTAG_ELEMENT;
112:       DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir);
113:       DMStagStencilToIndexLocal(dmc, dim, 1, &colc, &ic);
114:       MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES);
115:     }
116:   }
117:   return 0;
118: }

120: PETSC_INTERN PetscErrorCode DMStagPopulateInterpolation2d_0_a_b_Private(DM dmc, DM dmf, Mat A)
121: {
122:   PetscInt       exf, eyf, startexf, starteyf, nexf, neyf, nextraxf, nextrayf, startexc, starteyc, Nexf, Neyf;
123:   PetscInt       dof[3];
124:   const PetscInt dim = 2;

126:   DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], NULL);

130:   /* In 2D, each fine point can receive data from at most 4 coarse points, at most 3 of which could be off-process */
131:   MatSeqAIJSetPreallocation(A, 4, NULL);
132:   MatMPIAIJSetPreallocation(A, 4, NULL, 3, NULL);

134:   DMStagGetCorners(dmf, &startexf, &starteyf, NULL, &nexf, &neyf, NULL, &nextraxf, &nextrayf, NULL);
135:   DMStagGetCorners(dmc, &startexc, &starteyc, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
136:   DMStagGetGlobalSizes(dmf, &Nexf, &Neyf, NULL);
137:   for (eyf = starteyf; eyf < starteyf + neyf + nextrayf; ++eyf) {
138:     PetscInt eyc, eyf_local;

140:     eyf_local = eyf - starteyf;
141:     eyc       = starteyc + eyf_local / 2;
142:     for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
143:       PetscInt exc, exf_local;

145:       exf_local = exf - startexf;
146:       exc       = startexc + exf_local / 2;
147:       /* Left edges (excluding top "extra" dummy row) */
148:       if (eyf < starteyf + neyf) {
149:         DMStagStencil rowf, colc[4];
150:         PetscInt      ir, ic[4], nweight;
151:         PetscScalar   weight[4];

153:         rowf.i      = exf;
154:         rowf.j      = eyf;
155:         rowf.c      = 0;
156:         rowf.loc    = DMSTAG_LEFT;
157:         colc[0].i   = exc;
158:         colc[0].j   = eyc;
159:         colc[0].c   = 0;
160:         colc[0].loc = DMSTAG_LEFT;
161:         if (exf_local % 2 == 0) {
162:           if (eyf == Neyf - 1 || eyf == 0) {
163:             /* Note - this presumes something like a Neumann condition, assuming
164:                a ghost edge with the same value as the adjacent physical edge*/
165:             nweight   = 1;
166:             weight[0] = 1.0;
167:           } else {
168:             nweight   = 2;
169:             weight[0] = 0.75;
170:             weight[1] = 0.25;
171:             if (eyf_local % 2 == 0) {
172:               colc[1].i   = exc;
173:               colc[1].j   = eyc - 1;
174:               colc[1].c   = 0;
175:               colc[1].loc = DMSTAG_LEFT;
176:             } else {
177:               colc[1].i   = exc;
178:               colc[1].j   = eyc + 1;
179:               colc[1].c   = 0;
180:               colc[1].loc = DMSTAG_LEFT;
181:             }
182:           }
183:         } else {
184:           colc[1].i   = exc;
185:           colc[1].j   = eyc;
186:           colc[1].c   = 0;
187:           colc[1].loc = DMSTAG_RIGHT;
188:           if (eyf == Neyf - 1 || eyf == 0) {
189:             /* Note - this presumes something like a Neumann condition, assuming
190:                a ghost edge with the same value as the adjacent physical edge*/
191:             nweight   = 2;
192:             weight[0] = 0.5;
193:             weight[1] = 0.5;
194:           } else {
195:             nweight   = 4;
196:             weight[0] = 0.375;
197:             weight[1] = 0.375;
198:             weight[2] = 0.125;
199:             weight[3] = 0.125;
200:             if (eyf_local % 2 == 0) {
201:               colc[2].i   = exc;
202:               colc[2].j   = eyc - 1;
203:               colc[2].c   = 0;
204:               colc[2].loc = DMSTAG_LEFT;
205:               colc[3].i   = exc;
206:               colc[3].j   = eyc - 1;
207:               colc[3].c   = 0;
208:               colc[3].loc = DMSTAG_RIGHT;
209:             } else {
210:               colc[2].i   = exc;
211:               colc[2].j   = eyc + 1;
212:               colc[2].c   = 0;
213:               colc[2].loc = DMSTAG_LEFT;
214:               colc[3].i   = exc;
215:               colc[3].j   = eyc + 1;
216:               colc[3].c   = 0;
217:               colc[3].loc = DMSTAG_RIGHT;
218:             }
219:           }
220:         }
221:         DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir);
222:         DMStagStencilToIndexLocal(dmc, dim, nweight, colc, ic);
223:         MatSetValuesLocal(A, 1, &ir, nweight, ic, weight, INSERT_VALUES);
224:       }
225:       /* Down edges (excluding right "extra" dummy col) */
226:       if (exf < startexf + nexf) {
227:         DMStagStencil rowf, colc[4];
228:         PetscInt      ir, ic[4], nweight;
229:         PetscScalar   weight[4];

231:         rowf.i      = exf;
232:         rowf.j      = eyf;
233:         rowf.c      = 0;
234:         rowf.loc    = DMSTAG_DOWN;
235:         colc[0].i   = exc;
236:         colc[0].j   = eyc;
237:         colc[0].c   = 0;
238:         colc[0].loc = DMSTAG_DOWN;
239:         if (eyf_local % 2 == 0) {
240:           if (exf == Nexf - 1 || exf == 0) {
241:             /* Note - this presumes something like a Neumann condition, assuming
242:                a ghost edge with the same value as the adjacent physical edge*/
243:             nweight   = 1;
244:             weight[0] = 1.0;
245:           } else {
246:             nweight   = 2;
247:             weight[0] = 0.75;
248:             weight[1] = 0.25;
249:             if (exf_local % 2 == 0) {
250:               colc[1].i   = exc - 1;
251:               colc[1].j   = eyc;
252:               colc[1].c   = 0;
253:               colc[1].loc = DMSTAG_DOWN;
254:             } else {
255:               colc[1].i   = exc + 1;
256:               colc[1].j   = eyc;
257:               colc[1].c   = 0;
258:               colc[1].loc = DMSTAG_DOWN;
259:             }
260:           }
261:         } else {
262:           colc[1].i   = exc;
263:           colc[1].j   = eyc;
264:           colc[1].c   = 0;
265:           colc[1].loc = DMSTAG_UP;
266:           if (exf == Nexf - 1 || exf == 0) {
267:             /* Note - this presumes something like a Neumann condition, assuming
268:                a ghost edge with the same value as the adjacent physical edge*/
269:             nweight   = 2;
270:             weight[0] = 0.5;
271:             weight[1] = 0.5;
272:           } else {
273:             nweight   = 4;
274:             weight[0] = 0.375;
275:             weight[1] = 0.375;
276:             weight[2] = 0.125;
277:             weight[3] = 0.125;
278:             if (exf_local % 2 == 0) {
279:               colc[2].i   = exc - 1;
280:               colc[2].j   = eyc;
281:               colc[2].c   = 0;
282:               colc[2].loc = DMSTAG_DOWN;
283:               colc[3].i   = exc - 1;
284:               colc[3].j   = eyc;
285:               colc[3].c   = 0;
286:               colc[3].loc = DMSTAG_UP;
287:             } else {
288:               colc[2].i   = exc + 1;
289:               colc[2].j   = eyc;
290:               colc[2].c   = 0;
291:               colc[2].loc = DMSTAG_DOWN;
292:               colc[3].i   = exc + 1;
293:               colc[3].j   = eyc;
294:               colc[3].c   = 0;
295:               colc[3].loc = DMSTAG_UP;
296:             }
297:           }
298:         }
299:         DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir);
300:         DMStagStencilToIndexLocal(dmc, dim, nweight, colc, ic);
301:         MatSetValuesLocal(A, 1, &ir, nweight, ic, weight, INSERT_VALUES);
302:       }
303:       /* Elements (excluding "extra" dummy) */
304:       if (dof[2] > 0 && exf < startexf + nexf && eyf < starteyf + neyf) {
305:         DMStagStencil     rowf, colc;
306:         PetscInt          ir, ic;
307:         const PetscScalar weight = 1.0;

309:         rowf.i   = exf;
310:         rowf.j   = eyf;
311:         rowf.c   = 0;
312:         rowf.loc = DMSTAG_ELEMENT;
313:         colc.i   = exc;
314:         colc.j   = eyc;
315:         colc.c   = 0;
316:         colc.loc = DMSTAG_ELEMENT;
317:         DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir);
318:         DMStagStencilToIndexLocal(dmc, dim, 1, &colc, &ic);
319:         MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES);
320:       }
321:     }
322:   }
323:   return 0;
324: }

326: PETSC_INTERN PetscErrorCode DMStagPopulateInterpolation3d_0_0_a_b_Private(DM dmc, DM dmf, Mat A)
327: {
328:   PetscInt       exf, eyf, ezf, startexf, starteyf, startezf, nexf, neyf, nezf, nextraxf, nextrayf, nextrazf, startexc, starteyc, startezc, Nexf, Neyf, Nezf;
329:   PetscInt       dof[4];
330:   const PetscInt dim = 3;


333:   DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], &dof[3]);

337:   /* In 3D, each fine point can receive data from at most 8 coarse points, at most 7 of which could be off-process */
338:   MatSeqAIJSetPreallocation(A, 8, NULL);
339:   MatMPIAIJSetPreallocation(A, 8, NULL, 7, NULL);

341:   DMStagGetCorners(dmf, &startexf, &starteyf, &startezf, &nexf, &neyf, &nezf, &nextraxf, &nextrayf, &nextrazf);
342:   DMStagGetCorners(dmc, &startexc, &starteyc, &startezc, NULL, NULL, NULL, NULL, NULL, NULL);
343:   DMStagGetGlobalSizes(dmf, &Nexf, &Neyf, &Nezf);
344:   for (ezf = startezf; ezf < startezf + nezf + nextrayf; ++ezf) {
345:     const PetscInt  ezf_local  = ezf - startezf;
346:     const PetscInt  ezc        = startezc + ezf_local / 2;
347:     const PetscBool boundary_z = (PetscBool)(ezf == 0 || ezf == Nezf - 1);

349:     for (eyf = starteyf; eyf < starteyf + neyf + nextrayf; ++eyf) {
350:       const PetscInt  eyf_local  = eyf - starteyf;
351:       const PetscInt  eyc        = starteyc + eyf_local / 2;
352:       const PetscBool boundary_y = (PetscBool)(eyf == 0 || eyf == Neyf - 1);

354:       for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
355:         const PetscInt  exf_local  = exf - startexf;
356:         const PetscInt  exc        = startexc + exf_local / 2;
357:         const PetscBool boundary_x = (PetscBool)(exf == 0 || exf == Nexf - 1);

359:         /* Left faces (excluding top and front "extra" dummy layers) */
360:         if (eyf < starteyf + neyf && ezf < startezf + nezf) {
361:           DMStagStencil rowf, colc[8];
362:           PetscInt      ir, ic[8], nweight;
363:           PetscScalar   weight[8];

365:           rowf.i      = exf;
366:           rowf.j      = eyf;
367:           rowf.k      = ezf;
368:           rowf.c      = 0;
369:           rowf.loc    = DMSTAG_LEFT;
370:           colc[0].i   = exc;
371:           colc[0].j   = eyc;
372:           colc[0].k   = ezc;
373:           colc[0].c   = 0;
374:           colc[0].loc = DMSTAG_LEFT;
375:           if (exf_local % 2 == 0) {
376:             if (boundary_y) {
377:               if (boundary_z) {
378:                 nweight   = 1;
379:                 weight[0] = 1.0;
380:               } else {
381:                 const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;

383:                 nweight     = 2;
384:                 weight[0]   = 0.75;
385:                 weight[1]   = 0.25;
386:                 colc[1].i   = exc;
387:                 colc[1].j   = eyc;
388:                 colc[1].k   = ezc + ezc_offset;
389:                 colc[1].c   = 0;
390:                 colc[1].loc = DMSTAG_LEFT;
391:               }
392:             } else {
393:               const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;

395:               if (boundary_z) {
396:                 nweight     = 2;
397:                 weight[0]   = 0.75;
398:                 weight[1]   = 0.25;
399:                 colc[1].i   = exc;
400:                 colc[1].j   = eyc + eyc_offset;
401:                 colc[1].k   = ezc;
402:                 colc[1].c   = 0;
403:                 colc[1].loc = DMSTAG_LEFT;
404:               } else {
405:                 const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;

407:                 nweight   = 4;
408:                 weight[0] = 0.75 * 0.75;
409:                 weight[1] = weight[2] = 0.25 * 0.75;
410:                 weight[3]             = 0.25 * 0.25;
411:                 colc[1].i             = exc;
412:                 colc[1].j             = eyc;
413:                 colc[1].k             = ezc + ezc_offset;
414:                 colc[1].c             = 0;
415:                 colc[1].loc           = DMSTAG_LEFT;
416:                 colc[2].i             = exc;
417:                 colc[2].j             = eyc + eyc_offset;
418:                 colc[2].k             = ezc;
419:                 colc[2].c             = 0;
420:                 colc[2].loc           = DMSTAG_LEFT;
421:                 colc[3].i             = exc;
422:                 colc[3].j             = eyc + eyc_offset;
423:                 colc[3].k             = ezc + ezc_offset;
424:                 colc[3].c             = 0;
425:                 colc[3].loc           = DMSTAG_LEFT;
426:               }
427:             }
428:           } else {
429:             colc[1].i   = exc;
430:             colc[1].j   = eyc;
431:             colc[1].k   = ezc;
432:             colc[1].c   = 0;
433:             colc[1].loc = DMSTAG_RIGHT;
434:             if (boundary_y) {
435:               if (boundary_z) {
436:                 nweight   = 2;
437:                 weight[0] = weight[1] = 0.5;
438:               } else {
439:                 const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;

441:                 nweight   = 4;
442:                 weight[0] = weight[1] = 0.5 * 0.75;
443:                 weight[2] = weight[3] = 0.5 * 0.25;
444:                 colc[2].i             = exc;
445:                 colc[2].j             = eyc;
446:                 colc[2].k             = ezc + ezc_offset;
447:                 colc[2].c             = 0;
448:                 colc[2].loc           = DMSTAG_LEFT;
449:                 colc[3].i             = exc;
450:                 colc[3].j             = eyc;
451:                 colc[3].k             = ezc + ezc_offset;
452:                 colc[3].c             = 0;
453:                 colc[3].loc           = DMSTAG_RIGHT;
454:               }
455:             } else {
456:               const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;

458:               if (boundary_z) {
459:                 nweight   = 4;
460:                 weight[0] = weight[1] = 0.5 * 0.75;
461:                 weight[2] = weight[3] = 0.5 * 0.25;
462:                 colc[2].i             = exc;
463:                 colc[2].j             = eyc + eyc_offset;
464:                 colc[2].k             = ezc;
465:                 colc[2].c             = 0;
466:                 colc[2].loc           = DMSTAG_LEFT;
467:                 colc[3].i             = exc;
468:                 colc[3].j             = eyc + eyc_offset;
469:                 colc[3].k             = ezc;
470:                 colc[3].c             = 0;
471:                 colc[3].loc           = DMSTAG_RIGHT;
472:               } else {
473:                 const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;

475:                 nweight   = 8;
476:                 weight[0] = weight[1] = 0.5 * 0.75 * 0.75;
477:                 weight[2] = weight[3] = weight[4] = weight[5] = 0.5 * 0.25 * 0.75;
478:                 weight[6] = weight[7] = 0.5 * 0.25 * 0.25;
479:                 colc[2].i             = exc;
480:                 colc[2].j             = eyc;
481:                 colc[2].k             = ezc + ezc_offset;
482:                 colc[2].c             = 0;
483:                 colc[2].loc           = DMSTAG_LEFT;
484:                 colc[3].i             = exc;
485:                 colc[3].j             = eyc;
486:                 colc[3].k             = ezc + ezc_offset;
487:                 colc[3].c             = 0;
488:                 colc[3].loc           = DMSTAG_RIGHT;
489:                 colc[4].i             = exc;
490:                 colc[4].j             = eyc + eyc_offset;
491:                 colc[4].k             = ezc;
492:                 colc[4].c             = 0;
493:                 colc[4].loc           = DMSTAG_LEFT;
494:                 colc[5].i             = exc;
495:                 colc[5].j             = eyc + eyc_offset;
496:                 colc[5].k             = ezc;
497:                 colc[5].c             = 0;
498:                 colc[5].loc           = DMSTAG_RIGHT;
499:                 colc[6].i             = exc;
500:                 colc[6].j             = eyc + eyc_offset;
501:                 colc[6].k             = ezc + ezc_offset;
502:                 colc[6].c             = 0;
503:                 colc[6].loc           = DMSTAG_LEFT;
504:                 colc[7].i             = exc;
505:                 colc[7].j             = eyc + eyc_offset;
506:                 colc[7].k             = ezc + ezc_offset;
507:                 colc[7].c             = 0;
508:                 colc[7].loc           = DMSTAG_RIGHT;
509:               }
510:             }
511:           }
512:           DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir);
513:           DMStagStencilToIndexLocal(dmc, dim, nweight, colc, ic);
514:           MatSetValuesLocal(A, 1, &ir, nweight, ic, weight, INSERT_VALUES);
515:         }

517:         /* Bottom faces (excluding left and front "extra" dummy layers) */
518:         if (exf < startexf + nexf && ezf < startezf + nezf) {
519:           DMStagStencil rowf, colc[8];
520:           PetscInt      ir, ic[8], nweight;
521:           PetscScalar   weight[8];

523:           rowf.i      = exf;
524:           rowf.j      = eyf;
525:           rowf.k      = ezf;
526:           rowf.c      = 0;
527:           rowf.loc    = DMSTAG_DOWN;
528:           colc[0].i   = exc;
529:           colc[0].j   = eyc;
530:           colc[0].k   = ezc;
531:           colc[0].c   = 0;
532:           colc[0].loc = DMSTAG_DOWN;
533:           if (eyf_local % 2 == 0) {
534:             if (boundary_x) {
535:               if (boundary_z) {
536:                 nweight   = 1;
537:                 weight[0] = 1.0;
538:               } else {
539:                 const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;

541:                 nweight     = 2;
542:                 weight[0]   = 0.75;
543:                 weight[1]   = 0.25;
544:                 colc[1].i   = exc;
545:                 colc[1].j   = eyc;
546:                 colc[1].k   = ezc + ezc_offset;
547:                 colc[1].c   = 0;
548:                 colc[1].loc = DMSTAG_DOWN;
549:               }
550:             } else {
551:               const PetscInt exc_offset = exf_local % 2 == 0 ? -1 : 1;

553:               if (boundary_z) {
554:                 nweight     = 2;
555:                 weight[0]   = 0.75;
556:                 weight[1]   = 0.25;
557:                 colc[1].i   = exc + exc_offset;
558:                 colc[1].j   = eyc;
559:                 colc[1].k   = ezc;
560:                 colc[1].c   = 0;
561:                 colc[1].loc = DMSTAG_DOWN;
562:               } else {
563:                 const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;

565:                 nweight   = 4;
566:                 weight[0] = 0.75 * 0.75;
567:                 weight[1] = weight[2] = 0.25 * 0.75;
568:                 weight[3]             = 0.25 * 0.25;
569:                 colc[1].i             = exc;
570:                 colc[1].j             = eyc;
571:                 colc[1].k             = ezc + ezc_offset;
572:                 colc[1].c             = 0;
573:                 colc[1].loc           = DMSTAG_DOWN;
574:                 colc[2].i             = exc + exc_offset;
575:                 colc[2].j             = eyc;
576:                 colc[2].k             = ezc;
577:                 colc[2].c             = 0;
578:                 colc[2].loc           = DMSTAG_DOWN;
579:                 colc[3].i             = exc + exc_offset;
580:                 colc[3].j             = eyc;
581:                 colc[3].k             = ezc + ezc_offset;
582:                 colc[3].c             = 0;
583:                 colc[3].loc           = DMSTAG_DOWN;
584:               }
585:             }
586:           } else {
587:             colc[1].i   = exc;
588:             colc[1].j   = eyc;
589:             colc[1].k   = ezc;
590:             colc[1].c   = 0;
591:             colc[1].loc = DMSTAG_UP;
592:             if (boundary_x) {
593:               if (boundary_z) {
594:                 nweight   = 2;
595:                 weight[0] = weight[1] = 0.5;
596:               } else {
597:                 const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;

599:                 nweight   = 4;
600:                 weight[0] = weight[1] = 0.5 * 0.75;
601:                 weight[2] = weight[3] = 0.5 * 0.25;
602:                 colc[2].i             = exc;
603:                 colc[2].j             = eyc;
604:                 colc[2].k             = ezc + ezc_offset;
605:                 colc[2].c             = 0;
606:                 colc[2].loc           = DMSTAG_DOWN;
607:                 colc[3].i             = exc;
608:                 colc[3].j             = eyc;
609:                 colc[3].k             = ezc + ezc_offset;
610:                 colc[3].c             = 0;
611:                 colc[3].loc           = DMSTAG_UP;
612:               }
613:             } else {
614:               const PetscInt exc_offset = exf_local % 2 == 0 ? -1 : 1;

616:               if (boundary_z) {
617:                 nweight   = 4;
618:                 weight[0] = weight[1] = 0.5 * 0.75;
619:                 weight[2] = weight[3] = 0.5 * 0.25;
620:                 colc[2].i             = exc + exc_offset;
621:                 colc[2].j             = eyc;
622:                 colc[2].k             = ezc;
623:                 colc[2].c             = 0;
624:                 colc[2].loc           = DMSTAG_DOWN;
625:                 colc[3].i             = exc + exc_offset;
626:                 colc[3].j             = eyc;
627:                 colc[3].k             = ezc;
628:                 colc[3].c             = 0;
629:                 colc[3].loc           = DMSTAG_UP;
630:               } else {
631:                 const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;

633:                 nweight   = 8;
634:                 weight[0] = weight[1] = 0.5 * 0.75 * 0.75;
635:                 weight[2] = weight[3] = weight[4] = weight[5] = 0.5 * 0.25 * 0.75;
636:                 weight[6] = weight[7] = 0.5 * 0.25 * 0.25;
637:                 colc[2].i             = exc;
638:                 colc[2].j             = eyc;
639:                 colc[2].k             = ezc + ezc_offset;
640:                 colc[2].c             = 0;
641:                 colc[2].loc           = DMSTAG_DOWN;
642:                 colc[3].i             = exc;
643:                 colc[3].j             = eyc;
644:                 colc[3].k             = ezc + ezc_offset;
645:                 colc[3].c             = 0;
646:                 colc[3].loc           = DMSTAG_UP;
647:                 colc[4].i             = exc + exc_offset;
648:                 colc[4].j             = eyc;
649:                 colc[4].k             = ezc;
650:                 colc[4].c             = 0;
651:                 colc[4].loc           = DMSTAG_DOWN;
652:                 colc[5].i             = exc + exc_offset;
653:                 colc[5].j             = eyc;
654:                 colc[5].k             = ezc;
655:                 colc[5].c             = 0;
656:                 colc[5].loc           = DMSTAG_UP;
657:                 colc[6].i             = exc + exc_offset;
658:                 colc[6].j             = eyc;
659:                 colc[6].k             = ezc + ezc_offset;
660:                 colc[6].c             = 0;
661:                 colc[6].loc           = DMSTAG_DOWN;
662:                 colc[7].i             = exc + exc_offset;
663:                 colc[7].j             = eyc;
664:                 colc[7].k             = ezc + ezc_offset;
665:                 colc[7].c             = 0;
666:                 colc[7].loc           = DMSTAG_UP;
667:               }
668:             }
669:           }
670:           DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir);
671:           DMStagStencilToIndexLocal(dmc, dim, nweight, colc, ic);
672:           MatSetValuesLocal(A, 1, &ir, nweight, ic, weight, INSERT_VALUES);
673:         }

675:         /* Back faces (excluding left and top "extra" dummy layers) */
676:         if (exf < startexf + nexf && ezf < startezf + nezf) {
677:           DMStagStencil rowf, colc[8];
678:           PetscInt      ir, ic[8], nweight;
679:           PetscScalar   weight[8];

681:           rowf.i      = exf;
682:           rowf.j      = eyf;
683:           rowf.k      = ezf;
684:           rowf.c      = 0;
685:           rowf.loc    = DMSTAG_BACK;
686:           colc[0].i   = exc;
687:           colc[0].j   = eyc;
688:           colc[0].k   = ezc;
689:           colc[0].c   = 0;
690:           colc[0].loc = DMSTAG_BACK;
691:           if (ezf_local % 2 == 0) {
692:             if (boundary_x) {
693:               if (boundary_y) {
694:                 nweight   = 1;
695:                 weight[0] = 1.0;
696:               } else {
697:                 const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;

699:                 nweight     = 2;
700:                 weight[0]   = 0.75;
701:                 weight[1]   = 0.25;
702:                 colc[1].i   = exc;
703:                 colc[1].j   = eyc + eyc_offset;
704:                 colc[1].k   = ezc;
705:                 colc[1].c   = 0;
706:                 colc[1].loc = DMSTAG_BACK;
707:               }
708:             } else {
709:               const PetscInt exc_offset = exf_local % 2 == 0 ? -1 : 1;

711:               if (boundary_y) {
712:                 nweight     = 2;
713:                 weight[0]   = 0.75;
714:                 weight[1]   = 0.25;
715:                 colc[1].i   = exc + exc_offset;
716:                 colc[1].j   = eyc;
717:                 colc[1].k   = ezc;
718:                 colc[1].c   = 0;
719:                 colc[1].loc = DMSTAG_BACK;
720:               } else {
721:                 const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;

723:                 nweight   = 4;
724:                 weight[0] = 0.75 * 0.75;
725:                 weight[1] = weight[2] = 0.25 * 0.75;
726:                 weight[3]             = 0.25 * 0.25;
727:                 colc[1].i             = exc + exc_offset;
728:                 colc[1].j             = eyc;
729:                 colc[1].k             = ezc;
730:                 colc[1].c             = 0;
731:                 colc[1].loc           = DMSTAG_BACK;
732:                 colc[2].i             = exc;
733:                 colc[2].j             = eyc + eyc_offset;
734:                 colc[2].k             = ezc;
735:                 colc[2].c             = 0;
736:                 colc[2].loc           = DMSTAG_BACK;
737:                 colc[3].i             = exc + exc_offset;
738:                 colc[3].j             = eyc + eyc_offset;
739:                 colc[3].k             = ezc;
740:                 colc[3].c             = 0;
741:                 colc[3].loc           = DMSTAG_BACK;
742:               }
743:             }
744:           } else {
745:             colc[1].i   = exc;
746:             colc[1].j   = eyc;
747:             colc[1].k   = ezc;
748:             colc[1].c   = 0;
749:             colc[1].loc = DMSTAG_FRONT;
750:             if (boundary_x) {
751:               if (boundary_y) {
752:                 nweight   = 2;
753:                 weight[0] = weight[1] = 0.5;
754:               } else {
755:                 const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;

757:                 nweight   = 4;
758:                 weight[0] = weight[1] = 0.5 * 0.75;
759:                 weight[2] = weight[3] = 0.5 * 0.25;
760:                 colc[2].i             = exc;
761:                 colc[2].j             = eyc + eyc_offset;
762:                 colc[2].k             = ezc;
763:                 colc[2].c             = 0;
764:                 colc[2].loc           = DMSTAG_BACK;
765:                 colc[3].i             = exc;
766:                 colc[3].j             = eyc + eyc_offset;
767:                 colc[3].k             = ezc;
768:                 colc[3].c             = 0;
769:                 colc[3].loc           = DMSTAG_FRONT;
770:               }
771:             } else {
772:               const PetscInt exc_offset = exf_local % 2 == 0 ? -1 : 1;

774:               if (boundary_y) {
775:                 nweight   = 4;
776:                 weight[0] = weight[1] = 0.5 * 0.75;
777:                 weight[2] = weight[3] = 0.5 * 0.25;
778:                 colc[2].i             = exc + exc_offset;
779:                 colc[2].j             = eyc;
780:                 colc[2].k             = ezc;
781:                 colc[2].c             = 0;
782:                 colc[2].loc           = DMSTAG_BACK;
783:                 colc[3].i             = exc + exc_offset;
784:                 colc[3].j             = eyc;
785:                 colc[3].k             = ezc;
786:                 colc[3].c             = 0;
787:                 colc[3].loc           = DMSTAG_FRONT;
788:               } else {
789:                 const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;

791:                 nweight   = 8;
792:                 weight[0] = weight[1] = 0.5 * 0.75 * 0.75;
793:                 weight[2] = weight[3] = weight[4] = weight[5] = 0.5 * 0.25 * 0.75;
794:                 weight[6] = weight[7] = 0.5 * 0.25 * 0.25;
795:                 colc[2].i             = exc;
796:                 colc[2].j             = eyc + eyc_offset;
797:                 colc[2].k             = ezc;
798:                 colc[2].c             = 0;
799:                 colc[2].loc           = DMSTAG_BACK;
800:                 colc[3].i             = exc;
801:                 colc[3].j             = eyc + eyc_offset;
802:                 colc[3].k             = ezc;
803:                 colc[3].c             = 0;
804:                 colc[3].loc           = DMSTAG_FRONT;
805:                 colc[4].i             = exc + exc_offset;
806:                 colc[4].j             = eyc;
807:                 colc[4].k             = ezc;
808:                 colc[4].c             = 0;
809:                 colc[4].loc           = DMSTAG_BACK;
810:                 colc[5].i             = exc + exc_offset;
811:                 colc[5].j             = eyc;
812:                 colc[5].k             = ezc;
813:                 colc[5].c             = 0;
814:                 colc[5].loc           = DMSTAG_FRONT;
815:                 colc[6].i             = exc + exc_offset;
816:                 colc[6].j             = eyc + eyc_offset;
817:                 colc[6].k             = ezc;
818:                 colc[6].c             = 0;
819:                 colc[6].loc           = DMSTAG_BACK;
820:                 colc[7].i             = exc + exc_offset;
821:                 colc[7].j             = eyc + eyc_offset;
822:                 colc[7].k             = ezc;
823:                 colc[7].c             = 0;
824:                 colc[7].loc           = DMSTAG_FRONT;
825:               }
826:             }
827:           }
828:           DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir);
829:           DMStagStencilToIndexLocal(dmc, dim, nweight, colc, ic);
830:           MatSetValuesLocal(A, 1, &ir, nweight, ic, weight, INSERT_VALUES);
831:         }
832:         /* Elements */
833:         if (dof[3] == 1 && exf < startexf + nexf && eyf < starteyf + neyf && ezf < startezf + nezf) {
834:           DMStagStencil     rowf, colc;
835:           PetscInt          ir, ic;
836:           const PetscScalar weight = 1.0;

838:           rowf.i   = exf;
839:           rowf.j   = eyf;
840:           rowf.k   = ezf;
841:           rowf.c   = 0;
842:           rowf.loc = DMSTAG_ELEMENT;
843:           colc.i   = exc;
844:           colc.j   = eyc;
845:           colc.k   = ezc;
846:           colc.c   = 0;
847:           colc.loc = DMSTAG_ELEMENT;
848:           DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir);
849:           DMStagStencilToIndexLocal(dmc, dim, 1, &colc, &ic);
850:           MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES);
851:         }
852:       }
853:     }
854:   }
855:   return 0;
856: }

858: PETSC_INTERN PetscErrorCode DMStagPopulateRestriction1d_a_b_Private(DM dmc, DM dmf, Mat A)
859: {
860:   PetscInt       exf, startexf, nexf, nextraxf, startexc, Nexf;
861:   PetscInt       dof[2];
862:   const PetscInt dim = 1;

864:   DMStagGetDOF(dmc, &dof[0], &dof[1], NULL, NULL);

868:   /* In 1D, each coarse point can receive from up to 3 fine points, one of which may be off-rank */
869:   MatSeqAIJSetPreallocation(A, 3, NULL);
870:   MatMPIAIJSetPreallocation(A, 3, NULL, 1, NULL);

872:   DMStagGetCorners(dmf, &startexf, NULL, NULL, &nexf, NULL, NULL, &nextraxf, NULL, NULL);
873:   DMStagGetCorners(dmc, &startexc, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
874:   DMStagGetGlobalSizes(dmf, &Nexf, NULL, NULL);
875:   for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
876:     PetscInt exc, exf_local;

878:     exf_local = exf - startexf;
879:     exc       = startexc + exf_local / 2;
880:     /* "even" vertices contribute to the overlying coarse vertex, odd vertices to the two adjacent */
881:     if (exf_local % 2 == 0) {
882:       DMStagStencil colf, rowc;
883:       PetscInt      ir, ic;
884:       PetscScalar   weight;

886:       colf.i   = exf;
887:       colf.c   = 0;
888:       colf.loc = DMSTAG_LEFT;
889:       rowc.i   = exc;
890:       rowc.c   = 0;
891:       rowc.loc = DMSTAG_LEFT;
892:       DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir);
893:       DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic);
894:       weight = (exf == Nexf || exf == 0) ? 0.75 : 0.5; /* Assume a Neuman-type condition */
895:       MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES);
896:     } else {
897:       DMStagStencil     colf, rowc[2];
898:       PetscInt          ic, ir[2];
899:       const PetscScalar weight[2] = {0.25, 0.25};

901:       colf.i      = exf;
902:       colf.c      = 0;
903:       colf.loc    = DMSTAG_LEFT;
904:       rowc[0].i   = exc;
905:       rowc[0].c   = 0;
906:       rowc[0].loc = DMSTAG_LEFT;
907:       rowc[1].i   = exc;
908:       rowc[1].c   = 0;
909:       rowc[1].loc = DMSTAG_RIGHT;
910:       DMStagStencilToIndexLocal(dmc, dim, 2, rowc, ir);
911:       DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic);
912:       MatSetValuesLocal(A, 2, ir, 1, &ic, weight, INSERT_VALUES);
913:     }
914:     if (dof[1] > 0 && exf < startexf + nexf) {
915:       DMStagStencil     rowc, colf;
916:       PetscInt          ir, ic;
917:       const PetscScalar weight = 0.5;

919:       rowc.i   = exc;
920:       rowc.c   = 0;
921:       rowc.loc = DMSTAG_ELEMENT;
922:       colf.i   = exf;
923:       colf.c   = 0;
924:       colf.loc = DMSTAG_ELEMENT;
925:       DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir);
926:       DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic);
927:       MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES);
928:     }
929:   }
930:   return 0;
931: }

933: PETSC_INTERN PetscErrorCode DMStagPopulateRestriction2d_0_a_b_Private(DM dmc, DM dmf, Mat A)
934: {
935:   PetscInt       exf, eyf, startexf, starteyf, nexf, neyf, nextraxf, nextrayf, startexc, starteyc, Nexf, Neyf;
936:   PetscInt       dof[3];
937:   const PetscInt dim = 2;

939:   DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], NULL);

943:   /* In 2D, each coarse point can receive from up to 6 fine points,
944:      up to 2 of which may be off rank */
945:   MatSeqAIJSetPreallocation(A, 6, NULL);
946:   MatMPIAIJSetPreallocation(A, 6, NULL, 2, NULL);

948:   DMStagGetCorners(dmf, &startexf, &starteyf, NULL, &nexf, &neyf, NULL, &nextraxf, &nextrayf, NULL);
949:   DMStagGetCorners(dmc, &startexc, &starteyc, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
950:   DMStagGetGlobalSizes(dmf, &Nexf, &Neyf, NULL);
951:   for (eyf = starteyf; eyf < starteyf + neyf + nextrayf; ++eyf) {
952:     PetscInt eyc, eyf_local;
953:     eyf_local = eyf - starteyf;
954:     eyc       = starteyc + eyf_local / 2;
955:     for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
956:       PetscInt exc, exf_local;
957:       exf_local = exf - startexf;
958:       exc       = startexc + exf_local / 2;
959:       /* Left edges (excluding top "extra" dummy row) */
960:       if (eyf < starteyf + neyf) {
961:         DMStagStencil rowc[2], colf;
962:         PetscInt      ir[2], ic, nweight;
963:         PetscScalar   weight[2];

965:         colf.i      = exf;
966:         colf.j      = eyf;
967:         colf.c      = 0;
968:         colf.loc    = DMSTAG_LEFT;
969:         rowc[0].i   = exc;
970:         rowc[0].j   = eyc;
971:         rowc[0].c   = 0;
972:         rowc[0].loc = DMSTAG_LEFT;
973:         if (exf_local % 2 == 0) {
974:           nweight = 1;
975:           if (exf == Nexf || exf == 0) {
976:             /* Note - this presumes something like a Neumann condition, assuming
977:                a ghost edge with the same value as the adjacent physical edge*/
978:             weight[0] = 0.375;
979:           } else {
980:             weight[0] = 0.25;
981:           }
982:         } else {
983:           nweight     = 2;
984:           rowc[1].i   = exc;
985:           rowc[1].j   = eyc;
986:           rowc[1].c   = 0;
987:           rowc[1].loc = DMSTAG_RIGHT;
988:           weight[0]   = 0.125;
989:           weight[1]   = 0.125;
990:         }
991:         DMStagStencilToIndexLocal(dmc, dim, nweight, rowc, ir);
992:         DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic);
993:         MatSetValuesLocal(A, nweight, ir, 1, &ic, weight, INSERT_VALUES);
994:       }
995:       /* Down edges (excluding right "extra" dummy col) */
996:       if (exf < startexf + nexf) {
997:         DMStagStencil rowc[2], colf;
998:         PetscInt      ir[2], ic, nweight;
999:         PetscScalar   weight[2];

1001:         colf.i      = exf;
1002:         colf.j      = eyf;
1003:         colf.c      = 0;
1004:         colf.loc    = DMSTAG_DOWN;
1005:         rowc[0].i   = exc;
1006:         rowc[0].j   = eyc;
1007:         rowc[0].c   = 0;
1008:         rowc[0].loc = DMSTAG_DOWN;
1009:         if (eyf_local % 2 == 0) {
1010:           nweight = 1;
1011:           if (eyf == Neyf || eyf == 0) {
1012:             /* Note - this presumes something like a Neumann condition, assuming
1013:                a ghost edge with the same value as the adjacent physical edge*/
1014:             weight[0] = 0.375;
1015:           } else {
1016:             weight[0] = 0.25;
1017:           }
1018:         } else {
1019:           nweight     = 2;
1020:           rowc[1].i   = exc;
1021:           rowc[1].j   = eyc;
1022:           rowc[1].c   = 0;
1023:           rowc[1].loc = DMSTAG_UP;
1024:           weight[0]   = 0.125;
1025:           weight[1]   = 0.125;
1026:         }
1027:         DMStagStencilToIndexLocal(dmc, dim, nweight, rowc, ir);
1028:         DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic);
1029:         MatSetValuesLocal(A, nweight, ir, 1, &ic, weight, INSERT_VALUES);
1030:       }
1031:       /* Elements (excluding "extra" dummies) */
1032:       if (dof[2] > 0 && exf < startexf + nexf && eyf < starteyf + neyf) {
1033:         DMStagStencil     rowc, colf;
1034:         PetscInt          ir, ic;
1035:         const PetscScalar cellScale = 0.25;

1037:         rowc.i   = exc;
1038:         rowc.j   = eyc;
1039:         rowc.c   = 0;
1040:         rowc.loc = DMSTAG_ELEMENT;
1041:         colf.i   = exf;
1042:         colf.j   = eyf;
1043:         colf.c   = 0;
1044:         colf.loc = DMSTAG_ELEMENT;
1045:         DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir);
1046:         DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic);
1047:         MatSetValuesLocal(A, 1, &ir, 1, &ic, &cellScale, INSERT_VALUES);
1048:       }
1049:     }
1050:   }
1051:   return 0;
1052: }

1054: PETSC_INTERN PetscErrorCode DMStagPopulateRestriction3d_0_0_a_b_Private(DM dmc, DM dmf, Mat A)
1055: {
1056:   PetscInt       exf, eyf, ezf, startexf, starteyf, startezf, nexf, neyf, nezf, nextraxf, nextrayf, nextrazf, startexc, starteyc, startezc, Nexf, Neyf, Nezf;
1057:   PetscInt       dof[4];
1058:   const PetscInt dim = 3;


1061:   DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], &dof[3]);

1065:   /* In 3D, each coarse point can receive from up to 12 fine points,
1066:      up to 4 of which may be off rank */
1067:   MatSeqAIJSetPreallocation(A, 12, NULL);
1068:   MatMPIAIJSetPreallocation(A, 12, NULL, 4, NULL);

1070:   DMStagGetCorners(dmf, &startexf, &starteyf, &startezf, &nexf, &neyf, &nezf, &nextraxf, &nextrayf, &nextrazf);
1071:   DMStagGetCorners(dmc, &startexc, &starteyc, &startezc, NULL, NULL, NULL, NULL, NULL, NULL);
1072:   DMStagGetGlobalSizes(dmf, &Nexf, &Neyf, &Nezf);

1074:   for (ezf = startezf; ezf < startezf + nezf + nextrazf; ++ezf) {
1075:     const PetscInt ezf_local = ezf - startezf;
1076:     const PetscInt ezc       = startezc + ezf_local / 2;

1078:     for (eyf = starteyf; eyf < starteyf + neyf + nextrayf; ++eyf) {
1079:       const PetscInt eyf_local = eyf - starteyf;
1080:       const PetscInt eyc       = starteyc + eyf_local / 2;

1082:       for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
1083:         const PetscInt exf_local = exf - startexf;
1084:         const PetscInt exc       = startexc + exf_local / 2;

1086:         /* Left faces (excluding top and front "extra" dummy layers) */
1087:         if (eyf < starteyf + neyf && ezf < startezf + nezf) {
1088:           DMStagStencil rowc[2], colf;
1089:           PetscInt      ir[2], ic, nweight;
1090:           PetscScalar   weight[2];

1092:           colf.i      = exf;
1093:           colf.j      = eyf;
1094:           colf.k      = ezf;
1095:           colf.c      = 0;
1096:           colf.loc    = DMSTAG_LEFT;
1097:           rowc[0].i   = exc;
1098:           rowc[0].j   = eyc;
1099:           rowc[0].k   = ezc;
1100:           rowc[0].c   = 0;
1101:           rowc[0].loc = DMSTAG_LEFT;
1102:           if (exf_local % 2 == 0) {
1103:             nweight = 1;
1104:             if (exf == Nexf || exf == 0) {
1105:               /* Note - this presumes something like a Neumann condition, assuming
1106:                  a ghost edge with the same value as the adjacent physical edge*/
1107:               weight[0] = 0.1875;
1108:             } else {
1109:               weight[0] = 0.125;
1110:             }
1111:           } else {
1112:             nweight     = 2;
1113:             rowc[1].i   = exc;
1114:             rowc[1].j   = eyc;
1115:             rowc[1].k   = ezc;
1116:             rowc[1].c   = 0;
1117:             rowc[1].loc = DMSTAG_RIGHT;
1118:             weight[0]   = 0.0625;
1119:             weight[1]   = 0.0625;
1120:           }
1121:           DMStagStencilToIndexLocal(dmc, dim, nweight, rowc, ir);
1122:           DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic);
1123:           MatSetValuesLocal(A, nweight, ir, 1, &ic, weight, INSERT_VALUES);
1124:         }

1126:         /* Down faces (excluding right and front "extra" dummy layers) */
1127:         if (exf < startexf + nexf && ezf < startezf + nezf) {
1128:           DMStagStencil rowc[2], colf;
1129:           PetscInt      ir[2], ic, nweight;
1130:           PetscScalar   weight[2];

1132:           colf.i      = exf;
1133:           colf.j      = eyf;
1134:           colf.k      = ezf;
1135:           colf.c      = 0;
1136:           colf.loc    = DMSTAG_DOWN;
1137:           rowc[0].i   = exc;
1138:           rowc[0].j   = eyc;
1139:           rowc[0].k   = ezc;
1140:           rowc[0].c   = 0;
1141:           rowc[0].loc = DMSTAG_DOWN;
1142:           if (eyf_local % 2 == 0) {
1143:             nweight = 1;
1144:             if (eyf == Neyf || eyf == 0) {
1145:               /* Note - this presumes something like a Neumann condition, assuming
1146:                  a ghost edge with the same value as the adjacent physical edge*/
1147:               weight[0] = 0.1875;
1148:             } else {
1149:               weight[0] = 0.125;
1150:             }
1151:           } else {
1152:             nweight     = 2;
1153:             rowc[1].i   = exc;
1154:             rowc[1].j   = eyc;
1155:             rowc[1].k   = ezc;
1156:             rowc[1].c   = 0;
1157:             rowc[1].loc = DMSTAG_UP;
1158:             weight[0]   = 0.0625;
1159:             weight[1]   = 0.0625;
1160:           }
1161:           DMStagStencilToIndexLocal(dmc, dim, nweight, rowc, ir);
1162:           DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic);
1163:           MatSetValuesLocal(A, nweight, ir, 1, &ic, weight, INSERT_VALUES);
1164:         }

1166:         /* Back faces (excluding left and top "extra" dummy laers) */
1167:         if (exf < startexf + nexf && eyf < starteyf + neyf) {
1168:           DMStagStencil rowc[2], colf;
1169:           PetscInt      ir[2], ic, nweight;
1170:           PetscScalar   weight[2];

1172:           colf.i      = exf;
1173:           colf.j      = eyf;
1174:           colf.k      = ezf;
1175:           colf.c      = 0;
1176:           colf.loc    = DMSTAG_BACK;
1177:           rowc[0].i   = exc;
1178:           rowc[0].j   = eyc;
1179:           rowc[0].k   = ezc;
1180:           rowc[0].c   = 0;
1181:           rowc[0].loc = DMSTAG_BACK;
1182:           if (ezf_local % 2 == 0) {
1183:             nweight = 1;
1184:             if (ezf == Nezf || ezf == 0) {
1185:               /* Note - this presumes something like a Neumann condition, assuming
1186:                  a ghost edge with the same value as the adjacent physical edge*/
1187:               weight[0] = 0.1875;
1188:             } else {
1189:               weight[0] = 0.125;
1190:             }
1191:           } else {
1192:             nweight     = 2;
1193:             rowc[1].i   = exc;
1194:             rowc[1].j   = eyc;
1195:             rowc[1].k   = ezc;
1196:             rowc[1].c   = 0;
1197:             rowc[1].loc = DMSTAG_FRONT;
1198:             weight[0]   = 0.0625;
1199:             weight[1]   = 0.0625;
1200:           }
1201:           DMStagStencilToIndexLocal(dmc, dim, nweight, rowc, ir);
1202:           DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic);
1203:           MatSetValuesLocal(A, nweight, ir, 1, &ic, weight, INSERT_VALUES);
1204:         }
1205:         /* Elements */
1206:         if (dof[3] == 1 && exf < startexf + nexf && eyf < starteyf + neyf && ezf < startezf + nezf) {
1207:           DMStagStencil     rowc, colf;
1208:           PetscInt          ir, ic;
1209:           const PetscScalar weight = 0.125;

1211:           colf.i   = exf;
1212:           colf.j   = eyf;
1213:           colf.k   = ezf;
1214:           colf.c   = 0;
1215:           colf.loc = DMSTAG_ELEMENT;
1216:           rowc.i   = exc;
1217:           rowc.j   = eyc;
1218:           rowc.k   = ezc;
1219:           rowc.c   = 0;
1220:           rowc.loc = DMSTAG_ELEMENT;
1221:           DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir);
1222:           DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic);
1223:           MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES);
1224:         }
1225:       }
1226:     }
1227:   }
1228:   return 0;
1229: }