Actual source code: stagmulti.c

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

  4: /*@
  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 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: [](ch_stag), `DMSTAG`, `DM`, `DMRestrict()`, `DMCoarsen()`, `DMCreateInjection()`
 23: @*/
 24: PetscErrorCode DMStagRestrictSimple(DM dmf, Vec xf, DM dmc, Vec xc)
 25: {
 26:   PetscInt dim;

 28:   PetscFunctionBegin;
 29:   PetscCall(DMGetDimension(dmf, &dim));
 30:   if (PetscDefined(USE_DEBUG)) {
 31:     PetscInt d, nf[DMSTAG_MAX_DIM], nc[DMSTAG_MAX_DIM], doff[DMSTAG_MAX_STRATA], dofc[DMSTAG_MAX_STRATA];

 33:     PetscCall(DMStagGetLocalSizes(dmf, &nf[0], &nf[1], &nf[2]));
 34:     PetscCall(DMStagGetLocalSizes(dmc, &nc[0], &nc[1], &nc[2]));
 35:     for (d = 0; d < dim; ++d) PetscCheck(nf[d] % nc[d] == 0, PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Not implemented for non-integer refinement factor");
 36:     PetscCall(DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]));
 37:     PetscCall(DMStagGetDOF(dmc, &dofc[0], &dofc[1], &dofc[2], &dofc[3]));
 38:     for (d = 0; d < dim + 1; ++d) PetscCheck(doff[d] == dofc[d], PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Cannot transfer between DMStag objects with different dof on each stratum");
 39:     {
 40:       PetscInt size_local, entries_local;

 42:       PetscCall(DMStagGetEntriesLocal(dmf, &entries_local));
 43:       PetscCall(VecGetLocalSize(xf, &size_local));
 44:       PetscCheck(entries_local == size_local, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Fine vector must be a local vector of size %" PetscInt_FMT ", but a vector of size %" PetscInt_FMT " was supplied", entries_local, size_local);
 45:     }
 46:     {
 47:       PetscInt size_local, entries_local;

 49:       PetscCall(DMStagGetEntriesLocal(dmc, &entries_local));
 50:       PetscCall(VecGetLocalSize(xc, &size_local));
 51:       PetscCheck(entries_local == size_local, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Coarse vector must be a local vector of size %" PetscInt_FMT ", but a vector of size %" PetscInt_FMT " was supplied", entries_local, size_local);
 52:     }
 53:   }
 54:   switch (dim) {
 55:   case 1:
 56:     PetscCall(DMStagRestrictSimple_1d(dmf, xf, dmc, xc));
 57:     break;
 58:   case 2:
 59:     PetscCall(DMStagRestrictSimple_2d(dmf, xf, dmc, xc));
 60:     break;
 61:   case 3:
 62:     PetscCall(DMStagRestrictSimple_3d(dmf, xf, dmc, xc));
 63:     break;
 64:   default:
 65:     SETERRQ(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
 66:     break;
 67:   }
 68:   PetscFunctionReturn(PETSC_SUCCESS);
 69: }

 71: static PetscErrorCode SetInterpolationCoefficientVertex_Private(PetscInt index, PetscInt factor, PetscScalar *a)
 72: {
 73:   PetscFunctionBegin;
 74:   *a = (index % factor) / (PetscScalar)factor;
 75:   PetscFunctionReturn(PETSC_SUCCESS);
 76: }

 78: static PetscErrorCode SetInterpolationCoefficientCenter_Private(PetscBool belowHalf, PetscInt index, PetscInt factor, PetscScalar *a)
 79: {
 80:   PetscFunctionBegin;
 81:   if (belowHalf) *a = 0.5 + ((index % factor) + 0.5) / (PetscScalar)factor;
 82:   else *a = ((index % factor) + 0.5) / (PetscScalar)factor - 0.5;
 83:   PetscFunctionReturn(PETSC_SUCCESS);
 84: }

 86: static PetscErrorCode SetInterpolationWeight1d_Private(PetscScalar ax, PetscScalar weight[])
 87: {
 88:   PetscFunctionBegin;
 89:   weight[0] = 1.0 - ax;
 90:   weight[1] = ax;
 91:   PetscFunctionReturn(PETSC_SUCCESS);
 92: }

 94: static PetscErrorCode SetInterpolationWeight2d_Private(PetscScalar ax, PetscScalar ay, PetscScalar weight[])
 95: {
 96:   PetscFunctionBegin;
 97:   weight[0] = (1.0 - ax) * (1.0 - ay);
 98:   weight[1] = ax * (1.0 - ay);
 99:   weight[2] = (1.0 - ax) * ay;
100:   weight[3] = ax * ay;
101:   PetscFunctionReturn(PETSC_SUCCESS);
102: }

104: static PetscErrorCode MergeInterpolationWeightToLeft2d_Private(PetscScalar weight[])
105: {
106:   PetscFunctionBegin;
107:   weight[0] += weight[1];
108:   weight[2] += weight[3];
109:   weight[1] = weight[3] = 0.0;
110:   PetscFunctionReturn(PETSC_SUCCESS);
111: }

113: static PetscErrorCode MergeInterpolationWeightToRight2d_Private(PetscScalar weight[])
114: {
115:   PetscFunctionBegin;
116:   weight[1] += weight[0];
117:   weight[3] += weight[2];
118:   weight[0] = weight[2] = 0.0;
119:   PetscFunctionReturn(PETSC_SUCCESS);
120: }

122: static PetscErrorCode MergeInterpolationWeightToBottom2d_Private(PetscScalar weight[])
123: {
124:   PetscFunctionBegin;
125:   weight[0] += weight[2];
126:   weight[1] += weight[3];
127:   weight[2] = weight[3] = 0.0;
128:   PetscFunctionReturn(PETSC_SUCCESS);
129: }

131: static PetscErrorCode MergeInterpolationWeightToTop2d_Private(PetscScalar weight[])
132: {
133:   PetscFunctionBegin;
134:   weight[2] += weight[0];
135:   weight[3] += weight[1];
136:   weight[0] = weight[1] = 0.0;
137:   PetscFunctionReturn(PETSC_SUCCESS);
138: }

140: static PetscErrorCode SetInterpolationWeight3d_Private(PetscScalar ax, PetscScalar ay, PetscScalar az, PetscScalar weight[])
141: {
142:   PetscFunctionBegin;
143:   weight[0] = (1.0 - ax) * (1.0 - ay) * (1.0 - az);
144:   weight[1] = ax * (1.0 - ay) * (1.0 - az);
145:   weight[2] = (1.0 - ax) * ay * (1.0 - az);
146:   weight[3] = ax * ay * (1.0 - az);
147:   weight[4] = (1.0 - ax) * (1.0 - ay) * az;
148:   weight[5] = ax * (1.0 - ay) * az;
149:   weight[6] = (1.0 - ax) * ay * az;
150:   weight[7] = ax * ay * az;
151:   PetscFunctionReturn(PETSC_SUCCESS);
152: }

154: static PetscErrorCode MergeInterpolationWeightToLeft3d_Private(PetscScalar weight[])
155: {
156:   PetscFunctionBegin;
157:   weight[0] += weight[1];
158:   weight[2] += weight[3];
159:   weight[4] += weight[5];
160:   weight[6] += weight[7];
161:   weight[1] = weight[3] = weight[5] = weight[7] = 0.0;
162:   PetscFunctionReturn(PETSC_SUCCESS);
163: }

165: static PetscErrorCode MergeInterpolationWeightToRight3d_Private(PetscScalar weight[])
166: {
167:   PetscFunctionBegin;
168:   weight[1] += weight[0];
169:   weight[3] += weight[2];
170:   weight[5] += weight[4];
171:   weight[7] += weight[6];
172:   weight[0] = weight[2] = weight[4] = weight[6] = 0.0;
173:   PetscFunctionReturn(PETSC_SUCCESS);
174: }

176: static PetscErrorCode MergeInterplationWeightToBottom3d_Private(PetscScalar weight[])
177: {
178:   PetscFunctionBegin;
179:   weight[0] += weight[2];
180:   weight[1] += weight[3];
181:   weight[4] += weight[6];
182:   weight[5] += weight[7];
183:   weight[2] = weight[3] = weight[6] = weight[7] = 0.0;
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

187: static PetscErrorCode MergeInterpolationWeightToTop3d_Private(PetscScalar weight[])
188: {
189:   PetscFunctionBegin;
190:   weight[2] += weight[0];
191:   weight[3] += weight[1];
192:   weight[6] += weight[4];
193:   weight[7] += weight[5];
194:   weight[0] = weight[1] = weight[4] = weight[5] = 0.0;
195:   PetscFunctionReturn(PETSC_SUCCESS);
196: }

198: static PetscErrorCode MergeInterpolationWeightToBack3d_Private(PetscScalar weight[])
199: {
200:   PetscFunctionBegin;
201:   weight[0] += weight[4];
202:   weight[1] += weight[5];
203:   weight[2] += weight[6];
204:   weight[3] += weight[7];
205:   weight[4] = weight[5] = weight[6] = weight[7] = 0.0;
206:   PetscFunctionReturn(PETSC_SUCCESS);
207: }

209: static PetscErrorCode MergeInterpolationWeightToFront3d_Private(PetscScalar weight[])
210: {
211:   PetscFunctionBegin;
212:   weight[4] += weight[0];
213:   weight[5] += weight[1];
214:   weight[6] += weight[2];
215:   weight[7] += weight[3];
216:   weight[0] = weight[1] = weight[2] = weight[3] = 0.0;
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: static PetscErrorCode SetRestrictionCoefficientVertex_Private(PetscInt index, PetscInt factor, PetscScalar *a)
221: {
222:   PetscFunctionBegin;
223:   *a = (factor - PetscAbsInt(index)) / (PetscScalar)(factor * factor);
224:   PetscFunctionReturn(PETSC_SUCCESS);
225: }

227: static PetscErrorCode SetRestrictionCoefficientCenter_Private(PetscInt index, PetscInt factor, PetscScalar *a)
228: {
229:   PetscFunctionBegin;
230:   if (2 * index + 1 < factor) *a = 0.5 + (index + 0.5) / (PetscScalar)factor;
231:   else *a = 1.5 - (index + 0.5) / (PetscScalar)factor;
232:   /* Normalization depends on whether the restriction factor is even or odd */
233:   if (factor % 2 == 0) *a /= 0.75 * factor;
234:   else *a /= (3 * factor * factor + 1) / (PetscScalar)(4 * factor);
235:   PetscFunctionReturn(PETSC_SUCCESS);
236: }

238: static PetscErrorCode MergeRestrictionWeightToLeft1d_Private(PetscScalar weight[], PetscInt m)
239: {
240:   PetscInt       i;
241:   const PetscInt dst = m / 2;

243:   PetscFunctionBegin;
244:   PetscCheck(m % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
245:   for (i = m / 2 + 1; i < m; ++i) {
246:     weight[dst] += weight[i];
247:     weight[i] = 0.0;
248:   }
249:   PetscFunctionReturn(PETSC_SUCCESS);
250: }

252: static PetscErrorCode MergeRestrictionWeightToRight1d_Private(PetscScalar weight[], PetscInt m)
253: {
254:   PetscInt       i;
255:   const PetscInt dst = m / 2;

257:   PetscFunctionBegin;
258:   PetscCheck(m % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
259:   for (i = 0; i < m / 2; ++i) {
260:     weight[dst] += weight[i];
261:     weight[i] = 0.0;
262:   }
263:   PetscFunctionReturn(PETSC_SUCCESS);
264: }

266: static PetscErrorCode MergeRestrictionWeightToLeft2d_Private(PetscScalar weight[], PetscInt m, PetscInt n)
267: {
268:   PetscInt i, j, src, dst;

270:   PetscFunctionBegin;
271:   PetscCheck(m % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
272:   for (j = 0; j < n; ++j)
273:     for (i = m / 2 + 1; i < m; ++i) {
274:       src = m * j + i;
275:       dst = m * j + m / 2;
276:       weight[dst] += weight[src];
277:       weight[src] = 0.0;
278:     }
279:   PetscFunctionReturn(PETSC_SUCCESS);
280: }

282: static PetscErrorCode MergeRestrictionWeightToRight2d_Private(PetscScalar weight[], PetscInt m, PetscInt n)
283: {
284:   PetscInt i, j, src, dst;

286:   PetscFunctionBegin;
287:   PetscCheck(m % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
288:   for (j = 0; j < n; ++j)
289:     for (i = 0; i < m / 2; ++i) {
290:       src = m * j + i;
291:       dst = m * j + m / 2;
292:       weight[dst] += weight[src];
293:       weight[src] = 0.0;
294:     }
295:   PetscFunctionReturn(PETSC_SUCCESS);
296: }

298: static PetscErrorCode MergeRestrictionWeightToBottom2d_Private(PetscScalar weight[], PetscInt m, PetscInt n)
299: {
300:   PetscInt i, j, src, dst;

302:   PetscFunctionBegin;
303:   PetscCheck(n % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
304:   for (j = n / 2 + 1; j < n; ++j)
305:     for (i = 0; i < m; ++i) {
306:       src = m * j + i;
307:       dst = m * (n / 2) + i;
308:       weight[dst] += weight[src];
309:       weight[src] = 0.0;
310:     }
311:   PetscFunctionReturn(PETSC_SUCCESS);
312: }

314: static PetscErrorCode MergeRestrictionWeightToTop2d_Private(PetscScalar weight[], PetscInt m, PetscInt n)
315: {
316:   PetscInt i, j, src, dst;

318:   PetscFunctionBegin;
319:   PetscCheck(n % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
320:   for (j = 0; j < n / 2; ++j)
321:     for (i = 0; i < m; ++i) {
322:       src = m * j + i;
323:       dst = m * (n / 2) + i;
324:       weight[dst] += weight[src];
325:       weight[src] = 0.0;
326:     }
327:   PetscFunctionReturn(PETSC_SUCCESS);
328: }

330: static PetscErrorCode MergeRestrictionWeightToLeft3d_Private(PetscScalar weight[], PetscInt m, PetscInt n, PetscInt p)
331: {
332:   PetscInt i, j, k, src, dst;

334:   PetscFunctionBegin;
335:   PetscCheck(m % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
336:   for (k = 0; k < p; ++k)
337:     for (j = 0; j < n; ++j)
338:       for (i = m / 2 + 1; i < m; ++i) {
339:         src = m * n * k + m * j + i;
340:         dst = m * n * k + m * j + m / 2;
341:         weight[dst] += weight[src];
342:         weight[src] = 0.0;
343:       }
344:   PetscFunctionReturn(PETSC_SUCCESS);
345: }

347: static PetscErrorCode MergeRestrictionWeightToRight3d_Private(PetscScalar weight[], PetscInt m, PetscInt n, PetscInt p)
348: {
349:   PetscInt i, j, k, src, dst;

351:   PetscFunctionBegin;
352:   PetscCheck(m % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
353:   for (k = 0; k < p; ++k)
354:     for (j = 0; j < n; ++j)
355:       for (i = 0; i < m / 2; ++i) {
356:         src = m * n * k + m * j + i;
357:         dst = m * n * k + m * j + m / 2;
358:         weight[dst] += weight[src];
359:         weight[src] = 0.0;
360:       }
361:   PetscFunctionReturn(PETSC_SUCCESS);
362: }

364: static PetscErrorCode MergeRestrictionWeightToBottom3d_Private(PetscScalar weight[], PetscInt m, PetscInt n, PetscInt p)
365: {
366:   PetscInt i, j, k, src, dst;

368:   PetscFunctionBegin;
369:   PetscCheck(n % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
370:   for (k = 0; k < p; ++k)
371:     for (j = n / 2 + 1; j < n; ++j)
372:       for (i = 0; i < m; ++i) {
373:         src = m * n * k + m * j + i;
374:         dst = m * n * k + m * (n / 2) + i;
375:         weight[dst] += weight[src];
376:         weight[src] = 0.0;
377:       }
378:   PetscFunctionReturn(PETSC_SUCCESS);
379: }

381: static PetscErrorCode MergeRestrictionWeightToTop3d_Private(PetscScalar weight[], PetscInt m, PetscInt n, PetscInt p)
382: {
383:   PetscInt i, j, k, src, dst;

385:   PetscFunctionBegin;
386:   PetscCheck(n % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
387:   for (k = 0; k < p; ++k)
388:     for (j = 0; j < n / 2; ++j)
389:       for (i = 0; i < m; ++i) {
390:         src = m * n * k + m * j + i;
391:         dst = m * n * k + m * (n / 2) + i;
392:         weight[dst] += weight[src];
393:         weight[src] = 0.0;
394:       }
395:   PetscFunctionReturn(PETSC_SUCCESS);
396: }

398: static PetscErrorCode MergeRestrictionWeightToBack3d_Private(PetscScalar weight[], PetscInt m, PetscInt n, PetscInt p)
399: {
400:   PetscInt i, j, k, src, dst;

402:   PetscFunctionBegin;
403:   PetscCheck(p % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
404:   for (k = p / 2 + 1; k < p; ++k)
405:     for (j = 0; j < n; ++j)
406:       for (i = 0; i < m; ++i) {
407:         src = m * n * k + m * j + i;
408:         dst = m * n * (p / 2) + m * j + i;
409:         weight[dst] += weight[src];
410:         weight[src] = 0.0;
411:       }
412:   PetscFunctionReturn(PETSC_SUCCESS);
413: }

415: static PetscErrorCode MergeRestrictionWeightToFront3d_Private(PetscScalar weight[], PetscInt m, PetscInt n, PetscInt p)
416: {
417:   PetscInt i, j, k, src, dst;

419:   PetscFunctionBegin;
420:   PetscCheck(p % 2 == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dimension must be odd");
421:   for (k = 0; k < p / 2; ++k)
422:     for (j = 0; j < n; ++j)
423:       for (i = 0; i < m; ++i) {
424:         src = m * n * k + m * j + i;
425:         dst = m * n * (p / 2) + m * j + i;
426:         weight[dst] += weight[src];
427:         weight[src] = 0.0;
428:       }
429:   PetscFunctionReturn(PETSC_SUCCESS);
430: }

432: static PetscErrorCode RemoveZeroWeights_Private(PetscInt n, DMStagStencil colc[], PetscScalar weight[], PetscInt *count)
433: {
434:   PetscInt i;

436:   PetscFunctionBegin;
437:   *count = 0;
438:   for (i = 0; i < n; ++i)
439:     if (weight[i] != 0.0) {
440:       colc[*count]   = colc[i];
441:       weight[*count] = weight[i];
442:       ++(*count);
443:     }
444:   PetscFunctionReturn(PETSC_SUCCESS);
445: }

447: /* Code duplication note: the next two functions are nearly identical, save the inclusion of the element terms */
448: PETSC_INTERN PetscErrorCode DMStagPopulateInterpolation1d_Internal(DM dmc, DM dmf, Mat A)
449: {
450:   PetscInt       Mc, Mf, factorx, dof[2];
451:   PetscInt       xf, mf, nExtraxf, i, d, count;
452:   DMStagStencil  rowf, colc[2];
453:   PetscScalar    ax, weight[2];
454:   PetscInt       ir, ic[2];
455:   const PetscInt dim = 1;

457:   PetscFunctionBegin;
458:   PetscCall(DMStagGetGlobalSizes(dmc, &Mc, NULL, NULL));
459:   PetscCall(DMStagGetGlobalSizes(dmf, &Mf, NULL, NULL));
460:   factorx = Mf / Mc;
461:   PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], NULL, NULL));

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

467:   PetscCall(DMStagGetCorners(dmf, &xf, NULL, NULL, &mf, NULL, NULL, &nExtraxf, NULL, NULL));

469:   /* Linear interpolation for vertices */
470:   for (d = 0; d < dof[0]; ++d)
471:     for (i = xf; i < xf + mf + nExtraxf; ++i) {
472:       rowf.i   = i;
473:       rowf.c   = d;
474:       rowf.loc = DMSTAG_LEFT;
475:       for (count = 0; count < 2; ++count) {
476:         colc[count].i = i / factorx;
477:         colc[count].c = d;
478:       }
479:       colc[0].loc = DMSTAG_LEFT;
480:       colc[1].loc = DMSTAG_RIGHT;
481:       PetscCall(SetInterpolationCoefficientVertex_Private(i, factorx, &ax));
482:       PetscCall(SetInterpolationWeight1d_Private(ax, weight));

484:       PetscCall(RemoveZeroWeights_Private(2, colc, weight, &count));
485:       PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
486:       PetscCall(DMStagStencilToIndexLocal(dmc, dim, count, colc, ic));
487:       PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
488:     }

490:   /* Nearest neighbor for elements */
491:   for (d = 0; d < dof[1]; ++d)
492:     for (i = xf; i < xf + mf; ++i) {
493:       rowf.i      = i;
494:       rowf.c      = d;
495:       rowf.loc    = DMSTAG_ELEMENT;
496:       colc[0].i   = i / factorx;
497:       colc[0].c   = d;
498:       colc[0].loc = DMSTAG_ELEMENT;
499:       weight[0]   = 1.0;
500:       PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
501:       PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, colc, ic));
502:       PetscCall(MatSetValuesLocal(A, 1, &ir, 1, ic, weight, INSERT_VALUES));
503:     }
504:   PetscFunctionReturn(PETSC_SUCCESS);
505: }

507: PETSC_INTERN PetscErrorCode DMStagPopulateInterpolation2d_Internal(DM dmc, DM dmf, Mat A)
508: {
509:   PetscInt       Mc, Nc, Mf, Nf, factorx, factory, dof[3];
510:   PetscInt       xf, yf, mf, nf, nExtraxf, nExtrayf, i, j, d, count;
511:   DMStagStencil  rowf, colc[4];
512:   PetscScalar    ax, ay, weight[4];
513:   PetscInt       ir, ic[4];
514:   const PetscInt dim = 2;

516:   PetscFunctionBegin;
517:   PetscCall(DMStagGetGlobalSizes(dmc, &Mc, &Nc, NULL));
518:   PetscCall(DMStagGetGlobalSizes(dmf, &Mf, &Nf, NULL));
519:   factorx = Mf / Mc;
520:   factory = Nf / Nc;
521:   PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], NULL));

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

527:   PetscCall(DMStagGetCorners(dmf, &xf, &yf, NULL, &mf, &nf, NULL, &nExtraxf, &nExtrayf, NULL));

529:   /* Linear interpolation for vertices */
530:   for (d = 0; d < dof[0]; ++d)
531:     for (j = yf; j < yf + nf + nExtrayf; ++j)
532:       for (i = xf; i < xf + mf + nExtraxf; ++i) {
533:         rowf.i   = i;
534:         rowf.j   = j;
535:         rowf.c   = d;
536:         rowf.loc = DMSTAG_DOWN_LEFT;
537:         for (count = 0; count < 4; ++count) {
538:           colc[count].i = i / factorx;
539:           colc[count].j = j / factory;
540:           colc[count].c = d;
541:         }
542:         colc[0].loc = DMSTAG_DOWN_LEFT;
543:         colc[1].loc = DMSTAG_DOWN_RIGHT;
544:         colc[2].loc = DMSTAG_UP_LEFT;
545:         colc[3].loc = DMSTAG_UP_RIGHT;
546:         PetscCall(SetInterpolationCoefficientVertex_Private(i, factorx, &ax));
547:         PetscCall(SetInterpolationCoefficientVertex_Private(j, factory, &ay));
548:         PetscCall(SetInterpolationWeight2d_Private(ax, ay, weight));

550:         PetscCall(RemoveZeroWeights_Private(4, colc, weight, &count));
551:         PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
552:         PetscCall(DMStagStencilToIndexLocal(dmc, dim, count, colc, ic));
553:         PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
554:       }

556:   /* Linear interpolation for left edges */
557:   for (d = 0; d < dof[1]; ++d)
558:     for (j = yf; j < yf + nf; ++j)
559:       for (i = xf; i < xf + mf + nExtraxf; ++i) {
560:         const PetscBool belowHalfy = (PetscBool)(2 * (j % factory) + 1 < factory);

562:         rowf.i   = i;
563:         rowf.j   = j;
564:         rowf.c   = d;
565:         rowf.loc = DMSTAG_LEFT;
566:         for (count = 0; count < 4; ++count) {
567:           colc[count].i = i / factorx;
568:           colc[count].j = j / factory;
569:           if (belowHalfy) colc[count].j -= 1;
570:           if (count / 2 == 1) colc[count].j += 1;
571:           colc[count].c = d;
572:         }
573:         colc[0].loc = DMSTAG_LEFT;
574:         colc[1].loc = DMSTAG_RIGHT;
575:         colc[2].loc = DMSTAG_LEFT;
576:         colc[3].loc = DMSTAG_RIGHT;
577:         PetscCall(SetInterpolationCoefficientVertex_Private(i, factorx, &ax));
578:         PetscCall(SetInterpolationCoefficientCenter_Private(belowHalfy, j, factory, &ay));
579:         PetscCall(SetInterpolationWeight2d_Private(ax, ay, weight));
580:         /* Assume Neumann boundary condition */
581:         if (j / factory == 0 && belowHalfy) PetscCall(MergeInterpolationWeightToTop2d_Private(weight));
582:         else if (j / factory == Nc - 1 && !belowHalfy) PetscCall(MergeInterpolationWeightToBottom2d_Private(weight));

584:         PetscCall(RemoveZeroWeights_Private(4, colc, weight, &count));
585:         PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
586:         PetscCall(DMStagStencilToIndexLocal(dmc, dim, count, colc, ic));
587:         PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
588:       }

590:   /* Linear interpolation for down edges */
591:   for (d = 0; d < dof[1]; ++d)
592:     for (j = yf; j < yf + nf + nExtrayf; ++j)
593:       for (i = xf; i < xf + mf; ++i) {
594:         const PetscBool belowHalfx = (PetscBool)(2 * (i % factorx) + 1 < factorx);

596:         rowf.i   = i;
597:         rowf.j   = j;
598:         rowf.c   = d;
599:         rowf.loc = DMSTAG_DOWN;
600:         for (count = 0; count < 4; ++count) {
601:           colc[count].i = i / factorx;
602:           colc[count].j = j / factory;
603:           if (belowHalfx) colc[count].i -= 1;
604:           if (count % 2 == 1) colc[count].i += 1;
605:           colc[count].c = d;
606:         }
607:         colc[0].loc = DMSTAG_DOWN;
608:         colc[1].loc = DMSTAG_DOWN;
609:         colc[2].loc = DMSTAG_UP;
610:         colc[3].loc = DMSTAG_UP;
611:         PetscCall(SetInterpolationCoefficientCenter_Private(belowHalfx, i, factorx, &ax));
612:         PetscCall(SetInterpolationCoefficientVertex_Private(j, factory, &ay));
613:         PetscCall(SetInterpolationWeight2d_Private(ax, ay, weight));
614:         /* Assume Neumann boundary condition */
615:         if (i / factorx == 0 && belowHalfx) PetscCall(MergeInterpolationWeightToRight2d_Private(weight));
616:         else if (i / factorx == Mc - 1 && !belowHalfx) PetscCall(MergeInterpolationWeightToLeft2d_Private(weight));

618:         PetscCall(RemoveZeroWeights_Private(4, colc, weight, &count));
619:         PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
620:         PetscCall(DMStagStencilToIndexLocal(dmc, dim, count, colc, ic));
621:         PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
622:       }

624:   /* Nearest neighbor for elements */
625:   for (d = 0; d < dof[2]; ++d)
626:     for (j = yf; j < yf + nf; ++j)
627:       for (i = xf; i < xf + mf; ++i) {
628:         rowf.i      = i;
629:         rowf.j      = j;
630:         rowf.c      = d;
631:         rowf.loc    = DMSTAG_ELEMENT;
632:         colc[0].i   = i / factorx;
633:         colc[0].j   = j / factory;
634:         colc[0].c   = d;
635:         colc[0].loc = DMSTAG_ELEMENT;
636:         weight[0]   = 1.0;
637:         PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
638:         PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, colc, ic));
639:         PetscCall(MatSetValuesLocal(A, 1, &ir, 1, ic, weight, INSERT_VALUES));
640:       }
641:   PetscFunctionReturn(PETSC_SUCCESS);
642: }

644: PETSC_INTERN PetscErrorCode DMStagPopulateInterpolation3d_Internal(DM dmc, DM dmf, Mat A)
645: {
646:   PetscInt       Mc, Nc, Pc, Mf, Nf, Pf, factorx, factory, factorz, dof[4];
647:   PetscInt       xf, yf, zf, mf, nf, pf, nExtraxf, nExtrayf, nExtrazf, i, j, k, d, count;
648:   DMStagStencil  rowf, colc[8];
649:   PetscScalar    ax, ay, az, weight[8];
650:   PetscInt       ir, ic[8];
651:   const PetscInt dim = 3;

653:   PetscFunctionBegin;
654:   PetscCall(DMStagGetGlobalSizes(dmc, &Mc, &Nc, &Pc));
655:   PetscCall(DMStagGetGlobalSizes(dmf, &Mf, &Nf, &Pf));
656:   factorx = Mf / Mc;
657:   factory = Nf / Nc;
658:   factorz = Pf / Pc;
659:   PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], &dof[3]));

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

665:   PetscCall(DMStagGetCorners(dmf, &xf, &yf, &zf, &mf, &nf, &pf, &nExtraxf, &nExtrayf, &nExtrazf));

667:   /* Linear interpolation for vertices */
668:   for (d = 0; d < dof[0]; ++d)
669:     for (k = zf; k < zf + pf + nExtrazf; ++k)
670:       for (j = yf; j < yf + nf + nExtrayf; ++j)
671:         for (i = xf; i < xf + mf + nExtraxf; ++i) {
672:           rowf.i   = i;
673:           rowf.j   = j;
674:           rowf.k   = k;
675:           rowf.c   = d;
676:           rowf.loc = DMSTAG_BACK_DOWN_LEFT;
677:           for (count = 0; count < 8; ++count) {
678:             colc[count].i = i / factorx;
679:             colc[count].j = j / factory;
680:             colc[count].k = k / factorz;
681:             colc[count].c = d;
682:           }
683:           colc[0].loc = DMSTAG_BACK_DOWN_LEFT;
684:           colc[1].loc = DMSTAG_BACK_DOWN_RIGHT;
685:           colc[2].loc = DMSTAG_BACK_UP_LEFT;
686:           colc[3].loc = DMSTAG_BACK_UP_RIGHT;
687:           colc[4].loc = DMSTAG_FRONT_DOWN_LEFT;
688:           colc[5].loc = DMSTAG_FRONT_DOWN_RIGHT;
689:           colc[6].loc = DMSTAG_FRONT_UP_LEFT;
690:           colc[7].loc = DMSTAG_FRONT_UP_RIGHT;
691:           PetscCall(SetInterpolationCoefficientVertex_Private(i, factorx, &ax));
692:           PetscCall(SetInterpolationCoefficientVertex_Private(j, factory, &ay));
693:           PetscCall(SetInterpolationCoefficientVertex_Private(k, factorz, &az));
694:           PetscCall(SetInterpolationWeight3d_Private(ax, ay, az, weight));

696:           PetscCall(RemoveZeroWeights_Private(8, colc, weight, &count));
697:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
698:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, count, colc, ic));
699:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
700:         }

702:   /* Linear interpolation for down left edges */
703:   for (d = 0; d < dof[1]; ++d)
704:     for (k = zf; k < zf + pf; ++k)
705:       for (j = yf; j < yf + nf + nExtrayf; ++j)
706:         for (i = xf; i < xf + mf + nExtraxf; ++i) {
707:           const PetscBool belowHalfz = (PetscBool)(2 * (k % factorz) + 1 < factorz);

709:           rowf.i   = i;
710:           rowf.j   = j;
711:           rowf.k   = k;
712:           rowf.c   = d;
713:           rowf.loc = DMSTAG_DOWN_LEFT;
714:           for (count = 0; count < 8; ++count) {
715:             colc[count].i = i / factorx;
716:             colc[count].j = j / factory;
717:             colc[count].k = k / factorz;
718:             if (belowHalfz) colc[count].k -= 1;
719:             if (count / 4 == 1) colc[count].k += 1;
720:             colc[count].c = d;
721:           }
722:           colc[0].loc = DMSTAG_DOWN_LEFT;
723:           colc[1].loc = DMSTAG_DOWN_RIGHT;
724:           colc[2].loc = DMSTAG_UP_LEFT;
725:           colc[3].loc = DMSTAG_UP_RIGHT;
726:           colc[4].loc = DMSTAG_DOWN_LEFT;
727:           colc[5].loc = DMSTAG_DOWN_RIGHT;
728:           colc[6].loc = DMSTAG_UP_LEFT;
729:           colc[7].loc = DMSTAG_UP_RIGHT;
730:           PetscCall(SetInterpolationCoefficientVertex_Private(i, factorx, &ax));
731:           PetscCall(SetInterpolationCoefficientVertex_Private(j, factory, &ay));
732:           PetscCall(SetInterpolationCoefficientCenter_Private(belowHalfz, k, factorz, &az));
733:           PetscCall(SetInterpolationWeight3d_Private(ax, ay, az, weight));
734:           /* Assume Neumann boundary condition */
735:           if (k / factorz == 0 && belowHalfz) PetscCall(MergeInterpolationWeightToFront3d_Private(weight));
736:           else if (k / factorz == Pc - 1 && !belowHalfz) PetscCall(MergeInterpolationWeightToBack3d_Private(weight));

738:           PetscCall(RemoveZeroWeights_Private(8, colc, weight, &count));
739:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
740:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, count, colc, ic));
741:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
742:         }

744:   /* Linear interpolation for back left edges */
745:   for (d = 0; d < dof[1]; ++d)
746:     for (k = zf; k < zf + pf + nExtrazf; ++k)
747:       for (j = yf; j < yf + nf; ++j)
748:         for (i = xf; i < xf + mf + nExtraxf; ++i) {
749:           const PetscBool belowHalfy = (PetscBool)(2 * (j % factory) + 1 < factory);

751:           rowf.i   = i;
752:           rowf.j   = j;
753:           rowf.k   = k;
754:           rowf.c   = d;
755:           rowf.loc = DMSTAG_BACK_LEFT;
756:           for (count = 0; count < 8; ++count) {
757:             colc[count].i = i / factorx;
758:             colc[count].j = j / factory;
759:             colc[count].k = k / factorz;
760:             if (belowHalfy) colc[count].j -= 1;
761:             if ((count % 4) / 2 == 1) colc[count].j += 1;
762:             colc[count].c = d;
763:           }
764:           colc[0].loc = DMSTAG_BACK_LEFT;
765:           colc[1].loc = DMSTAG_BACK_RIGHT;
766:           colc[2].loc = DMSTAG_BACK_LEFT;
767:           colc[3].loc = DMSTAG_BACK_RIGHT;
768:           colc[4].loc = DMSTAG_FRONT_LEFT;
769:           colc[5].loc = DMSTAG_FRONT_RIGHT;
770:           colc[6].loc = DMSTAG_FRONT_LEFT;
771:           colc[7].loc = DMSTAG_FRONT_RIGHT;
772:           PetscCall(SetInterpolationCoefficientVertex_Private(i, factorx, &ax));
773:           PetscCall(SetInterpolationCoefficientCenter_Private(belowHalfy, j, factory, &ay));
774:           PetscCall(SetInterpolationCoefficientVertex_Private(k, factorz, &az));
775:           PetscCall(SetInterpolationWeight3d_Private(ax, ay, az, weight));
776:           /* Assume Neumann boundary condition */
777:           if (j / factory == 0 && belowHalfy) PetscCall(MergeInterpolationWeightToTop3d_Private(weight));
778:           else if (j / factory == Nc - 1 && !belowHalfy) PetscCall(MergeInterplationWeightToBottom3d_Private(weight));

780:           PetscCall(RemoveZeroWeights_Private(8, colc, weight, &count));
781:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
782:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, count, colc, ic));
783:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
784:         }

