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: }