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