786:   /* Linear interpolation for down back edges */
787:   for (d = 0; d < dof[1]; ++d)
788:     for (k = zf; k < zf + pf + nExtrazf; ++k)
789:       for (j = yf; j < yf + nf + nExtrayf; ++j)
790:         for (i = xf; i < xf + mf; ++i) {
791:           const PetscBool belowHalfx = (PetscBool)(2 * (i % factorx) + 1 < factorx);

793:           rowf.i   = i;
794:           rowf.j   = j;
795:           rowf.k   = k;
796:           rowf.c   = d;
797:           rowf.loc = DMSTAG_BACK_DOWN;
798:           for (count = 0; count < 8; ++count) {
799:             colc[count].i = i / factorx;
800:             colc[count].j = j / factory;
801:             colc[count].k = k / factorz;
802:             if (belowHalfx) colc[count].i -= 1;
803:             if (count % 2 == 1) colc[count].i += 1;
804:             colc[count].c = d;
805:           }
806:           colc[0].loc = DMSTAG_BACK_DOWN;
807:           colc[1].loc = DMSTAG_BACK_DOWN;
808:           colc[2].loc = DMSTAG_BACK_UP;
809:           colc[3].loc = DMSTAG_BACK_UP;
810:           colc[4].loc = DMSTAG_FRONT_DOWN;
811:           colc[5].loc = DMSTAG_FRONT_DOWN;
812:           colc[6].loc = DMSTAG_FRONT_UP;
813:           colc[7].loc = DMSTAG_FRONT_UP;
814:           PetscCall(SetInterpolationCoefficientCenter_Private(belowHalfx, i, factorx, &ax));
815:           PetscCall(SetInterpolationCoefficientVertex_Private(j, factory, &ay));
816:           PetscCall(SetInterpolationCoefficientVertex_Private(k, factorz, &az));
817:           PetscCall(SetInterpolationWeight3d_Private(ax, ay, az, weight));
818:           /* Assume Neumann boundary condition */
819:           if (i / factorx == 0 && belowHalfx) PetscCall(MergeInterpolationWeightToRight3d_Private(weight));
820:           else if (i / factorx == Mc - 1 && !belowHalfx) PetscCall(MergeInterpolationWeightToLeft3d_Private(weight));

822:           PetscCall(RemoveZeroWeights_Private(8, colc, weight, &count));
823:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
824:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, count, colc, ic));
825:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
826:         }

