Actual source code: stagmulti.c
1: /* Internal and DMStag-specific functions related to multigrid */
2: #include <petsc/private/dmstagimpl.h>
4: /*@C
5: DMStagRestrictSimple - restricts data from a fine to a coarse `DMSTAG`, in the simplest way
7: Values on coarse cells are averages of all fine cells that they cover.
8: Thus, values on vertices are injected, values on edges are averages
9: of the underlying two fine edges, and and values on elements in
10: d dimensions are averages of $2^d$ underlying elements.
12: Input Parameters:
13: + dmf - fine `DM`
14: . xf - data on fine `DM`
15: - dmc - coarse `DM`
17: Output Parameter:
18: . xc - data on coarse `DM`
20: Level: advanced
22: .seealso: [](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: switch (dim) {
31: case 1:
32: PetscCall(DMStagRestrictSimple_1d(dmf, xf, dmc, xc));
33: break;
34: case 2:
35: PetscCall(DMStagRestrictSimple_2d(dmf, xf, dmc, xc));
36: break;
37: case 3:
38: PetscCall(DMStagRestrictSimple_3d(dmf, xf, dmc, xc));
39: break;
40: default:
41: SETERRQ(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT "", dim);
42: break;
43: }
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: /* Code duplication note: the next two functions are nearly identical, save the inclusion of the element terms */
48: PETSC_INTERN PetscErrorCode DMStagPopulateInterpolation1d_a_b_Private(DM dmc, DM dmf, Mat A)
49: {
50: PetscInt exf, startexf, nexf, nextraxf, startexc;
51: PetscInt dof[2];
52: const PetscInt dim = 1;
54: PetscFunctionBegin;
55: PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], NULL, NULL));
56: PetscCheck(dof[0] == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented except for one dof per vertex");
57: PetscCheck(dof[1] <= 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented for more than one dof per element");
59: /* In 1D, each fine point can receive data from at most 2 coarse points, at most one of which could be off-process */
60: PetscCall(MatSeqAIJSetPreallocation(A, 2, NULL));
61: PetscCall(MatMPIAIJSetPreallocation(A, 2, NULL, 1, NULL));
63: PetscCall(DMStagGetCorners(dmf, &startexf, NULL, NULL, &nexf, NULL, NULL, &nextraxf, NULL, NULL));
64: PetscCall(DMStagGetCorners(dmc, &startexc, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
65: for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
66: PetscInt exc, exf_local;
67: exf_local = exf - startexf;
68: exc = startexc + exf_local / 2;
69: /* "even" vertices are just injected, odd vertices averaged */
70: if (exf_local % 2 == 0) {
71: DMStagStencil rowf, colc;
72: PetscInt ir, ic;
73: const PetscScalar one = 1.0;
75: rowf.i = exf;
76: rowf.c = 0;
77: rowf.loc = DMSTAG_LEFT;
78: colc.i = exc;
79: colc.c = 0;
80: colc.loc = DMSTAG_LEFT;
81: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
82: PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &colc, &ic));
83: PetscCall(MatSetValuesLocal(A, 1, &ir, 1, &ic, &one, INSERT_VALUES));
84: } else {
85: DMStagStencil rowf, colc[2];
86: PetscInt ir, ic[2];
87: const PetscScalar weight[2] = {0.5, 0.5};
89: rowf.i = exf;
90: rowf.c = 0;
91: rowf.loc = DMSTAG_LEFT;
92: colc[0].i = exc;
93: colc[0].c = 0;
94: colc[0].loc = DMSTAG_LEFT;
95: colc[1].i = exc;
96: colc[1].c = 0;
97: colc[1].loc = DMSTAG_RIGHT;
98: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
99: PetscCall(DMStagStencilToIndexLocal(dmc, dim, 2, colc, ic));
100: PetscCall(MatSetValuesLocal(A, 1, &ir, 2, ic, weight, INSERT_VALUES));
101: }
102: /* Elements (excluding "extra" dummies) */
103: if (dof[1] > 0 && exf < startexf + nexf) {
104: DMStagStencil rowf, colc;
105: PetscInt ir, ic;
106: const PetscScalar weight = 1.0;
108: rowf.i = exf;
109: rowf.c = 0;
110: rowf.loc = DMSTAG_ELEMENT; /* Note that this assumes only 1 dof */
111: colc.i = exc;
112: colc.c = 0;
113: colc.loc = DMSTAG_ELEMENT;
114: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
115: PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &colc, &ic));
116: PetscCall(MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES));
117: }
118: }
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: PETSC_INTERN PetscErrorCode DMStagPopulateInterpolation2d_0_a_b_Private(DM dmc, DM dmf, Mat A)
123: {
124: PetscInt exf, eyf, startexf, starteyf, nexf, neyf, nextraxf, nextrayf, startexc, starteyc, Nexf, Neyf;
125: PetscInt dof[3];
126: const PetscInt dim = 2;
128: PetscFunctionBegin;
129: PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], NULL));
130: PetscCheck(dof[1] == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented except for one dof per face");
131: PetscCheck(dof[2] <= 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented for more than one dof per element");
133: /* In 2D, each fine point can receive data from at most 4 coarse points, at most 3 of which could be off-process */
134: PetscCall(MatSeqAIJSetPreallocation(A, 4, NULL));
135: PetscCall(MatMPIAIJSetPreallocation(A, 4, NULL, 3, NULL));
137: PetscCall(DMStagGetCorners(dmf, &startexf, &starteyf, NULL, &nexf, &neyf, NULL, &nextraxf, &nextrayf, NULL));
138: PetscCall(DMStagGetCorners(dmc, &startexc, &starteyc, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
139: PetscCall(DMStagGetGlobalSizes(dmf, &Nexf, &Neyf, NULL));
140: for (eyf = starteyf; eyf < starteyf + neyf + nextrayf; ++eyf) {
141: PetscInt eyc, eyf_local;
143: eyf_local = eyf - starteyf;
144: eyc = starteyc + eyf_local / 2;
145: for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
146: PetscInt exc, exf_local;
148: exf_local = exf - startexf;
149: exc = startexc + exf_local / 2;
150: /* Left edges (excluding top "extra" dummy row) */
151: if (eyf < starteyf + neyf) {
152: DMStagStencil rowf, colc[4];
153: PetscInt ir, ic[4], nweight;
154: PetscScalar weight[4];
156: rowf.i = exf;
157: rowf.j = eyf;
158: rowf.c = 0;
159: rowf.loc = DMSTAG_LEFT;
160: colc[0].i = exc;
161: colc[0].j = eyc;
162: colc[0].c = 0;
163: colc[0].loc = DMSTAG_LEFT;
164: if (exf_local % 2 == 0) {
165: if (eyf == Neyf - 1 || eyf == 0) {
166: /* Note - this presumes something like a Neumann condition, assuming
167: a ghost edge with the same value as the adjacent physical edge*/
168: nweight = 1;
169: weight[0] = 1.0;
170: } else {
171: nweight = 2;
172: weight[0] = 0.75;
173: weight[1] = 0.25;
174: if (eyf_local % 2 == 0) {
175: colc[1].i = exc;
176: colc[1].j = eyc - 1;
177: colc[1].c = 0;
178: colc[1].loc = DMSTAG_LEFT;
179: } else {
180: colc[1].i = exc;
181: colc[1].j = eyc + 1;
182: colc[1].c = 0;
183: colc[1].loc = DMSTAG_LEFT;
184: }
185: }
186: } else {
187: colc[1].i = exc;
188: colc[1].j = eyc;
189: colc[1].c = 0;
190: colc[1].loc = DMSTAG_RIGHT;
191: if (eyf == Neyf - 1 || eyf == 0) {
192: /* Note - this presumes something like a Neumann condition, assuming
193: a ghost edge with the same value as the adjacent physical edge*/
194: nweight = 2;
195: weight[0] = 0.5;
196: weight[1] = 0.5;
197: } else {
198: nweight = 4;
199: weight[0] = 0.375;
200: weight[1] = 0.375;
201: weight[2] = 0.125;
202: weight[3] = 0.125;
203: if (eyf_local % 2 == 0) {
204: colc[2].i = exc;
205: colc[2].j = eyc - 1;
206: colc[2].c = 0;
207: colc[2].loc = DMSTAG_LEFT;
208: colc[3].i = exc;
209: colc[3].j = eyc - 1;
210: colc[3].c = 0;
211: colc[3].loc = DMSTAG_RIGHT;
212: } else {
213: colc[2].i = exc;
214: colc[2].j = eyc + 1;
215: colc[2].c = 0;
216: colc[2].loc = DMSTAG_LEFT;
217: colc[3].i = exc;
218: colc[3].j = eyc + 1;
219: colc[3].c = 0;
220: colc[3].loc = DMSTAG_RIGHT;
221: }
222: }
223: }
224: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
225: PetscCall(DMStagStencilToIndexLocal(dmc, dim, nweight, colc, ic));
226: PetscCall(MatSetValuesLocal(A, 1, &ir, nweight, ic, weight, INSERT_VALUES));
227: }
228: /* Down edges (excluding right "extra" dummy col) */
229: if (exf < startexf + nexf) {
230: DMStagStencil rowf, colc[4];
231: PetscInt ir, ic[4], nweight;
232: PetscScalar weight[4];
234: rowf.i = exf;
235: rowf.j = eyf;
236: rowf.c = 0;
237: rowf.loc = DMSTAG_DOWN;
238: colc[0].i = exc;
239: colc[0].j = eyc;
240: colc[0].c = 0;
241: colc[0].loc = DMSTAG_DOWN;
242: if (eyf_local % 2 == 0) {
243: if (exf == Nexf - 1 || exf == 0) {
244: /* Note - this presumes something like a Neumann condition, assuming
245: a ghost edge with the same value as the adjacent physical edge*/
246: nweight = 1;
247: weight[0] = 1.0;
248: } else {
249: nweight = 2;
250: weight[0] = 0.75;
251: weight[1] = 0.25;
252: if (exf_local % 2 == 0) {
253: colc[1].i = exc - 1;
254: colc[1].j = eyc;
255: colc[1].c = 0;
256: colc[1].loc = DMSTAG_DOWN;
257: } else {
258: colc[1].i = exc + 1;
259: colc[1].j = eyc;
260: colc[1].c = 0;
261: colc[1].loc = DMSTAG_DOWN;
262: }
263: }
264: } else {
265: colc[1].i = exc;
266: colc[1].j = eyc;
267: colc[1].c = 0;
268: colc[1].loc = DMSTAG_UP;
269: if (exf == Nexf - 1 || exf == 0) {
270: /* Note - this presumes something like a Neumann condition, assuming
271: a ghost edge with the same value as the adjacent physical edge*/
272: nweight = 2;
273: weight[0] = 0.5;
274: weight[1] = 0.5;
275: } else {
276: nweight = 4;
277: weight[0] = 0.375;
278: weight[1] = 0.375;
279: weight[2] = 0.125;
280: weight[3] = 0.125;
281: if (exf_local % 2 == 0) {
282: colc[2].i = exc - 1;
283: colc[2].j = eyc;
284: colc[2].c = 0;
285: colc[2].loc = DMSTAG_DOWN;
286: colc[3].i = exc - 1;
287: colc[3].j = eyc;
288: colc[3].c = 0;
289: colc[3].loc = DMSTAG_UP;
290: } else {
291: colc[2].i = exc + 1;
292: colc[2].j = eyc;
293: colc[2].c = 0;
294: colc[2].loc = DMSTAG_DOWN;
295: colc[3].i = exc + 1;
296: colc[3].j = eyc;
297: colc[3].c = 0;
298: colc[3].loc = DMSTAG_UP;
299: }
300: }
301: }
302: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
303: PetscCall(DMStagStencilToIndexLocal(dmc, dim, nweight, colc, ic));
304: PetscCall(MatSetValuesLocal(A, 1, &ir, nweight, ic, weight, INSERT_VALUES));
305: }
306: /* Elements (excluding "extra" dummy) */
307: if (dof[2] > 0 && exf < startexf + nexf && eyf < starteyf + neyf) {
308: DMStagStencil rowf, colc;
309: PetscInt ir, ic;
310: const PetscScalar weight = 1.0;
312: rowf.i = exf;
313: rowf.j = eyf;
314: rowf.c = 0;
315: rowf.loc = DMSTAG_ELEMENT;
316: colc.i = exc;
317: colc.j = eyc;
318: colc.c = 0;
319: colc.loc = DMSTAG_ELEMENT;
320: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
321: PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &colc, &ic));
322: PetscCall(MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES));
323: }
324: }
325: }
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: PETSC_INTERN PetscErrorCode DMStagPopulateInterpolation3d_0_0_a_b_Private(DM dmc, DM dmf, Mat A)
330: {
331: PetscInt exf, eyf, ezf, startexf, starteyf, startezf, nexf, neyf, nezf, nextraxf, nextrayf, nextrazf, startexc, starteyc, startezc, Nexf, Neyf, Nezf;
332: PetscInt dof[4];
333: const PetscInt dim = 3;
335: PetscFunctionBegin;
337: PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], &dof[3]));
338: PetscCheck(dof[2] == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented except for one dof per face");
339: PetscCheck(dof[3] <= 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented for more than one dof per element");
341: /* In 3D, each fine point can receive data from at most 8 coarse points, at most 7 of which could be off-process */
342: PetscCall(MatSeqAIJSetPreallocation(A, 8, NULL));
343: PetscCall(MatMPIAIJSetPreallocation(A, 8, NULL, 7, NULL));
345: PetscCall(DMStagGetCorners(dmf, &startexf, &starteyf, &startezf, &nexf, &neyf, &nezf, &nextraxf, &nextrayf, &nextrazf));
346: PetscCall(DMStagGetCorners(dmc, &startexc, &starteyc, &startezc, NULL, NULL, NULL, NULL, NULL, NULL));
347: PetscCall(DMStagGetGlobalSizes(dmf, &Nexf, &Neyf, &Nezf));
348: for (ezf = startezf; ezf < startezf + nezf + nextrayf; ++ezf) {
349: const PetscInt ezf_local = ezf - startezf;
350: const PetscInt ezc = startezc + ezf_local / 2;
351: const PetscBool boundary_z = (PetscBool)(ezf == 0 || ezf == Nezf - 1);
353: for (eyf = starteyf; eyf < starteyf + neyf + nextrayf; ++eyf) {
354: const PetscInt eyf_local = eyf - starteyf;
355: const PetscInt eyc = starteyc + eyf_local / 2;
356: const PetscBool boundary_y = (PetscBool)(eyf == 0 || eyf == Neyf - 1);
358: for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
359: const PetscInt exf_local = exf - startexf;
360: const PetscInt exc = startexc + exf_local / 2;
361: const PetscBool boundary_x = (PetscBool)(exf == 0 || exf == Nexf - 1);
363: /* Left faces (excluding top and front "extra" dummy layers) */
364: if (eyf < starteyf + neyf && ezf < startezf + nezf) {
365: DMStagStencil rowf, colc[8];
366: PetscInt ir, ic[8], nweight;
367: PetscScalar weight[8];
369: rowf.i = exf;
370: rowf.j = eyf;
371: rowf.k = ezf;
372: rowf.c = 0;
373: rowf.loc = DMSTAG_LEFT;
374: colc[0].i = exc;
375: colc[0].j = eyc;
376: colc[0].k = ezc;
377: colc[0].c = 0;
378: colc[0].loc = DMSTAG_LEFT;
379: if (exf_local % 2 == 0) {
380: if (boundary_y) {
381: if (boundary_z) {
382: nweight = 1;
383: weight[0] = 1.0;
384: } else {
385: const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;
387: nweight = 2;
388: weight[0] = 0.75;
389: weight[1] = 0.25;
390: colc[1].i = exc;
391: colc[1].j = eyc;
392: colc[1].k = ezc + ezc_offset;
393: colc[1].c = 0;
394: colc[1].loc = DMSTAG_LEFT;
395: }
396: } else {
397: const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;
399: if (boundary_z) {
400: nweight = 2;
401: weight[0] = 0.75;
402: weight[1] = 0.25;
403: colc[1].i = exc;
404: colc[1].j = eyc + eyc_offset;
405: colc[1].k = ezc;
406: colc[1].c = 0;
407: colc[1].loc = DMSTAG_LEFT;
408: } else {
409: const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;
411: nweight = 4;
412: weight[0] = 0.75 * 0.75;
413: weight[1] = weight[2] = 0.25 * 0.75;
414: weight[3] = 0.25 * 0.25;
415: colc[1].i = exc;
416: colc[1].j = eyc;
417: colc[1].k = ezc + ezc_offset;
418: colc[1].c = 0;
419: colc[1].loc = DMSTAG_LEFT;
420: colc[2].i = exc;
421: colc[2].j = eyc + eyc_offset;
422: colc[2].k = ezc;
423: colc[2].c = 0;
424: colc[2].loc = DMSTAG_LEFT;
425: colc[3].i = exc;
426: colc[3].j = eyc + eyc_offset;
427: colc[3].k = ezc + ezc_offset;
428: colc[3].c = 0;
429: colc[3].loc = DMSTAG_LEFT;
430: }
431: }
432: } else {
433: colc[1].i = exc;
434: colc[1].j = eyc;
435: colc[1].k = ezc;
436: colc[1].c = 0;
437: colc[1].loc = DMSTAG_RIGHT;
438: if (boundary_y) {
439: if (boundary_z) {
440: nweight = 2;
441: weight[0] = weight[1] = 0.5;
442: } else {
443: const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;
445: nweight = 4;
446: weight[0] = weight[1] = 0.5 * 0.75;
447: weight[2] = weight[3] = 0.5 * 0.25;
448: colc[2].i = exc;
449: colc[2].j = eyc;
450: colc[2].k = ezc + ezc_offset;
451: colc[2].c = 0;
452: colc[2].loc = DMSTAG_LEFT;
453: colc[3].i = exc;
454: colc[3].j = eyc;
455: colc[3].k = ezc + ezc_offset;
456: colc[3].c = 0;
457: colc[3].loc = DMSTAG_RIGHT;
458: }
459: } else {
460: const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;
462: if (boundary_z) {
463: nweight = 4;
464: weight[0] = weight[1] = 0.5 * 0.75;
465: weight[2] = weight[3] = 0.5 * 0.25;
466: colc[2].i = exc;
467: colc[2].j = eyc + eyc_offset;
468: colc[2].k = ezc;
469: colc[2].c = 0;
470: colc[2].loc = DMSTAG_LEFT;
471: colc[3].i = exc;
472: colc[3].j = eyc + eyc_offset;
473: colc[3].k = ezc;
474: colc[3].c = 0;
475: colc[3].loc = DMSTAG_RIGHT;
476: } else {
477: const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;
479: nweight = 8;
480: weight[0] = weight[1] = 0.5 * 0.75 * 0.75;
481: weight[2] = weight[3] = weight[4] = weight[5] = 0.5 * 0.25 * 0.75;
482: weight[6] = weight[7] = 0.5 * 0.25 * 0.25;
483: colc[2].i = exc;
484: colc[2].j = eyc;
485: colc[2].k = ezc + ezc_offset;
486: colc[2].c = 0;
487: colc[2].loc = DMSTAG_LEFT;
488: colc[3].i = exc;
489: colc[3].j = eyc;
490: colc[3].k = ezc + ezc_offset;
491: colc[3].c = 0;
492: colc[3].loc = DMSTAG_RIGHT;
493: colc[4].i = exc;
494: colc[4].j = eyc + eyc_offset;
495: colc[4].k = ezc;
496: colc[4].c = 0;
497: colc[4].loc = DMSTAG_LEFT;
498: colc[5].i = exc;
499: colc[5].j = eyc + eyc_offset;
500: colc[5].k = ezc;
501: colc[5].c = 0;
502: colc[5].loc = DMSTAG_RIGHT;
503: colc[6].i = exc;
504: colc[6].j = eyc + eyc_offset;
505: colc[6].k = ezc + ezc_offset;
506: colc[6].c = 0;
507: colc[6].loc = DMSTAG_LEFT;
508: colc[7].i = exc;
509: colc[7].j = eyc + eyc_offset;
510: colc[7].k = ezc + ezc_offset;
511: colc[7].c = 0;
512: colc[7].loc = DMSTAG_RIGHT;
513: }
514: }
515: }
516: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
517: PetscCall(DMStagStencilToIndexLocal(dmc, dim, nweight, colc, ic));
518: PetscCall(MatSetValuesLocal(A, 1, &ir, nweight, ic, weight, INSERT_VALUES));
519: }
521: /* Bottom faces (excluding left and front "extra" dummy layers) */
522: if (exf < startexf + nexf && ezf < startezf + nezf) {
523: DMStagStencil rowf, colc[8];
524: PetscInt ir, ic[8], nweight;
525: PetscScalar weight[8];
527: rowf.i = exf;
528: rowf.j = eyf;
529: rowf.k = ezf;
530: rowf.c = 0;
531: rowf.loc = DMSTAG_DOWN;
532: colc[0].i = exc;
533: colc[0].j = eyc;
534: colc[0].k = ezc;
535: colc[0].c = 0;
536: colc[0].loc = DMSTAG_DOWN;
537: if (eyf_local % 2 == 0) {
538: if (boundary_x) {
539: if (boundary_z) {
540: nweight = 1;
541: weight[0] = 1.0;
542: } else {
543: const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;
545: nweight = 2;
546: weight[0] = 0.75;
547: weight[1] = 0.25;
548: colc[1].i = exc;
549: colc[1].j = eyc;
550: colc[1].k = ezc + ezc_offset;
551: colc[1].c = 0;
552: colc[1].loc = DMSTAG_DOWN;
553: }
554: } else {
555: const PetscInt exc_offset = exf_local % 2 == 0 ? -1 : 1;
557: if (boundary_z) {
558: nweight = 2;
559: weight[0] = 0.75;
560: weight[1] = 0.25;
561: colc[1].i = exc + exc_offset;
562: colc[1].j = eyc;
563: colc[1].k = ezc;
564: colc[1].c = 0;
565: colc[1].loc = DMSTAG_DOWN;
566: } else {
567: const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;
569: nweight = 4;
570: weight[0] = 0.75 * 0.75;
571: weight[1] = weight[2] = 0.25 * 0.75;
572: weight[3] = 0.25 * 0.25;
573: colc[1].i = exc;
574: colc[1].j = eyc;
575: colc[1].k = ezc + ezc_offset;
576: colc[1].c = 0;
577: colc[1].loc = DMSTAG_DOWN;
578: colc[2].i = exc + exc_offset;
579: colc[2].j = eyc;
580: colc[2].k = ezc;
581: colc[2].c = 0;
582: colc[2].loc = DMSTAG_DOWN;
583: colc[3].i = exc + exc_offset;
584: colc[3].j = eyc;
585: colc[3].k = ezc + ezc_offset;
586: colc[3].c = 0;
587: colc[3].loc = DMSTAG_DOWN;
588: }
589: }
590: } else {
591: colc[1].i = exc;
592: colc[1].j = eyc;
593: colc[1].k = ezc;
594: colc[1].c = 0;
595: colc[1].loc = DMSTAG_UP;
596: if (boundary_x) {
597: if (boundary_z) {
598: nweight = 2;
599: weight[0] = weight[1] = 0.5;
600: } else {
601: const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;
603: nweight = 4;
604: weight[0] = weight[1] = 0.5 * 0.75;
605: weight[2] = weight[3] = 0.5 * 0.25;
606: colc[2].i = exc;
607: colc[2].j = eyc;
608: colc[2].k = ezc + ezc_offset;
609: colc[2].c = 0;
610: colc[2].loc = DMSTAG_DOWN;
611: colc[3].i = exc;
612: colc[3].j = eyc;
613: colc[3].k = ezc + ezc_offset;
614: colc[3].c = 0;
615: colc[3].loc = DMSTAG_UP;
616: }
617: } else {
618: const PetscInt exc_offset = exf_local % 2 == 0 ? -1 : 1;
620: if (boundary_z) {
621: nweight = 4;
622: weight[0] = weight[1] = 0.5 * 0.75;
623: weight[2] = weight[3] = 0.5 * 0.25;
624: colc[2].i = exc + exc_offset;
625: colc[2].j = eyc;
626: colc[2].k = ezc;
627: colc[2].c = 0;
628: colc[2].loc = DMSTAG_DOWN;
629: colc[3].i = exc + exc_offset;
630: colc[3].j = eyc;
631: colc[3].k = ezc;
632: colc[3].c = 0;
633: colc[3].loc = DMSTAG_UP;
634: } else {
635: const PetscInt ezc_offset = ezf_local % 2 == 0 ? -1 : 1;
637: nweight = 8;
638: weight[0] = weight[1] = 0.5 * 0.75 * 0.75;
639: weight[2] = weight[3] = weight[4] = weight[5] = 0.5 * 0.25 * 0.75;
640: weight[6] = weight[7] = 0.5 * 0.25 * 0.25;
641: colc[2].i = exc;
642: colc[2].j = eyc;
643: colc[2].k = ezc + ezc_offset;
644: colc[2].c = 0;
645: colc[2].loc = DMSTAG_DOWN;
646: colc[3].i = exc;
647: colc[3].j = eyc;
648: colc[3].k = ezc + ezc_offset;
649: colc[3].c = 0;
650: colc[3].loc = DMSTAG_UP;
651: colc[4].i = exc + exc_offset;
652: colc[4].j = eyc;
653: colc[4].k = ezc;
654: colc[4].c = 0;
655: colc[4].loc = DMSTAG_DOWN;
656: colc[5].i = exc + exc_offset;
657: colc[5].j = eyc;
658: colc[5].k = ezc;
659: colc[5].c = 0;
660: colc[5].loc = DMSTAG_UP;
661: colc[6].i = exc + exc_offset;
662: colc[6].j = eyc;
663: colc[6].k = ezc + ezc_offset;
664: colc[6].c = 0;
665: colc[6].loc = DMSTAG_DOWN;
666: colc[7].i = exc + exc_offset;
667: colc[7].j = eyc;
668: colc[7].k = ezc + ezc_offset;
669: colc[7].c = 0;
670: colc[7].loc = DMSTAG_UP;
671: }
672: }
673: }
674: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
675: PetscCall(DMStagStencilToIndexLocal(dmc, dim, nweight, colc, ic));
676: PetscCall(MatSetValuesLocal(A, 1, &ir, nweight, ic, weight, INSERT_VALUES));
677: }
679: /* Back faces (excluding left and top "extra" dummy layers) */
680: if (exf < startexf + nexf && ezf < startezf + nezf) {
681: DMStagStencil rowf, colc[8];
682: PetscInt ir, ic[8], nweight;
683: PetscScalar weight[8];
685: rowf.i = exf;
686: rowf.j = eyf;
687: rowf.k = ezf;
688: rowf.c = 0;
689: rowf.loc = DMSTAG_BACK;
690: colc[0].i = exc;
691: colc[0].j = eyc;
692: colc[0].k = ezc;
693: colc[0].c = 0;
694: colc[0].loc = DMSTAG_BACK;
695: if (ezf_local % 2 == 0) {
696: if (boundary_x) {
697: if (boundary_y) {
698: nweight = 1;
699: weight[0] = 1.0;
700: } else {
701: const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;
703: nweight = 2;
704: weight[0] = 0.75;
705: weight[1] = 0.25;
706: colc[1].i = exc;
707: colc[1].j = eyc + eyc_offset;
708: colc[1].k = ezc;
709: colc[1].c = 0;
710: colc[1].loc = DMSTAG_BACK;
711: }
712: } else {
713: const PetscInt exc_offset = exf_local % 2 == 0 ? -1 : 1;
715: if (boundary_y) {
716: nweight = 2;
717: weight[0] = 0.75;
718: weight[1] = 0.25;
719: colc[1].i = exc + exc_offset;
720: colc[1].j = eyc;
721: colc[1].k = ezc;
722: colc[1].c = 0;
723: colc[1].loc = DMSTAG_BACK;
724: } else {
725: const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;
727: nweight = 4;
728: weight[0] = 0.75 * 0.75;
729: weight[1] = weight[2] = 0.25 * 0.75;
730: weight[3] = 0.25 * 0.25;
731: colc[1].i = exc + exc_offset;
732: colc[1].j = eyc;
733: colc[1].k = ezc;
734: colc[1].c = 0;
735: colc[1].loc = DMSTAG_BACK;
736: colc[2].i = exc;
737: colc[2].j = eyc + eyc_offset;
738: colc[2].k = ezc;
739: colc[2].c = 0;
740: colc[2].loc = DMSTAG_BACK;
741: colc[3].i = exc + exc_offset;
742: colc[3].j = eyc + eyc_offset;
743: colc[3].k = ezc;
744: colc[3].c = 0;
745: colc[3].loc = DMSTAG_BACK;
746: }
747: }
748: } else {
749: colc[1].i = exc;
750: colc[1].j = eyc;
751: colc[1].k = ezc;
752: colc[1].c = 0;
753: colc[1].loc = DMSTAG_FRONT;
754: if (boundary_x) {
755: if (boundary_y) {
756: nweight = 2;
757: weight[0] = weight[1] = 0.5;
758: } else {
759: const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;
761: nweight = 4;
762: weight[0] = weight[1] = 0.5 * 0.75;
763: weight[2] = weight[3] = 0.5 * 0.25;
764: colc[2].i = exc;
765: colc[2].j = eyc + eyc_offset;
766: colc[2].k = ezc;
767: colc[2].c = 0;
768: colc[2].loc = DMSTAG_BACK;
769: colc[3].i = exc;
770: colc[3].j = eyc + eyc_offset;
771: colc[3].k = ezc;
772: colc[3].c = 0;
773: colc[3].loc = DMSTAG_FRONT;
774: }
775: } else {
776: const PetscInt exc_offset = exf_local % 2 == 0 ? -1 : 1;
778: if (boundary_y) {
779: nweight = 4;
780: weight[0] = weight[1] = 0.5 * 0.75;
781: weight[2] = weight[3] = 0.5 * 0.25;
782: colc[2].i = exc + exc_offset;
783: colc[2].j = eyc;
784: colc[2].k = ezc;
785: colc[2].c = 0;
786: colc[2].loc = DMSTAG_BACK;
787: colc[3].i = exc + exc_offset;
788: colc[3].j = eyc;
789: colc[3].k = ezc;
790: colc[3].c = 0;
791: colc[3].loc = DMSTAG_FRONT;
792: } else {
793: const PetscInt eyc_offset = eyf_local % 2 == 0 ? -1 : 1;
795: nweight = 8;
796: weight[0] = weight[1] = 0.5 * 0.75 * 0.75;
797: weight[2] = weight[3] = weight[4] = weight[5] = 0.5 * 0.25 * 0.75;
798: weight[6] = weight[7] = 0.5 * 0.25 * 0.25;
799: colc[2].i = exc;
800: colc[2].j = eyc + eyc_offset;
801: colc[2].k = ezc;
802: colc[2].c = 0;
803: colc[2].loc = DMSTAG_BACK;
804: colc[3].i = exc;
805: colc[3].j = eyc + eyc_offset;
806: colc[3].k = ezc;
807: colc[3].c = 0;
808: colc[3].loc = DMSTAG_FRONT;
809: colc[4].i = exc + exc_offset;
810: colc[4].j = eyc;
811: colc[4].k = ezc;
812: colc[4].c = 0;
813: colc[4].loc = DMSTAG_BACK;
814: colc[5].i = exc + exc_offset;
815: colc[5].j = eyc;
816: colc[5].k = ezc;
817: colc[5].c = 0;
818: colc[5].loc = DMSTAG_FRONT;
819: colc[6].i = exc + exc_offset;
820: colc[6].j = eyc + eyc_offset;
821: colc[6].k = ezc;
822: colc[6].c = 0;
823: colc[6].loc = DMSTAG_BACK;
824: colc[7].i = exc + exc_offset;
825: colc[7].j = eyc + eyc_offset;
826: colc[7].k = ezc;
827: colc[7].c = 0;
828: colc[7].loc = DMSTAG_FRONT;
829: }
830: }
831: }
832: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
833: PetscCall(DMStagStencilToIndexLocal(dmc, dim, nweight, colc, ic));
834: PetscCall(MatSetValuesLocal(A, 1, &ir, nweight, ic, weight, INSERT_VALUES));
835: }
836: /* Elements */
837: if (dof[3] == 1 && exf < startexf + nexf && eyf < starteyf + neyf && ezf < startezf + nezf) {
838: DMStagStencil rowf, colc;
839: PetscInt ir, ic;
840: const PetscScalar weight = 1.0;
842: rowf.i = exf;
843: rowf.j = eyf;
844: rowf.k = ezf;
845: rowf.c = 0;
846: rowf.loc = DMSTAG_ELEMENT;
847: colc.i = exc;
848: colc.j = eyc;
849: colc.k = ezc;
850: colc.c = 0;
851: colc.loc = DMSTAG_ELEMENT;
852: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &rowf, &ir));
853: PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &colc, &ic));
854: PetscCall(MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES));
855: }
856: }
857: }
858: }
859: PetscFunctionReturn(PETSC_SUCCESS);
860: }
862: PETSC_INTERN PetscErrorCode DMStagPopulateRestriction1d_a_b_Private(DM dmc, DM dmf, Mat A)
863: {
864: PetscInt exf, startexf, nexf, nextraxf, startexc, Nexf;
865: PetscInt dof[2];
866: const PetscInt dim = 1;
868: PetscFunctionBegin;
869: PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], NULL, NULL));
870: PetscCheck(dof[0] == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented except for one dof per vertex");
871: PetscCheck(dof[1] <= 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented for more than one dof per element");
873: /* In 1D, each coarse point can receive from up to 3 fine points, one of which may be off-rank */
874: PetscCall(MatSeqAIJSetPreallocation(A, 3, NULL));
875: PetscCall(MatMPIAIJSetPreallocation(A, 3, NULL, 1, NULL));
877: PetscCall(DMStagGetCorners(dmf, &startexf, NULL, NULL, &nexf, NULL, NULL, &nextraxf, NULL, NULL));
878: PetscCall(DMStagGetCorners(dmc, &startexc, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
879: PetscCall(DMStagGetGlobalSizes(dmf, &Nexf, NULL, NULL));
880: for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
881: PetscInt exc, exf_local;
883: exf_local = exf - startexf;
884: exc = startexc + exf_local / 2;
885: /* "even" vertices contribute to the overlying coarse vertex, odd vertices to the two adjacent */
886: if (exf_local % 2 == 0) {
887: DMStagStencil colf, rowc;
888: PetscInt ir, ic;
889: PetscScalar weight;
891: colf.i = exf;
892: colf.c = 0;
893: colf.loc = DMSTAG_LEFT;
894: rowc.i = exc;
895: rowc.c = 0;
896: rowc.loc = DMSTAG_LEFT;
897: PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
898: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic));
899: weight = (exf == Nexf || exf == 0) ? 0.75 : 0.5; /* Assume a Neumann-type condition */
900: PetscCall(MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES));
901: } else {
902: DMStagStencil colf, rowc[2];
903: PetscInt ic, ir[2];
904: const PetscScalar weight[2] = {0.25, 0.25};
906: colf.i = exf;
907: colf.c = 0;
908: colf.loc = DMSTAG_LEFT;
909: rowc[0].i = exc;
910: rowc[0].c = 0;
911: rowc[0].loc = DMSTAG_LEFT;
912: rowc[1].i = exc;
913: rowc[1].c = 0;
914: rowc[1].loc = DMSTAG_RIGHT;
915: PetscCall(DMStagStencilToIndexLocal(dmc, dim, 2, rowc, ir));
916: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic));
917: PetscCall(MatSetValuesLocal(A, 2, ir, 1, &ic, weight, INSERT_VALUES));
918: }
919: if (dof[1] > 0 && exf < startexf + nexf) {
920: DMStagStencil rowc, colf;
921: PetscInt ir, ic;
922: const PetscScalar weight = 0.5;
924: rowc.i = exc;
925: rowc.c = 0;
926: rowc.loc = DMSTAG_ELEMENT;
927: colf.i = exf;
928: colf.c = 0;
929: colf.loc = DMSTAG_ELEMENT;
930: PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
931: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic));
932: PetscCall(MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES));
933: }
934: }
935: PetscFunctionReturn(PETSC_SUCCESS);
936: }
938: PETSC_INTERN PetscErrorCode DMStagPopulateRestriction2d_0_a_b_Private(DM dmc, DM dmf, Mat A)
939: {
940: PetscInt exf, eyf, startexf, starteyf, nexf, neyf, nextraxf, nextrayf, startexc, starteyc, Nexf, Neyf;
941: PetscInt dof[3];
942: const PetscInt dim = 2;
944: PetscFunctionBegin;
945: PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], NULL));
946: PetscCheck(dof[1] == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented except for one dof per face");
947: PetscCheck(dof[2] <= 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented for more than one dof per element");
949: /* In 2D, each coarse point can receive from up to 6 fine points,
950: up to 2 of which may be off rank */
951: PetscCall(MatSeqAIJSetPreallocation(A, 6, NULL));
952: PetscCall(MatMPIAIJSetPreallocation(A, 6, NULL, 2, NULL));
954: PetscCall(DMStagGetCorners(dmf, &startexf, &starteyf, NULL, &nexf, &neyf, NULL, &nextraxf, &nextrayf, NULL));
955: PetscCall(DMStagGetCorners(dmc, &startexc, &starteyc, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
956: PetscCall(DMStagGetGlobalSizes(dmf, &Nexf, &Neyf, NULL));
957: for (eyf = starteyf; eyf < starteyf + neyf + nextrayf; ++eyf) {
958: PetscInt eyc, eyf_local;
959: eyf_local = eyf - starteyf;
960: eyc = starteyc + eyf_local / 2;
961: for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
962: PetscInt exc, exf_local;
963: exf_local = exf - startexf;
964: exc = startexc + exf_local / 2;
965: /* Left edges (excluding top "extra" dummy row) */
966: if (eyf < starteyf + neyf) {
967: DMStagStencil rowc[2], colf;
968: PetscInt ir[2], ic, nweight;
969: PetscScalar weight[2];
971: colf.i = exf;
972: colf.j = eyf;
973: colf.c = 0;
974: colf.loc = DMSTAG_LEFT;
975: rowc[0].i = exc;
976: rowc[0].j = eyc;
977: rowc[0].c = 0;
978: rowc[0].loc = DMSTAG_LEFT;
979: if (exf_local % 2 == 0) {
980: nweight = 1;
981: if (exf == Nexf || exf == 0) {
982: /* Note - this presumes something like a Neumann condition, assuming
983: a ghost edge with the same value as the adjacent physical edge*/
984: weight[0] = 0.375;
985: } else {
986: weight[0] = 0.25;
987: }
988: } else {
989: nweight = 2;
990: rowc[1].i = exc;
991: rowc[1].j = eyc;
992: rowc[1].c = 0;
993: rowc[1].loc = DMSTAG_RIGHT;
994: weight[0] = 0.125;
995: weight[1] = 0.125;
996: }
997: PetscCall(DMStagStencilToIndexLocal(dmc, dim, nweight, rowc, ir));
998: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic));
999: PetscCall(MatSetValuesLocal(A, nweight, ir, 1, &ic, weight, INSERT_VALUES));
1000: }
1001: /* Down edges (excluding right "extra" dummy col) */
1002: if (exf < startexf + nexf) {
1003: DMStagStencil rowc[2], colf;
1004: PetscInt ir[2], ic, nweight;
1005: PetscScalar weight[2];
1007: colf.i = exf;
1008: colf.j = eyf;
1009: colf.c = 0;
1010: colf.loc = DMSTAG_DOWN;
1011: rowc[0].i = exc;
1012: rowc[0].j = eyc;
1013: rowc[0].c = 0;
1014: rowc[0].loc = DMSTAG_DOWN;
1015: if (eyf_local % 2 == 0) {
1016: nweight = 1;
1017: if (eyf == Neyf || eyf == 0) {
1018: /* Note - this presumes something like a Neumann condition, assuming
1019: a ghost edge with the same value as the adjacent physical edge*/
1020: weight[0] = 0.375;
1021: } else {
1022: weight[0] = 0.25;
1023: }
1024: } else {
1025: nweight = 2;
1026: rowc[1].i = exc;
1027: rowc[1].j = eyc;
1028: rowc[1].c = 0;
1029: rowc[1].loc = DMSTAG_UP;
1030: weight[0] = 0.125;
1031: weight[1] = 0.125;
1032: }
1033: PetscCall(DMStagStencilToIndexLocal(dmc, dim, nweight, rowc, ir));
1034: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic));
1035: PetscCall(MatSetValuesLocal(A, nweight, ir, 1, &ic, weight, INSERT_VALUES));
1036: }
1037: /* Elements (excluding "extra" dummies) */
1038: if (dof[2] > 0 && exf < startexf + nexf && eyf < starteyf + neyf) {
1039: DMStagStencil rowc, colf;
1040: PetscInt ir, ic;
1041: const PetscScalar cellScale = 0.25;
1043: rowc.i = exc;
1044: rowc.j = eyc;
1045: rowc.c = 0;
1046: rowc.loc = DMSTAG_ELEMENT;
1047: colf.i = exf;
1048: colf.j = eyf;
1049: colf.c = 0;
1050: colf.loc = DMSTAG_ELEMENT;
1051: PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1052: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic));
1053: PetscCall(MatSetValuesLocal(A, 1, &ir, 1, &ic, &cellScale, INSERT_VALUES));
1054: }
1055: }
1056: }
1057: PetscFunctionReturn(PETSC_SUCCESS);
1058: }
1060: PETSC_INTERN PetscErrorCode DMStagPopulateRestriction3d_0_0_a_b_Private(DM dmc, DM dmf, Mat A)
1061: {
1062: PetscInt exf, eyf, ezf, startexf, starteyf, startezf, nexf, neyf, nezf, nextraxf, nextrayf, nextrazf, startexc, starteyc, startezc, Nexf, Neyf, Nezf;
1063: PetscInt dof[4];
1064: const PetscInt dim = 3;
1066: PetscFunctionBegin;
1068: PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], &dof[3]));
1069: PetscCheck(dof[2] == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented except for one dof per face");
1070: PetscCheck(dof[3] <= 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "Not Implemented for more than one dof per element");
1072: /* In 3D, each coarse point can receive from up to 12 fine points,
1073: up to 4 of which may be off rank */
1074: PetscCall(MatSeqAIJSetPreallocation(A, 12, NULL));
1075: PetscCall(MatMPIAIJSetPreallocation(A, 12, NULL, 4, NULL));
1077: PetscCall(DMStagGetCorners(dmf, &startexf, &starteyf, &startezf, &nexf, &neyf, &nezf, &nextraxf, &nextrayf, &nextrazf));
1078: PetscCall(DMStagGetCorners(dmc, &startexc, &starteyc, &startezc, NULL, NULL, NULL, NULL, NULL, NULL));
1079: PetscCall(DMStagGetGlobalSizes(dmf, &Nexf, &Neyf, &Nezf));
1081: for (ezf = startezf; ezf < startezf + nezf + nextrazf; ++ezf) {
1082: const PetscInt ezf_local = ezf - startezf;
1083: const PetscInt ezc = startezc + ezf_local / 2;
1085: for (eyf = starteyf; eyf < starteyf + neyf + nextrayf; ++eyf) {
1086: const PetscInt eyf_local = eyf - starteyf;
1087: const PetscInt eyc = starteyc + eyf_local / 2;
1089: for (exf = startexf; exf < startexf + nexf + nextraxf; ++exf) {
1090: const PetscInt exf_local = exf - startexf;
1091: const PetscInt exc = startexc + exf_local / 2;
1093: /* Left faces (excluding top and front "extra" dummy layers) */
1094: if (eyf < starteyf + neyf && ezf < startezf + nezf) {
1095: DMStagStencil rowc[2], colf;
1096: PetscInt ir[2], ic, nweight;
1097: PetscScalar weight[2];
1099: colf.i = exf;
1100: colf.j = eyf;
1101: colf.k = ezf;
1102: colf.c = 0;
1103: colf.loc = DMSTAG_LEFT;
1104: rowc[0].i = exc;
1105: rowc[0].j = eyc;
1106: rowc[0].k = ezc;
1107: rowc[0].c = 0;
1108: rowc[0].loc = DMSTAG_LEFT;
1109: if (exf_local % 2 == 0) {
1110: nweight = 1;
1111: if (exf == Nexf || exf == 0) {
1112: /* Note - this presumes something like a Neumann condition, assuming
1113: a ghost edge with the same value as the adjacent physical edge*/
1114: weight[0] = 0.1875;
1115: } else {
1116: weight[0] = 0.125;
1117: }
1118: } else {
1119: nweight = 2;
1120: rowc[1].i = exc;
1121: rowc[1].j = eyc;
1122: rowc[1].k = ezc;
1123: rowc[1].c = 0;
1124: rowc[1].loc = DMSTAG_RIGHT;
1125: weight[0] = 0.0625;
1126: weight[1] = 0.0625;
1127: }
1128: PetscCall(DMStagStencilToIndexLocal(dmc, dim, nweight, rowc, ir));
1129: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic));
1130: PetscCall(MatSetValuesLocal(A, nweight, ir, 1, &ic, weight, INSERT_VALUES));
1131: }
1133: /* Down faces (excluding right and front "extra" dummy layers) */
1134: if (exf < startexf + nexf && ezf < startezf + nezf) {
1135: DMStagStencil rowc[2], colf;
1136: PetscInt ir[2], ic, nweight;
1137: PetscScalar weight[2];
1139: colf.i = exf;
1140: colf.j = eyf;
1141: colf.k = ezf;
1142: colf.c = 0;
1143: colf.loc = DMSTAG_DOWN;
1144: rowc[0].i = exc;
1145: rowc[0].j = eyc;
1146: rowc[0].k = ezc;
1147: rowc[0].c = 0;
1148: rowc[0].loc = DMSTAG_DOWN;
1149: if (eyf_local % 2 == 0) {
1150: nweight = 1;
1151: if (eyf == Neyf || eyf == 0) {
1152: /* Note - this presumes something like a Neumann condition, assuming
1153: a ghost edge with the same value as the adjacent physical edge*/
1154: weight[0] = 0.1875;
1155: } else {
1156: weight[0] = 0.125;
1157: }
1158: } else {
1159: nweight = 2;
1160: rowc[1].i = exc;
1161: rowc[1].j = eyc;
1162: rowc[1].k = ezc;
1163: rowc[1].c = 0;
1164: rowc[1].loc = DMSTAG_UP;
1165: weight[0] = 0.0625;
1166: weight[1] = 0.0625;
1167: }
1168: PetscCall(DMStagStencilToIndexLocal(dmc, dim, nweight, rowc, ir));
1169: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic));
1170: PetscCall(MatSetValuesLocal(A, nweight, ir, 1, &ic, weight, INSERT_VALUES));
1171: }
1173: /* Back faces (excluding left and top "extra" dummy laers) */
1174: if (exf < startexf + nexf && eyf < starteyf + neyf) {
1175: DMStagStencil rowc[2], colf;
1176: PetscInt ir[2], ic, nweight;
1177: PetscScalar weight[2];
1179: colf.i = exf;
1180: colf.j = eyf;
1181: colf.k = ezf;
1182: colf.c = 0;
1183: colf.loc = DMSTAG_BACK;
1184: rowc[0].i = exc;
1185: rowc[0].j = eyc;
1186: rowc[0].k = ezc;
1187: rowc[0].c = 0;
1188: rowc[0].loc = DMSTAG_BACK;
1189: if (ezf_local % 2 == 0) {
1190: nweight = 1;
1191: if (ezf == Nezf || ezf == 0) {
1192: /* Note - this presumes something like a Neumann condition, assuming
1193: a ghost edge with the same value as the adjacent physical edge*/
1194: weight[0] = 0.1875;
1195: } else {
1196: weight[0] = 0.125;
1197: }
1198: } else {
1199: nweight = 2;
1200: rowc[1].i = exc;
1201: rowc[1].j = eyc;
1202: rowc[1].k = ezc;
1203: rowc[1].c = 0;
1204: rowc[1].loc = DMSTAG_FRONT;
1205: weight[0] = 0.0625;
1206: weight[1] = 0.0625;
1207: }
1208: PetscCall(DMStagStencilToIndexLocal(dmc, dim, nweight, rowc, ir));
1209: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic));
1210: PetscCall(MatSetValuesLocal(A, nweight, ir, 1, &ic, weight, INSERT_VALUES));
1211: }
1212: /* Elements */
1213: if (dof[3] == 1 && exf < startexf + nexf && eyf < starteyf + neyf && ezf < startezf + nezf) {
1214: DMStagStencil rowc, colf;
1215: PetscInt ir, ic;
1216: const PetscScalar weight = 0.125;
1218: colf.i = exf;
1219: colf.j = eyf;
1220: colf.k = ezf;
1221: colf.c = 0;
1222: colf.loc = DMSTAG_ELEMENT;
1223: rowc.i = exc;
1224: rowc.j = eyc;
1225: rowc.k = ezc;
1226: rowc.c = 0;
1227: rowc.loc = DMSTAG_ELEMENT;
1228: PetscCall(DMStagStencilToIndexLocal(dmc, dim, 1, &rowc, &ir));
1229: PetscCall(DMStagStencilToIndexLocal(dmf, dim, 1, &colf, &ic));
1230: PetscCall(MatSetValuesLocal(A, 1, &ir, 1, &ic, &weight, INSERT_VALUES));
1231: }
1232: }
1233: }
1234: }
1235: PetscFunctionReturn(PETSC_SUCCESS);
1236: }