828:   /* Linear interpolation for left faces */
829:   for (d = 0; d < dof[2]; ++d)
830:     for (k = zf; k < zf + pf; ++k)
831:       for (j = yf; j < yf + nf; ++j)
832:         for (i = xf; i < xf + mf + nExtraxf; ++i) {
833:           const PetscBool belowHalfy = (PetscBool)(2 * (j % factory) + 1 < factory);
834:           const PetscBool belowHalfz = (PetscBool)(2 * (k % factorz) + 1 < factorz);

836:           rowf.i   = i;
837:           rowf.j   = j;
838:           rowf.k   = k;
839:           rowf.c   = d;
840:           rowf.loc = DMSTAG_LEFT;
841:           for (count = 0; count < 8; ++count) {
842:             colc[count].i = i / factorx;
843:             colc[count].j = j / factory;
844:             colc[count].k = k / factorz;
845:             if (belowHalfy) colc[count].j -= 1;
846:             if ((count % 4) / 2 == 1) colc[count].j += 1;
847:             if (belowHalfz) colc[count].k -= 1;
848:             if (count / 4 == 1) colc[count].k += 1;
849:             colc[count].c = d;
850:           }
851:           colc[0].loc = DMSTAG_LEFT;
852:           colc[1].loc = DMSTAG_RIGHT;
853:           colc[2].loc = DMSTAG_LEFT;
854:           colc[3].loc = DMSTAG_RIGHT;
855:           colc[4].loc = DMSTAG_LEFT;
856:           colc[5].loc = DMSTAG_RIGHT;
857:           colc[6].loc = DMSTAG_LEFT;
858:           colc[7].loc = DMSTAG_RIGHT;
859:           PetscCall(SetInterpolationCoefficientVertex_Private(i, factorx, &ax));
860:           PetscCall(SetInterpolationCoefficientCenter_Private(belowHalfy, j, factory, &ay));
861:           PetscCall(SetInterpolationCoefficientCenter_Private(belowHalfz, k, factorz, &az));
862:           PetscCall(SetInterpolationWeight3d_Private(ax, ay, az, weight));
863:           /* Assume Neumann boundary condition */
864:           if (j / factory == 0 && belowHalfy) PetscCall(MergeInterpolationWeightToTop3d_Private(weight));
865:           else if (j / factory == Nc - 1 && !belowHalfy) PetscCall(MergeInterplationWeightToBottom3d_Private(weight));
866:           if (k / factorz == 0 && belowHalfz) PetscCall(MergeInterpolationWeightToFront3d_Private(weight));
867:           else if (k / factorz == Pc - 1 && !belowHalfz) PetscCall(MergeInterpolationWeightToBack3d_Private(weight));

869:           PetscCall(RemoveZeroWeights_Private(8, colc, weight, &count));
870:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
871:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, count, colc, ic));
872:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
873:         }

875:   /* Linear interpolation for down faces */
876:   for (d = 0; d < dof[2]; ++d)
877:     for (k = zf; k < zf + pf; ++k)
878:       for (j = yf; j < yf + nf + nExtrayf; ++j)
879:         for (i = xf; i < xf + mf; ++i) {
880:           const PetscBool belowHalfx = (PetscBool)(2 * (i % factorx) + 1 < factorx);
881:           const PetscBool belowHalfz = (PetscBool)(2 * (k % factorz) + 1 < factorz);

883:           rowf.i   = i;
884:           rowf.j   = j;
885:           rowf.k   = k;
886:           rowf.c   = d;
887:           rowf.loc = DMSTAG_DOWN;
888:           for (count = 0; count < 8; ++count) {
889:             colc[count].i = i / factorx;
890:             colc[count].j = j / factory;
891:             colc[count].k = k / factorz;
892:             if (belowHalfx) colc[count].i -= 1;
893:             if (count % 2 == 1) colc[count].i += 1;
894:             if (belowHalfz) colc[count].k -= 1;
895:             if (count / 4 == 1) colc[count].k += 1;
896:             colc[count].c = d;
897:           }
898:           colc[0].loc = DMSTAG_DOWN;
899:           colc[1].loc = DMSTAG_DOWN;
900:           colc[2].loc = DMSTAG_UP;
901:           colc[3].loc = DMSTAG_UP;
902:           colc[4].loc = DMSTAG_DOWN;
903:           colc[5].loc = DMSTAG_DOWN;
904:           colc[6].loc = DMSTAG_UP;
905:           colc[7].loc = DMSTAG_UP;
906:           PetscCall(SetInterpolationCoefficientCenter_Private(belowHalfx, i, factorx, &ax));
907:           PetscCall(SetInterpolationCoefficientVertex_Private(j, factory, &ay));
908:           PetscCall(SetInterpolationCoefficientCenter_Private(belowHalfz, k, factorz, &az));
909:           PetscCall(SetInterpolationWeight3d_Private(ax, ay, az, weight));
910:           /* Assume Neumann boundary condition */
911:           if (i / factorx == 0 && belowHalfx) PetscCall(MergeInterpolationWeightToRight3d_Private(weight));
912:           else if (i / factorx == Mc - 1 && !belowHalfx) PetscCall(MergeInterpolationWeightToLeft3d_Private(weight));
913:           if (k / factorz == 0 && belowHalfz) PetscCall(MergeInterpolationWeightToFront3d_Private(weight));
914:           else if (k / factorz == Pc - 1 && !belowHalfz) PetscCall(MergeInterpolationWeightToBack3d_Private(weight));

916:           PetscCall(RemoveZeroWeights_Private(8, colc, weight, &count));
917:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
918:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, count, colc, ic));
919:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
920:         }

922:   /* Linear interpolation for back faces */
923:   for (d = 0; d < dof[2]; ++d)
924:     for (k = zf; k < zf + pf + nExtrazf; ++k)
925:       for (j = yf; j < yf + nf; ++j)
926:         for (i = xf; i < xf + mf; ++i) {
927:           const PetscBool belowHalfx = (PetscBool)(2 * (i % factorx) + 1 < factorx);
928:           const PetscBool belowHalfy = (PetscBool)(2 * (j % factory) + 1 < factory);

930:           rowf.i   = i;
931:           rowf.j   = j;
932:           rowf.k   = k;
933:           rowf.c   = d;
934:           rowf.loc = DMSTAG_BACK;
935:           for (count = 0; count < 8; ++count) {
936:             colc[count].i = i / factorx;
937:             colc[count].j = j / factory;
938:             colc[count].k = k / factorz;
939:             if (belowHalfx) colc[count].i -= 1;
940:             if (count % 2 == 1) colc[count].i += 1;
941:             if (belowHalfy) colc[count].j -= 1;
942:             if ((count % 4) / 2 == 1) colc[count].j += 1;
943:             colc[count].c = d;
944:           }
945:           colc[0].loc = DMSTAG_BACK;
946:           colc[1].loc = DMSTAG_BACK;
947:           colc[2].loc = DMSTAG_BACK;
948:           colc[3].loc = DMSTAG_BACK;
949:           colc[4].loc = DMSTAG_FRONT;
950:           colc[5].loc = DMSTAG_FRONT;
951:           colc[6].loc = DMSTAG_FRONT;
952:           colc[7].loc = DMSTAG_FRONT;
953:           PetscCall(SetInterpolationCoefficientCenter_Private(belowHalfx, i, factorx, &ax));
954:           PetscCall(SetInterpolationCoefficientCenter_Private(belowHalfy, j, factory, &ay));
955:           PetscCall(SetInterpolationCoefficientVertex_Private(k, factorz, &az));
956:           PetscCall(SetInterpolationWeight3d_Private(ax, ay, az, weight));
957:           /* Assume Neumann boundary condition */
958:           if (i / factorx == 0 && belowHalfx) PetscCall(MergeInterpolationWeightToRight3d_Private(weight));
959:           else if (i / factorx == Mc - 1 && !belowHalfx) PetscCall(MergeInterpolationWeightToLeft3d_Private(weight));
960:           if (j / factory == 0 && belowHalfy) PetscCall(MergeInterpolationWeightToTop3d_Private(weight));
961:           else if (j / factory == Nc - 1 && !belowHalfy) PetscCall(MergeInterplationWeightToBottom3d_Private(weight));

963:           PetscCall(RemoveZeroWeights_Private(8, colc, weight, &count));
964:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
965:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, count, colc, ic));
966:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
967:         }

969:   /* Nearest neighbor for elements */
970:   for (d = 0; d < dof[3]; ++d)
971:     for (k = zf; k < zf + pf; ++k)
972:       for (j = yf; j < yf + nf; ++j)
973:         for (i = xf; i < xf + mf; ++i) {
974:           rowf.i      = i;
975:           rowf.j      = j;
976:           rowf.k      = k;
977:           rowf.c      = d;
978:           rowf.loc    = DMSTAG_ELEMENT;
979:           colc[0].i   = i / factorx;
980:           colc[0].j   = j / factory;
981:           colc[0].k   = k / factorz;
982:           colc[0].c   = d;
983:           colc[0].loc = DMSTAG_ELEMENT;
984:           weight[0]   = 1.0;
985:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
986:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, colc, ic));
987:           PetscCall(MatSetValuesLocal(A, 1, &ir, 1, ic, weight, INSERT_VALUES));
988:         }
989:   PetscFunctionReturn(PETSC_SUCCESS);
990: }

992: PETSC_INTERN PetscErrorCode DMStagPopulateRestriction1d_Internal(DM dmc, DM dmf, Mat A)
993: {
994:   PetscInt       Mc, Mf, factorx, dof[2];
995:   PetscInt       xc, mc, nExtraxc, i, d, ii, count;
996:   PetscInt       maxFinePoints, maxOffRankFinePoints;
997:   DMStagStencil  rowc;
998:   DMStagStencil *colf;
999:   PetscScalar   *weight;
1000:   PetscInt       ir;
1001:   PetscInt      *ic;
1002:   const PetscInt dim = 1;

1004:   PetscFunctionBegin;
1005:   PetscCall(DMStagGetGlobalSizes(dmc, &Mc, NULL, NULL));
1006:   PetscCall(DMStagGetGlobalSizes(dmf, &Mf, NULL, NULL));
1007:   factorx = Mf / Mc;
1008:   PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], NULL, NULL));

1010:   /* In 1D, each coarse point can receive from up to (2 * factorx - 1) fine points, (factorx - 1) of which may be off-rank */
1011:   maxFinePoints        = 2 * factorx - 1;
1012:   maxOffRankFinePoints = maxFinePoints - factorx;
1013:   PetscCall(MatSeqAIJSetPreallocation(A, maxFinePoints, NULL));
1014:   PetscCall(MatMPIAIJSetPreallocation(A, maxFinePoints, NULL, maxOffRankFinePoints, NULL));
1015:   PetscCall(PetscMalloc3(maxFinePoints, &colf, maxFinePoints, &weight, maxFinePoints, &ic));

1017:   PetscCall(DMStagGetCorners(dmc, &xc, NULL, NULL, &mc, NULL, NULL, &nExtraxc, NULL, NULL));

1019:   for (d = 0; d < dof[0]; ++d)
1020:     for (i = xc; i < xc + mc + nExtraxc; ++i) {
1021:       rowc.i   = i;
1022:       rowc.c   = d;
1023:       rowc.loc = DMSTAG_LEFT;
1024:       count    = 0;
1025:       for (ii = -(factorx - 1); ii <= factorx - 1; ++ii) {
1026:         colf[count].i   = i * factorx + ii;
1027:         colf[count].c   = d;
1028:         colf[count].loc = DMSTAG_LEFT;
1029:         PetscCall(SetRestrictionCoefficientVertex_Private(ii, factorx, &weight[count]));
1030:         ++count;
1031:       }
1032:       if (i == 0) PetscCall(MergeRestrictionWeightToRight1d_Private(weight, 2 * factorx - 1));
1033:       else if (i == Mc) PetscCall(MergeRestrictionWeightToLeft1d_Private(weight, 2 * factorx - 1));

1035:       PetscCall(RemoveZeroWeights_Private(2 * factorx - 1, colf, weight, &count));
1036:       PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1037:       PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1038:       PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1039:     }

1041:   for (d = 0; d < dof[1]; ++d)
1042:     for (i = xc; i < xc + mc; ++i) {
1043:       rowc.i   = i;
1044:       rowc.c   = d;
1045:       rowc.loc = DMSTAG_ELEMENT;
1046:       count    = 0;
1047:       for (ii = 0; ii < factorx; ++ii) {
1048:         colf[count].i   = i * factorx + ii;
1049:         colf[count].c   = d;
1050:         colf[count].loc = DMSTAG_ELEMENT;
1051:         weight[count]   = 1 / (PetscScalar)factorx;
1052:         ++count;
1053:       }

1055:       PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1056:       PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1057:       PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1058:     }

1060:   PetscCall(PetscFree3(colf, weight, ic));
1061:   PetscFunctionReturn(PETSC_SUCCESS);
1062: }

1064: PETSC_INTERN PetscErrorCode DMStagPopulateRestriction2d_Internal(DM dmc, DM dmf, Mat A)
1065: {
1066:   PetscInt       Mc, Nc, Mf, Nf, factorx, factory, dof[3];
1067:   PetscInt       xc, yc, mc, nc, nExtraxc, nExtrayc, i, j, d, ii, jj, count;
1068:   PetscInt       maxFinePoints, maxOffRankFinePoints;
1069:   DMStagStencil  rowc;
1070:   DMStagStencil *colf;
1071:   PetscScalar    ax, ay;
1072:   PetscScalar   *weight;
1073:   PetscInt       ir;
1074:   PetscInt      *ic;
1075:   const PetscInt dim = 2;

1077:   PetscFunctionBegin;
1078:   PetscCall(DMStagGetGlobalSizes(dmc, &Mc, &Nc, NULL));
1079:   PetscCall(DMStagGetGlobalSizes(dmf, &Mf, &Nf, NULL));
1080:   factorx = Mf / Mc;
1081:   factory = Nf / Nc;
1082:   PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], NULL));

1084:   /* In 2D, each coarse point can receive from up to ((2 * factorx - 1) * (2 * factory - 1)) fine points,
1085:      up to ((2 * factorx - 1) * (2 * factory - 1) - factorx * factory) of which may be off rank */
1086:   maxFinePoints        = (2 * factorx - 1) * (2 * factory - 1);
1087:   maxOffRankFinePoints = maxFinePoints - factorx * factory;
1088:   PetscCall(MatSeqAIJSetPreallocation(A, maxFinePoints, NULL));
1089:   PetscCall(MatMPIAIJSetPreallocation(A, maxFinePoints, NULL, maxOffRankFinePoints, NULL));
1090:   PetscCall(PetscMalloc3(maxFinePoints, &colf, maxFinePoints, &weight, maxFinePoints, &ic));

1092:   PetscCall(DMStagGetCorners(dmc, &xc, &yc, NULL, &mc, &nc, NULL, &nExtraxc, &nExtrayc, NULL));

1094:   for (d = 0; d < dof[0]; ++d)
1095:     for (j = yc; j < yc + nc + nExtrayc; ++j)
1096:       for (i = xc; i < xc + mc + nExtraxc; ++i) {
1097:         rowc.i   = i;
1098:         rowc.j   = j;
1099:         rowc.c   = d;
1100:         rowc.loc = DMSTAG_DOWN_LEFT;
1101:         count    = 0;
1102:         for (jj = -(factory - 1); jj <= factory - 1; ++jj)
1103:           for (ii = -(factorx - 1); ii <= factorx - 1; ++ii) {
1104:             colf[count].i   = i * factorx + ii;
1105:             colf[count].j   = j * factory + jj;
1106:             colf[count].c   = d;
1107:             colf[count].loc = DMSTAG_DOWN_LEFT;
1108:             PetscCall(SetRestrictionCoefficientVertex_Private(ii, factorx, &ax));
1109:             PetscCall(SetRestrictionCoefficientVertex_Private(jj, factory, &ay));
1110:             weight[count] = ax * ay;
1111:             ++count;
1112:           }
1113:         if (i == 0) PetscCall(MergeRestrictionWeightToRight2d_Private(weight, 2 * factorx - 1, 2 * factory - 1));
1114:         else if (i == Mc) PetscCall(MergeRestrictionWeightToLeft2d_Private(weight, 2 * factorx - 1, 2 * factory - 1));
1115:         if (j == 0) PetscCall(MergeRestrictionWeightToTop2d_Private(weight, 2 * factorx - 1, 2 * factory - 1));
1116:         else if (j == Nc) PetscCall(MergeRestrictionWeightToBottom2d_Private(weight, 2 * factorx - 1, 2 * factory - 1));

1118:         PetscCall(RemoveZeroWeights_Private((2 * factorx - 1) * (2 * factory - 1), colf, weight, &count));
1119:         PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1120:         PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1121:         PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1122:       }

1124:   for (d = 0; d < dof[1]; ++d)
1125:     for (j = yc; j < yc + nc; ++j)
1126:       for (i = xc; i < xc + mc + nExtraxc; ++i) {
1127:         rowc.i   = i;
1128:         rowc.j   = j;
1129:         rowc.c   = d;
1130:         rowc.loc = DMSTAG_LEFT;
1131:         count    = 0;
1132:         for (jj = 0; jj < factory; ++jj)
1133:           for (ii = -(factorx - 1); ii <= factorx - 1; ++ii) {
1134:             colf[count].i   = i * factorx + ii;
1135:             colf[count].j   = j * factory + jj;
1136:             colf[count].c   = d;
1137:             colf[count].loc = DMSTAG_LEFT;
1138:             PetscCall(SetRestrictionCoefficientVertex_Private(ii, factorx, &ax));
1139:             PetscCall(SetRestrictionCoefficientCenter_Private(jj, factory, &ay));
1140:             weight[count] = ax * ay;
1141:             ++count;
1142:           }
1143:         if (i == 0) PetscCall(MergeRestrictionWeightToRight2d_Private(weight, 2 * factorx - 1, factory));
1144:         else if (i == Mc) PetscCall(MergeRestrictionWeightToLeft2d_Private(weight, 2 * factorx - 1, factory));

1146:         PetscCall(RemoveZeroWeights_Private((2 * factorx - 1) * factory, colf, weight, &count));
1147:         PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1148:         PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1149:         PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1150:       }

1152:   for (d = 0; d < dof[1]; ++d)
1153:     for (j = yc; j < yc + nc + nExtrayc; ++j)
1154:       for (i = xc; i < xc + mc; ++i) {
1155:         rowc.i   = i;
1156:         rowc.j   = j;
1157:         rowc.c   = d;
1158:         rowc.loc = DMSTAG_DOWN;
1159:         count    = 0;
1160:         for (jj = -(factory - 1); jj <= factory - 1; ++jj)
1161:           for (ii = 0; ii < factorx; ++ii) {
1162:             colf[count].i   = i * factorx + ii;
1163:             colf[count].j   = j * factory + jj;
1164:             colf[count].c   = d;
1165:             colf[count].loc = DMSTAG_DOWN;
1166:             PetscCall(SetRestrictionCoefficientCenter_Private(ii, factorx, &ax));
1167:             PetscCall(SetRestrictionCoefficientVertex_Private(jj, factory, &ay));
1168:             weight[count] = ax * ay;
1169:             ++count;
1170:           }
1171:         if (j == 0) PetscCall(MergeRestrictionWeightToTop2d_Private(weight, factorx, 2 * factory - 1));
1172:         else if (j == Nc) PetscCall(MergeRestrictionWeightToBottom2d_Private(weight, factorx, 2 * factory - 1));

1174:         PetscCall(RemoveZeroWeights_Private(factorx * (2 * factory - 1), colf, weight, &count));
1175:         PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1176:         PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1177:         PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1178:       }

1180:   for (d = 0; d < dof[2]; ++d)
1181:     for (j = yc; j < yc + nc; ++j)
1182:       for (i = xc; i < xc + mc; ++i) {
1183:         rowc.i   = i;
1184:         rowc.j   = j;
1185:         rowc.c   = d;
1186:         rowc.loc = DMSTAG_ELEMENT;
1187:         count    = 0;
1188:         for (jj = 0; jj < factory; ++jj)
1189:           for (ii = 0; ii < factorx; ++ii) {
1190:             colf[count].i   = i * factorx + ii;
1191:             colf[count].j   = j * factory + jj;
1192:             colf[count].c   = d;
1193:             colf[count].loc = DMSTAG_ELEMENT;
1194:             weight[count]   = 1 / (PetscScalar)(factorx * factory);
1195:             ++count;
1196:           }

1198:         PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1199:         PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1200:         PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1201:       }

1203:   PetscCall(PetscFree3(colf, weight, ic));
1204:   PetscFunctionReturn(PETSC_SUCCESS);
1205: }

1207: PETSC_INTERN PetscErrorCode DMStagPopulateRestriction3d_Internal(DM dmc, DM dmf, Mat A)
1208: {
1209:   PetscInt       Mc, Nc, Pc, Mf, Nf, Pf, factorx, factory, factorz, dof[4];
1210:   PetscInt       xc, yc, zc, mc, nc, pc, nExtraxc, nExtrayc, nExtrazc, i, j, k, d, ii, jj, kk, count;
1211:   PetscInt       maxFinePoints, maxOffRankFinePoints;
1212:   DMStagStencil  rowc;
1213:   DMStagStencil *colf;
1214:   PetscScalar    ax, ay, az;
1215:   PetscScalar   *weight;
1216:   PetscInt       ir;
1217:   PetscInt      *ic;
1218:   const PetscInt dim = 3;

1220:   PetscFunctionBegin;
1221:   PetscCall(DMStagGetGlobalSizes(dmc, &Mc, &Nc, &Pc));
1222:   PetscCall(DMStagGetGlobalSizes(dmf, &Mf, &Nf, &Pf));
1223:   factorx = Mf / Mc;
1224:   factory = Nf / Nc;
1225:   factorz = Pf / Pc;
1226:   PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], &dof[3]));

1228:   /* In 2D, each coarse point can receive from up to ((2 * factorx - 1) * (2 * factory - 1) * (2 * factorz - 1)) fine points,
1229:      up to ((2 * factorx - 1) * (2 * factory - 1) * (2 * factorz - 1) - factorx * factory * factorz) of which may be off rank */
1230:   maxFinePoints        = (2 * factorx - 1) * (2 * factory - 1) * (2 * factorz - 1);
1231:   maxOffRankFinePoints = maxFinePoints - factorx * factory * factorz;
1232:   PetscCall(MatSeqAIJSetPreallocation(A, maxFinePoints, NULL));
1233:   PetscCall(MatMPIAIJSetPreallocation(A, maxFinePoints, NULL, maxOffRankFinePoints, NULL));
1234:   PetscCall(PetscMalloc3(maxFinePoints, &colf, maxFinePoints, &weight, maxFinePoints, &ic));

1236:   PetscCall(DMStagGetCorners(dmc, &xc, &yc, &zc, &mc, &nc, &pc, &nExtraxc, &nExtrayc, &nExtrazc));

1238:   for (d = 0; d < dof[0]; ++d)
1239:     for (k = zc; k < zc + pc + nExtrazc; ++k)
1240:       for (j = yc; j < yc + nc + nExtrayc; ++j)
1241:         for (i = xc; i < xc + mc + nExtraxc; ++i) {
1242:           rowc.i   = i;
1243:           rowc.j   = j;
1244:           rowc.k   = k;
1245:           rowc.c   = d;
1246:           rowc.loc = DMSTAG_BACK_DOWN_LEFT;
1247:           count    = 0;
1248:           for (kk = -(factorz - 1); kk <= factorz - 1; ++kk)
1249:             for (jj = -(factory - 1); jj <= factory - 1; ++jj)
1250:               for (ii = -(factorx - 1); ii <= factorx - 1; ++ii) {
1251:                 colf[count].i   = i * factorx + ii;
1252:                 colf[count].j   = j * factory + jj;
1253:                 colf[count].k   = k * factorz + kk;
1254:                 colf[count].c   = d;
1255:                 colf[count].loc = DMSTAG_BACK_DOWN_LEFT;
1256:                 PetscCall(SetRestrictionCoefficientVertex_Private(ii, factorx, &ax));
1257:                 PetscCall(SetRestrictionCoefficientVertex_Private(jj, factory, &ay));
1258:                 PetscCall(SetRestrictionCoefficientVertex_Private(kk, factorz, &az));
1259:                 weight[count] = ax * ay * az;
1260:                 ++count;
1261:               }
1262:           if (i == 0) PetscCall(MergeRestrictionWeightToRight3d_Private(weight, 2 * factorx - 1, 2 * factory - 1, 2 * factorz - 1));
1263:           else if (i == Mc) PetscCall(MergeRestrictionWeightToLeft3d_Private(weight, 2 * factorx - 1, 2 * factory - 1, 2 * factorz - 1));
1264:           if (j == 0) PetscCall(MergeRestrictionWeightToTop3d_Private(weight, 2 * factorx - 1, 2 * factory - 1, 2 * factorz - 1));
1265:           else if (j == Nc) PetscCall(MergeRestrictionWeightToBottom3d_Private(weight, 2 * factorx - 1, 2 * factory - 1, 2 * factorz - 1));
1266:           if (k == 0) PetscCall(MergeRestrictionWeightToFront3d_Private(weight, 2 * factorx - 1, 2 * factory - 1, 2 * factorz - 1));
1267:           else if (k == Pc) PetscCall(MergeRestrictionWeightToBack3d_Private(weight, 2 * factorx - 1, 2 * factory - 1, 2 * factorz - 1));

1269:           PetscCall(RemoveZeroWeights_Private((2 * factorx - 1) * (2 * factory - 1) * (2 * factorz - 1), colf, weight, &count));
1270:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1271:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1272:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1273:         }

1275:   for (d = 0; d < dof[1]; ++d)
1276:     for (k = zc; k < zc + pc; ++k)
1277:       for (j = yc; j < yc + nc + nExtrayc; ++j)
1278:         for (i = xc; i < xc + mc + nExtraxc; ++i) {
1279:           rowc.i   = i;
1280:           rowc.j   = j;
1281:           rowc.k   = k;
1282:           rowc.c   = d;
1283:           rowc.loc = DMSTAG_DOWN_LEFT;
1284:           count    = 0;
1285:           for (kk = 0; kk < factorz; ++kk)
1286:             for (jj = -(factory - 1); jj <= factory - 1; ++jj)
1287:               for (ii = -(factorx - 1); ii <= factorx - 1; ++ii) {
1288:                 colf[count].i   = i * factorx + ii;
1289:                 colf[count].j   = j * factory + jj;
1290:                 colf[count].k   = k * factorz + kk;
1291:                 colf[count].c   = d;
1292:                 colf[count].loc = DMSTAG_DOWN_LEFT;
1293:                 PetscCall(SetRestrictionCoefficientVertex_Private(ii, factorx, &ax));
1294:                 PetscCall(SetRestrictionCoefficientVertex_Private(jj, factory, &ay));
1295:                 PetscCall(SetRestrictionCoefficientCenter_Private(kk, factorz, &az));
1296:                 weight[count] = ax * ay * az;
1297:                 ++count;
1298:               }
1299:           if (i == 0) PetscCall(MergeRestrictionWeightToRight3d_Private(weight, 2 * factorx - 1, 2 * factory - 1, factorz));
1300:           else if (i == Mc) PetscCall(MergeRestrictionWeightToLeft3d_Private(weight, 2 * factorx - 1, 2 * factory - 1, factorz));
1301:           if (j == 0) PetscCall(MergeRestrictionWeightToTop3d_Private(weight, 2 * factorx - 1, 2 * factory - 1, factorz));
1302:           else if (j == Nc) PetscCall(MergeRestrictionWeightToBottom3d_Private(weight, 2 * factorx - 1, 2 * factory - 1, factorz));

1304:           PetscCall(RemoveZeroWeights_Private((2 * factorx - 1) * (2 * factory - 1) * factorz, colf, weight, &count));
1305:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1306:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1307:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1308:         }

1310:   for (d = 0; d < dof[1]; ++d)
1311:     for (k = zc; k < zc + pc + nExtrazc; ++k)
1312:       for (j = yc; j < yc + nc; ++j)
1313:         for (i = xc; i < xc + mc + nExtraxc; ++i) {
1314:           rowc.i   = i;
1315:           rowc.j   = j;
1316:           rowc.k   = k;
1317:           rowc.c   = d;
1318:           rowc.loc = DMSTAG_BACK_LEFT;
1319:           count    = 0;
1320:           for (kk = -(factorz - 1); kk <= factorz - 1; ++kk)
1321:             for (jj = 0; jj < factory; ++jj)
1322:               for (ii = -(factorx - 1); ii <= factorx - 1; ++ii) {
1323:                 colf[count].i   = i * factorx + ii;
1324:                 colf[count].j   = j * factory + jj;
1325:                 colf[count].k   = k * factorz + kk;
1326:                 colf[count].c   = d;
1327:                 colf[count].loc = DMSTAG_BACK_LEFT;
1328:                 PetscCall(SetRestrictionCoefficientVertex_Private(ii, factorx, &ax));
1329:                 PetscCall(SetRestrictionCoefficientCenter_Private(jj, factory, &ay));
1330:                 PetscCall(SetRestrictionCoefficientVertex_Private(kk, factorz, &az));
1331:                 weight[count] = ax * ay * az;
1332:                 ++count;
1333:               }
1334:           if (i == 0) PetscCall(MergeRestrictionWeightToRight3d_Private(weight, 2 * factorx - 1, factory, 2 * factorz - 1));
1335:           else if (i == Mc) PetscCall(MergeRestrictionWeightToLeft3d_Private(weight, 2 * factorx - 1, factory, 2 * factorz - 1));
1336:           if (k == 0) PetscCall(MergeRestrictionWeightToFront3d_Private(weight, 2 * factorx - 1, factory, 2 * factorz - 1));
1337:           else if (k == Pc) PetscCall(MergeRestrictionWeightToBack3d_Private(weight, 2 * factorx - 1, factory, 2 * factorz - 1));
1338:         }

1340:   for (d = 0; d < dof[1]; ++d)
1341:     for (k = zc; k < zc + pc + nExtrazc; ++k)
1342:       for (j = yc; j < yc + nc + nExtrayc; ++j)
1343:         for (i = xc; i < xc + mc; ++i) {
1344:           rowc.i   = i;
1345:           rowc.j   = j;
1346:           rowc.k   = k;
1347:           rowc.c   = d;
1348:           rowc.loc = DMSTAG_BACK_DOWN;
1349:           count    = 0;
1350:           for (kk = -(factorz - 1); kk <= factorz - 1; ++kk)
1351:             for (jj = -(factory - 1); jj <= factory - 1; ++jj)
1352:               for (ii = 0; ii < factorx; ++ii) {
1353:                 colf[count].i   = i * factorx + ii;
1354:                 colf[count].j   = j * factory + jj;
1355:                 colf[count].k   = k * factorz + kk;
1356:                 colf[count].c   = d;
1357:                 colf[count].loc = DMSTAG_BACK_DOWN;
1358:                 PetscCall(SetRestrictionCoefficientCenter_Private(ii, factorx, &ax));
1359:                 PetscCall(SetRestrictionCoefficientVertex_Private(jj, factory, &ay));
1360:                 PetscCall(SetRestrictionCoefficientVertex_Private(kk, factorz, &az));
1361:                 weight[count] = ax * ay * az;
1362:                 ++count;
1363:               }
1364:           if (j == 0) PetscCall(MergeRestrictionWeightToTop3d_Private(weight, factorx, 2 * factory - 1, 2 * factorz - 1));
1365:           else if (j == Nc) PetscCall(MergeRestrictionWeightToBottom3d_Private(weight, factorx, 2 * factory - 1, 2 * factorz - 1));
1366:           if (k == 0) PetscCall(MergeRestrictionWeightToFront3d_Private(weight, factorx, 2 * factory - 1, 2 * factorz - 1));
1367:           else if (k == Pc) PetscCall(MergeRestrictionWeightToBack3d_Private(weight, factorx, 2 * factory - 1, 2 * factorz - 1));

1369:           PetscCall(RemoveZeroWeights_Private(factorx * (2 * factory - 1) * (2 * factorz - 1), colf, weight, &count));
1370:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1371:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1372:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1373:         }

1375:   for (d = 0; d < dof[2]; ++d)
1376:     for (k = zc; k < zc + pc; ++k)
1377:       for (j = yc; j < yc + nc; ++j)
1378:         for (i = xc; i < xc + mc + nExtraxc; ++i) {
1379:           rowc.i   = i;
1380:           rowc.j   = j;
1381:           rowc.k   = k;
1382:           rowc.c   = d;
1383:           rowc.loc = DMSTAG_LEFT;
1384:           count    = 0;
1385:           for (kk = 0; kk < factorz; ++kk)
1386:             for (jj = 0; jj < factory; ++jj)
1387:               for (ii = -(factorx - 1); ii <= factorx - 1; ++ii) {
1388:                 colf[count].i   = i * factorx + ii;
1389:                 colf[count].j   = j * factory + jj;
1390:                 colf[count].k   = k * factorz + kk;
1391:                 colf[count].c   = d;
1392:                 colf[count].loc = DMSTAG_LEFT;
1393:                 PetscCall(SetRestrictionCoefficientVertex_Private(ii, factorx, &ax));
1394:                 PetscCall(SetRestrictionCoefficientCenter_Private(jj, factory, &ay));
1395:                 PetscCall(SetRestrictionCoefficientCenter_Private(kk, factorz, &az));
1396:                 weight[count] = ax * ay * az;
1397:                 ++count;
1398:               }
1399:           if (i == 0) PetscCall(MergeRestrictionWeightToRight3d_Private(weight, 2 * factorx - 1, factory, factorz));
1400:           else if (i == Mc) PetscCall(MergeRestrictionWeightToLeft3d_Private(weight, 2 * factorx - 1, factory, factorz));

1402:           PetscCall(RemoveZeroWeights_Private((2 * factorx - 1) * factory * factorz, colf, weight, &count));
1403:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1404:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1405:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1406:         }

1408:   for (d = 0; d < dof[2]; ++d)
1409:     for (k = zc; k < zc + pc; ++k)
1410:       for (j = yc; j < yc + nc + nExtrayc; ++j)
1411:         for (i = xc; i < xc + mc; ++i) {
1412:           rowc.i   = i;
1413:           rowc.j   = j;
1414:           rowc.k   = k;
1415:           rowc.c   = d;
1416:           rowc.loc = DMSTAG_DOWN;
1417:           count    = 0;
1418:           for (kk = 0; kk < factorz; ++kk)
1419:             for (jj = -(factory - 1); jj <= factory - 1; ++jj)
1420:               for (ii = 0; ii < factorx; ++ii) {
1421:                 colf[count].i   = i * factorx + ii;
1422:                 colf[count].j   = j * factory + jj;
1423:                 colf[count].k   = k * factorz + kk;
1424:                 colf[count].c   = d;
1425:                 colf[count].loc = DMSTAG_DOWN;
1426:                 PetscCall(SetRestrictionCoefficientCenter_Private(ii, factorx, &ax));
1427:                 PetscCall(SetRestrictionCoefficientVertex_Private(jj, factory, &ay));
1428:                 PetscCall(SetRestrictionCoefficientCenter_Private(kk, factorz, &az));
1429:                 weight[count] = ax * ay * az;
1430:                 ++count;
1431:               }
1432:           if (j == 0) PetscCall(MergeRestrictionWeightToTop3d_Private(weight, factorx, 2 * factory - 1, factorz));
1433:           else if (j == Nc) PetscCall(MergeRestrictionWeightToBottom3d_Private(weight, factorx, 2 * factory - 1, factorz));

1435:           PetscCall(RemoveZeroWeights_Private(factorx * (2 * factory - 1) * factorz, colf, weight, &count));
1436:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1437:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1438:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1439:         }

1441:   for (d = 0; d < dof[2]; ++d)
1442:     for (k = zc; k < zc + pc + nExtrazc; ++k)
1443:       for (j = yc; j < yc + nc; ++j)
1444:         for (i = xc; i < xc + mc; ++i) {
1445:           rowc.i   = i;
1446:           rowc.j   = j;
1447:           rowc.k   = k;
1448:           rowc.c   = d;
1449:           rowc.loc = DMSTAG_BACK;
1450:           count    = 0;
1451:           for (kk = -(factorz - 1); kk <= factorz - 1; ++kk)
1452:             for (jj = 0; jj < factory; ++jj)
1453:               for (ii = 0; ii < factorx; ++ii) {
1454:                 colf[count].i   = i * factorx + ii;
1455:                 colf[count].j   = j * factory + jj;
1456:                 colf[count].k   = k * factorz + kk;
1457:                 colf[count].c   = d;
1458:                 colf[count].loc = DMSTAG_BACK;
1459:                 PetscCall(SetRestrictionCoefficientCenter_Private(ii, factorx, &ax));
1460:                 PetscCall(SetRestrictionCoefficientCenter_Private(jj, factory, &ay));
1461:                 PetscCall(SetRestrictionCoefficientVertex_Private(kk, factorz, &az));
1462:                 weight[count] = ax * ay * az;
1463:                 ++count;
1464:               }
1465:           if (k == 0) PetscCall(MergeRestrictionWeightToFront3d_Private(weight, factorx, factory, 2 * factorz - 1));
1466:           else if (k == Pc) PetscCall(MergeRestrictionWeightToBack3d_Private(weight, factorx, factory, 2 * factorz - 1));

1468:           PetscCall(RemoveZeroWeights_Private(factorx * factory * (2 * factorz - 1), colf, weight, &count));
1469:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1470:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1471:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1472:         }

1474:   for (d = 0; d < dof[3]; ++d)
1475:     for (k = zc; k < zc + pc; ++k)
1476:       for (j = yc; j < yc + nc; ++j)
1477:         for (i = xc; i < xc + mc; ++i) {
1478:           rowc.i   = i;
1479:           rowc.j   = j;
1480:           rowc.k   = k;
1481:           rowc.c   = d;
1482:           rowc.loc = DMSTAG_ELEMENT;
1483:           count    = 0;
1484:           for (kk = 0; kk < factorz; ++kk)
1485:             for (jj = 0; jj < factory; ++jj)
1486:               for (ii = 0; ii < factorx; ++ii) {
1487:                 colf[count].i   = i * factorx + ii;
1488:                 colf[count].j   = j * factory + jj;
1489:                 colf[count].k   = k * factorz + kk;
1490:                 colf[count].c   = d;
1491:                 colf[count].loc = DMSTAG_ELEMENT;
1492:                 weight[count]   = 1 / (PetscScalar)(factorx * factory * factorz);
1493:                 ++count;
1494:               }

1496:           PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1497:           PetscCall(DMStagStencilToIndexLocal(dmf, dim, count, colf, ic));
1498:           PetscCall(MatSetValuesLocal(A, 1, &ir, count, ic, weight, INSERT_VALUES));
1499:         }

1501:   PetscCall(PetscFree3(colf, weight, ic));
1502:   PetscFunctionReturn(PETSC_SUCCESS);
1503: }