Actual source code: stag1d.c
1: /*
2: Functions specific to the 1-dimensional implementation of DMStag
3: */
4: #include <petsc/private/dmstagimpl.h>
6: /*@
7: DMStagCreate1d - Create an object to manage data living on the elements and vertices of a parallelized regular 1D grid.
9: Collective
11: Input Parameters:
12: + comm - MPI communicator
13: . bndx - boundary type: `DM_BOUNDARY_NONE`, `DM_BOUNDARY_PERIODIC`, or `DM_BOUNDARY_GHOSTED`
14: . M - global number of elements
15: . dof0 - number of degrees of freedom per vertex/0-cell
16: . dof1 - number of degrees of freedom per element/1-cell
17: . stencilType - ghost/halo region type: `DMSTAG_STENCIL_BOX` or `DMSTAG_STENCIL_NONE`
18: . stencilWidth - width, in elements, of halo/ghost region
19: - lx - array of local sizes, of length equal to the comm size, summing to `M` or `NULL`
21: Output Parameter:
22: . dm - the new `DMSTAG` object
24: Options Database Keys:
25: + -dm_view - calls `DMViewFromOptions()` at the conclusion of `DMSetUp()`
26: . -stag_grid_x <nx> - number of elements in the x direction
27: . -stag_ghost_stencil_width - width of ghost region, in elements
28: - -stag_boundary_type_x <none,ghosted,periodic> - `DMBoundaryType` value
30: Level: beginner
32: Notes:
33: You must call `DMSetUp()` after this call before using the `DM`.
34: If you wish to use the options database (see the keys above) to change values in the `DMSTAG`, you must call
35: `DMSetFromOptions()` after this function but before `DMSetUp()`.
37: .seealso: [](ch_stag), `DMSTAG`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMDestroy()`, `DMView()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMLocalToGlobalBegin()`, `DMDACreate1d()`
38: @*/
39: PetscErrorCode DMStagCreate1d(MPI_Comm comm, DMBoundaryType bndx, PetscInt M, PetscInt dof0, PetscInt dof1, DMStagStencilType stencilType, PetscInt stencilWidth, const PetscInt lx[], DM *dm)
40: {
41: PetscMPIInt size;
43: PetscFunctionBegin;
44: PetscCallMPI(MPI_Comm_size(comm, &size));
45: PetscCall(DMCreate(comm, dm));
46: PetscCall(DMSetDimension(*dm, 1));
47: PetscCall(DMStagInitialize(bndx, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, M, 0, 0, size, 0, 0, dof0, dof1, 0, 0, stencilType, stencilWidth, lx, NULL, NULL, *dm));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: PETSC_INTERN PetscErrorCode DMStagRestrictSimple_1d(DM dmf, Vec xf_local, DM dmc, Vec xc_local)
52: {
53: PetscInt Mf, Mc, factorx, dof[2];
54: PetscInt xc, mc, nExtraxc, i, d;
55: PetscInt ileftf, ielemf, ileftc, ielemc;
56: const PetscScalar **arrf;
57: PetscScalar **arrc;
59: PetscFunctionBegin;
60: PetscCall(DMStagGetGlobalSizes(dmf, &Mf, NULL, NULL));
61: PetscCall(DMStagGetGlobalSizes(dmc, &Mc, NULL, NULL));
62: factorx = Mf / Mc;
63: PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], NULL, NULL));
65: PetscCall(DMStagGetCorners(dmc, &xc, NULL, NULL, &mc, NULL, NULL, &nExtraxc, NULL, NULL));
66: PetscCall(VecZeroEntries(xc_local));
67: PetscCall(DMStagVecGetArray(dmf, xf_local, &arrf));
68: PetscCall(DMStagVecGetArray(dmc, xc_local, &arrc));
69: PetscCall(DMStagGetLocationSlot(dmf, DMSTAG_LEFT, 0, &ileftf));
70: PetscCall(DMStagGetLocationSlot(dmf, DMSTAG_ELEMENT, 0, &ielemf));
71: PetscCall(DMStagGetLocationSlot(dmc, DMSTAG_LEFT, 0, &ileftc));
72: PetscCall(DMStagGetLocationSlot(dmc, DMSTAG_ELEMENT, 0, &ielemc));
74: for (d = 0; d < dof[0]; ++d)
75: for (i = xc; i < xc + mc + nExtraxc; ++i) arrc[i][ileftc + d] = arrf[factorx * i][ileftf + d];
77: for (d = 0; d < dof[1]; ++d)
78: for (i = xc; i < xc + mc; ++i) {
79: if (factorx % 2 == 0) arrc[i][ielemc + d] = 0.5 * (arrf[factorx * i + factorx / 2 - 1][ielemf + d] + arrf[factorx * i + factorx / 2][ielemf + d]);
80: else arrc[i][ielemc + d] = arrf[factorx * i + factorx / 2][ielemf + d];
81: }
83: PetscCall(DMStagVecRestoreArray(dmf, xf_local, &arrf));
84: PetscCall(DMStagVecRestoreArray(dmc, xc_local, &arrc));
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: PETSC_INTERN PetscErrorCode DMStagSetUniformCoordinatesExplicit_1d(DM dm, PetscReal xmin, PetscReal xmax)
89: {
90: DM_Stag *stagCoord;
91: DM dmCoord;
92: Vec coordLocal;
93: PetscReal h, min;
94: PetscScalar **arr;
95: PetscInt start_ghost, n_ghost, s;
96: PetscInt ileft, ielement;
98: PetscFunctionBegin;
99: PetscCall(DMGetCoordinateDM(dm, &dmCoord));
100: stagCoord = (DM_Stag *)dmCoord->data;
101: for (s = 0; s < 2; ++s) {
102: PetscCheck(stagCoord->dof[s] == 0 || stagCoord->dof[s] == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Coordinate DM in 1 dimensions must have 0 or 1 dof on each stratum, but stratum %" PetscInt_FMT " has %" PetscInt_FMT " dof", s,
103: stagCoord->dof[s]);
104: }
105: PetscCall(DMCreateLocalVector(dmCoord, &coordLocal));
107: PetscCall(DMStagVecGetArray(dmCoord, coordLocal, &arr));
108: if (stagCoord->dof[0]) PetscCall(DMStagGetLocationSlot(dmCoord, DMSTAG_LEFT, 0, &ileft));
109: if (stagCoord->dof[1]) PetscCall(DMStagGetLocationSlot(dmCoord, DMSTAG_ELEMENT, 0, &ielement));
110: PetscCall(DMStagGetGhostCorners(dmCoord, &start_ghost, NULL, NULL, &n_ghost, NULL, NULL));
112: min = xmin;
113: h = (xmax - xmin) / stagCoord->N[0];
115: for (PetscInt ind = start_ghost; ind < start_ghost + n_ghost; ++ind) {
116: if (stagCoord->dof[0]) {
117: const PetscReal off = 0.0;
118: arr[ind][ileft] = min + ((PetscReal)ind + off) * h;
119: }
120: if (stagCoord->dof[1]) {
121: const PetscReal off = 0.5;
122: arr[ind][ielement] = min + ((PetscReal)ind + off) * h;
123: }
124: }
125: PetscCall(DMStagVecRestoreArray(dmCoord, coordLocal, &arr));
126: PetscCall(DMSetCoordinatesLocal(dm, coordLocal));
127: PetscCall(VecDestroy(&coordLocal));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /* Helper functions used in DMSetUp_Stag() */
132: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM);
134: PETSC_INTERN PetscErrorCode DMSetUp_Stag_1d(DM dm)
135: {
136: DM_Stag *const stag = (DM_Stag *)dm->data;
137: PetscMPIInt size, rank;
138: MPI_Comm comm;
139: PetscInt j;
141: PetscFunctionBegin;
142: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
143: PetscCallMPI(MPI_Comm_size(comm, &size));
144: PetscCallMPI(MPI_Comm_rank(comm, &rank));
146: /* Check Global size */
147: PetscCheck(stag->N[0] >= 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Global grid size of %" PetscInt_FMT " < 1 specified", stag->N[0]);
149: /* Local sizes */
150: PetscCheck(stag->N[0] >= size, comm, PETSC_ERR_ARG_OUTOFRANGE, "More ranks (%d) than elements (%" PetscInt_FMT ") specified", size, stag->N[0]);
151: if (!stag->l[0]) {
152: /* Divide equally, giving an extra elements to higher ranks */
153: PetscCall(PetscMalloc1(stag->nRanks[0], &stag->l[0]));
154: for (j = 0; j < stag->nRanks[0]; ++j) stag->l[0][j] = stag->N[0] / stag->nRanks[0] + (stag->N[0] % stag->nRanks[0] > j ? 1 : 0);
155: }
156: {
157: PetscInt Nchk = 0;
158: for (j = 0; j < size; ++j) Nchk += stag->l[0][j];
159: PetscCheck(Nchk == stag->N[0], comm, PETSC_ERR_ARG_OUTOFRANGE, "Sum of specified local sizes (%" PetscInt_FMT ") is not equal to global size (%" PetscInt_FMT ")", Nchk, stag->N[0]);
160: }
161: stag->n[0] = stag->l[0][rank];
163: /* Rank (trivial in 1d) */
164: stag->rank[0] = rank;
165: stag->firstRank[0] = (PetscBool)(rank == 0);
166: stag->lastRank[0] = (PetscBool)(rank == size - 1);
168: /* Local (unghosted) numbers of entries */
169: stag->entriesPerElement = stag->dof[0] + stag->dof[1];
170: switch (stag->boundaryType[0]) {
171: case DM_BOUNDARY_NONE:
172: case DM_BOUNDARY_GHOSTED:
173: stag->entries = stag->n[0] * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0);
174: break;
175: case DM_BOUNDARY_PERIODIC:
176: stag->entries = stag->n[0] * stag->entriesPerElement;
177: break;
178: default:
179: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
180: }
182: /* Starting element */
183: stag->start[0] = 0;
184: for (j = 0; j < stag->rank[0]; ++j) stag->start[0] += stag->l[0][j];
186: /* Local/ghosted size and starting element */
187: switch (stag->boundaryType[0]) {
188: case DM_BOUNDARY_NONE:
189: switch (stag->stencilType) {
190: case DMSTAG_STENCIL_NONE: /* Only dummy cells on the right */
191: stag->startGhost[0] = stag->start[0];
192: stag->nGhost[0] = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
193: break;
194: case DMSTAG_STENCIL_STAR:
195: case DMSTAG_STENCIL_BOX:
196: stag->startGhost[0] = stag->firstRank[0] ? stag->start[0] : stag->start[0] - stag->stencilWidth;
197: stag->nGhost[0] = stag->n[0];
198: stag->nGhost[0] += stag->firstRank[0] ? 0 : stag->stencilWidth;
199: stag->nGhost[0] += stag->lastRank[0] ? 1 : stag->stencilWidth;
200: break;
201: default:
202: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
203: }
204: break;
205: case DM_BOUNDARY_GHOSTED:
206: switch (stag->stencilType) {
207: case DMSTAG_STENCIL_NONE:
208: stag->startGhost[0] = stag->start[0];
209: stag->nGhost[0] = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
210: break;
211: case DMSTAG_STENCIL_STAR:
212: case DMSTAG_STENCIL_BOX:
213: stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
214: stag->nGhost[0] = stag->n[0] + 2 * stag->stencilWidth + (stag->lastRank[0] && stag->stencilWidth == 0 ? 1 : 0);
215: break;
216: default:
217: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
218: }
219: break;
220: case DM_BOUNDARY_PERIODIC:
221: switch (stag->stencilType) {
222: case DMSTAG_STENCIL_NONE:
223: stag->startGhost[0] = stag->start[0];
224: stag->nGhost[0] = stag->n[0];
225: break;
226: case DMSTAG_STENCIL_STAR:
227: case DMSTAG_STENCIL_BOX:
228: stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
229: stag->nGhost[0] = stag->n[0] + 2 * stag->stencilWidth;
230: break;
231: default:
232: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
233: }
234: break;
235: default:
236: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
237: }
239: /* Total size of ghosted/local representation */
240: stag->entriesGhost = stag->nGhost[0] * stag->entriesPerElement;
242: /* Define neighbors */
243: PetscCall(PetscMalloc1(3, &stag->neighbors));
244: if (stag->firstRank[0]) {
245: switch (stag->boundaryType[0]) {
246: case DM_BOUNDARY_GHOSTED:
247: case DM_BOUNDARY_NONE:
248: stag->neighbors[0] = -1;
249: break;
250: case DM_BOUNDARY_PERIODIC:
251: stag->neighbors[0] = stag->nRanks[0] - 1;
252: break;
253: default:
254: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
255: }
256: } else {
257: stag->neighbors[0] = stag->rank[0] - 1;
258: }
259: stag->neighbors[1] = stag->rank[0];
260: if (stag->lastRank[0]) {
261: switch (stag->boundaryType[0]) {
262: case DM_BOUNDARY_GHOSTED:
263: case DM_BOUNDARY_NONE:
264: stag->neighbors[2] = -1;
265: break;
266: case DM_BOUNDARY_PERIODIC:
267: stag->neighbors[2] = 0;
268: break;
269: default:
270: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
271: }
272: } else {
273: stag->neighbors[2] = stag->rank[0] + 1;
274: }
276: PetscCheck(stag->n[0] >= stag->stencilWidth, PETSC_COMM_SELF, PETSC_ERR_SUP, "DMStag 1d setup does not support local sizes (%" PetscInt_FMT ") smaller than the elementwise stencil width (%" PetscInt_FMT ")", stag->n[0], stag->stencilWidth);
278: /* Create global->local VecScatter and ISLocalToGlobalMapping */
279: {
280: PetscInt *idxLocal, *idxGlobal, *idxGlobalAll;
281: PetscInt i, iLocal, d, entriesToTransferTotal, ghostOffsetStart, ghostOffsetEnd, nNonDummyGhost;
282: IS isLocal, isGlobal;
284: /* The offset on the right (may not be equal to the stencil width, as we
285: always have at least one ghost element, to account for the boundary
286: point, and may with ghosted boundaries), and the number of non-dummy ghost elements */
287: ghostOffsetStart = stag->start[0] - stag->startGhost[0];
288: ghostOffsetEnd = stag->startGhost[0] + stag->nGhost[0] - (stag->start[0] + stag->n[0]);
289: nNonDummyGhost = stag->nGhost[0] - (stag->lastRank[0] ? ghostOffsetEnd : 0) - (stag->firstRank[0] ? ghostOffsetStart : 0);
291: /* Compute the number of non-dummy entries in the local representation
292: This is equal to the number of non-dummy elements in the local (ghosted) representation,
293: plus some extra entries on the right boundary on the last rank*/
294: switch (stag->boundaryType[0]) {
295: case DM_BOUNDARY_GHOSTED:
296: case DM_BOUNDARY_NONE:
297: entriesToTransferTotal = nNonDummyGhost * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0);
298: break;
299: case DM_BOUNDARY_PERIODIC:
300: entriesToTransferTotal = stag->entriesGhost; /* No dummy points */
301: break;
302: default:
303: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
304: }
306: PetscCall(PetscMalloc1(entriesToTransferTotal, &idxLocal));
307: PetscCall(PetscMalloc1(entriesToTransferTotal, &idxGlobal));
308: PetscCall(PetscMalloc1(stag->entriesGhost, &idxGlobalAll));
309: if (stag->boundaryType[0] == DM_BOUNDARY_NONE) {
310: PetscInt count = 0, countAll = 0;
311: /* Left ghost points and native points */
312: for (i = stag->startGhost[0], iLocal = 0; iLocal < nNonDummyGhost; ++i, ++iLocal) {
313: for (d = 0; d < stag->entriesPerElement; ++d, ++count, ++countAll) {
314: idxLocal[count] = iLocal * stag->entriesPerElement + d;
315: idxGlobal[count] = i * stag->entriesPerElement + d;
316: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
317: }
318: }
319: /* Ghost points on the right
320: Special case for last (partial dummy) element on the last rank */
321: if (stag->lastRank[0]) {
322: i = stag->N[0];
323: iLocal = (stag->nGhost[0] - ghostOffsetEnd);
324: /* Only vertex (0-cell) dofs in global representation */
325: for (d = 0; d < stag->dof[0]; ++d, ++count, ++countAll) {
326: idxGlobal[count] = i * stag->entriesPerElement + d;
327: idxLocal[count] = iLocal * stag->entriesPerElement + d;
328: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
329: }
330: for (d = stag->dof[0]; d < stag->entriesPerElement; ++d, ++countAll) { /* Additional dummy entries */
331: idxGlobalAll[countAll] = -1;
332: }
333: }
334: } else if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC) {
335: PetscInt count = 0, iLocal = 0; /* No dummy points, so idxGlobal and idxGlobalAll are identical */
336: const PetscInt iMin = stag->firstRank[0] ? stag->start[0] : stag->startGhost[0];
337: const PetscInt iMax = stag->lastRank[0] ? stag->startGhost[0] + stag->nGhost[0] - stag->stencilWidth : stag->startGhost[0] + stag->nGhost[0];
338: /* Ghost points on the left */
339: if (stag->firstRank[0]) {
340: for (i = stag->N[0] - stag->stencilWidth; iLocal < stag->stencilWidth; ++i, ++iLocal) {
341: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
342: idxGlobal[count] = i * stag->entriesPerElement + d;
343: idxLocal[count] = iLocal * stag->entriesPerElement + d;
344: idxGlobalAll[count] = idxGlobal[count];
345: }
346: }
347: }
348: /* Native points */
349: for (i = iMin; i < iMax; ++i, ++iLocal) {
350: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
351: idxGlobal[count] = i * stag->entriesPerElement + d;
352: idxLocal[count] = iLocal * stag->entriesPerElement + d;
353: idxGlobalAll[count] = idxGlobal[count];
354: }
355: }
356: /* Ghost points on the right */
357: if (stag->lastRank[0]) {
358: for (i = 0; iLocal < stag->nGhost[0]; ++i, ++iLocal) {
359: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
360: idxGlobal[count] = i * stag->entriesPerElement + d;
361: idxLocal[count] = iLocal * stag->entriesPerElement + d;
362: idxGlobalAll[count] = idxGlobal[count];
363: }
364: }
365: }
366: } else if (stag->boundaryType[0] == DM_BOUNDARY_GHOSTED) {
367: PetscInt count = 0, countAll = 0;
368: /* Dummy elements on the left, on the first rank */
369: if (stag->firstRank[0]) {
370: for (iLocal = 0; iLocal < ghostOffsetStart; ++iLocal) {
371: /* Complete elements full of dummy entries */
372: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
373: }
374: i = 0; /* nonDummy entries start with global entry 0 */
375: } else {
376: /* nonDummy entries start as usual */
377: i = stag->startGhost[0];
378: iLocal = 0;
379: }
381: /* non-Dummy entries */
382: {
383: PetscInt iLocalNonDummyMax = stag->firstRank[0] ? nNonDummyGhost + ghostOffsetStart : nNonDummyGhost;
384: for (; iLocal < iLocalNonDummyMax; ++i, ++iLocal) {
385: for (d = 0; d < stag->entriesPerElement; ++d, ++count, ++countAll) {
386: idxLocal[count] = iLocal * stag->entriesPerElement + d;
387: idxGlobal[count] = i * stag->entriesPerElement + d;
388: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
389: }
390: }
391: }
393: /* (partial) dummy elements on the right, on the last rank */
394: if (stag->lastRank[0]) {
395: /* First one is partial dummy */
396: i = stag->N[0];
397: iLocal = (stag->nGhost[0] - ghostOffsetEnd);
398: for (d = 0; d < stag->dof[0]; ++d, ++count, ++countAll) { /* Only vertex (0-cell) dofs in global representation */
399: idxLocal[count] = iLocal * stag->entriesPerElement + d;
400: idxGlobal[count] = i * stag->entriesPerElement + d;
401: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
402: }
403: for (d = stag->dof[0]; d < stag->entriesPerElement; ++d, ++countAll) { /* Additional dummy entries */
404: idxGlobalAll[countAll] = -1;
405: }
406: for (iLocal = stag->nGhost[0] - ghostOffsetEnd + 1; iLocal < stag->nGhost[0]; ++iLocal) {
407: /* Additional dummy elements */
408: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
409: }
410: }
411: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
413: /* Create Local IS (transferring pointer ownership) */
414: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), entriesToTransferTotal, idxLocal, PETSC_OWN_POINTER, &isLocal));
416: /* Create Global IS (transferring pointer ownership) */
417: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), entriesToTransferTotal, idxGlobal, PETSC_OWN_POINTER, &isGlobal));
419: /* Create stag->gtol, which doesn't include dummy entries */
420: {
421: Vec local, global;
422: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm), 1, stag->entries, PETSC_DECIDE, NULL, &global));
423: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, stag->entriesPerElement, stag->entriesGhost, NULL, &local));
424: PetscCall(VecScatterCreate(global, isGlobal, local, isLocal, &stag->gtol));
425: PetscCall(VecDestroy(&global));
426: PetscCall(VecDestroy(&local));
427: }
429: /* In special cases, create a dedicated injective local-to-global map */
430: if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) PetscCall(DMStagPopulateLocalToGlobalInjective(dm));
432: /* Destroy ISs */
433: PetscCall(ISDestroy(&isLocal));
434: PetscCall(ISDestroy(&isGlobal));
436: /* Create local-to-global map (transferring pointer ownership) */
437: PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stag->entriesGhost, idxGlobalAll, PETSC_OWN_POINTER, &dm->ltogmap));
438: }
440: /* Precompute location offsets */
441: PetscCall(DMStagComputeLocationOffsets_1d(dm));
443: /* View from Options */
444: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM dm)
449: {
450: DM_Stag *const stag = (DM_Stag *)dm->data;
451: const PetscInt epe = stag->entriesPerElement;
453: PetscFunctionBegin;
454: PetscCall(PetscMalloc1(DMSTAG_NUMBER_LOCATIONS, &stag->locationOffsets));
455: stag->locationOffsets[DMSTAG_LEFT] = 0;
456: stag->locationOffsets[DMSTAG_ELEMENT] = stag->locationOffsets[DMSTAG_LEFT] + stag->dof[0];
457: stag->locationOffsets[DMSTAG_RIGHT] = stag->locationOffsets[DMSTAG_LEFT] + epe;
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_1d(DM dm)
462: {
463: DM_Stag *const stag = (DM_Stag *)dm->data;
464: PetscInt *idxLocal, *idxGlobal;
465: PetscInt i, iLocal, d, count;
466: IS isLocal, isGlobal;
468: PetscFunctionBegin;
469: PetscCall(PetscMalloc1(stag->entries, &idxLocal));
470: PetscCall(PetscMalloc1(stag->entries, &idxGlobal));
471: count = 0;
472: iLocal = stag->start[0] - stag->startGhost[0];
473: for (i = stag->start[0]; i < stag->start[0] + stag->n[0]; ++i, ++iLocal) {
474: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
475: idxGlobal[count] = i * stag->entriesPerElement + d;
476: idxLocal[count] = iLocal * stag->entriesPerElement + d;
477: }
478: }
479: if (stag->lastRank[0] && stag->boundaryType[0] != DM_BOUNDARY_PERIODIC) {
480: i = stag->start[0] + stag->n[0];
481: iLocal = stag->start[0] - stag->startGhost[0] + stag->n[0];
482: for (d = 0; d < stag->dof[0]; ++d, ++count) {
483: idxGlobal[count] = i * stag->entriesPerElement + d;
484: idxLocal[count] = iLocal * stag->entriesPerElement + d;
485: }
486: }
487: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), stag->entries, idxLocal, PETSC_OWN_POINTER, &isLocal));
488: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), stag->entries, idxGlobal, PETSC_OWN_POINTER, &isGlobal));
489: {
490: Vec local, global;
491: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm), 1, stag->entries, PETSC_DECIDE, NULL, &global));
492: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, stag->entriesPerElement, stag->entriesGhost, NULL, &local));
493: PetscCall(VecScatterCreate(local, isLocal, global, isGlobal, &stag->ltog_injective));
494: PetscCall(VecDestroy(&global));
495: PetscCall(VecDestroy(&local));
496: }
497: PetscCall(ISDestroy(&isLocal));
498: PetscCall(ISDestroy(&isGlobal));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToLocal1d_Internal(DM dm)
503: {
504: DM_Stag *const stag = (DM_Stag *)dm->data;
505: PetscInt *idxRemap;
506: PetscInt i, leftGhostEntries;
508: PetscFunctionBegin;
509: PetscCall(VecScatterCopy(stag->gtol, &stag->ltol));
510: PetscCall(PetscMalloc1(stag->entries, &idxRemap));
512: leftGhostEntries = (stag->start[0] - stag->startGhost[0]) * stag->entriesPerElement;
513: for (i = 0; i < stag->entries; ++i) idxRemap[i] = leftGhostEntries + i;
515: PetscCall(VecScatterRemap(stag->ltol, idxRemap, NULL));
516: PetscCall(PetscFree(idxRemap));
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: PETSC_INTERN PetscErrorCode DMCreateMatrix_Stag_1D_AIJ_Assemble(DM dm, Mat A)
521: {
522: DMStagStencilType stencil_type;
523: PetscInt dof[2], start, n, n_extra, stencil_width, N, epe;
524: DMBoundaryType boundary_type_x;
526: PetscFunctionBegin;
527: PetscCall(DMStagGetDOF(dm, &dof[0], &dof[1], NULL, NULL));
528: PetscCall(DMStagGetStencilType(dm, &stencil_type));
529: PetscCall(DMStagGetStencilWidth(dm, &stencil_width));
530: PetscCall(DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, &n_extra, NULL, NULL));
531: PetscCall(DMStagGetGlobalSizes(dm, &N, NULL, NULL));
532: PetscCall(DMStagGetEntriesPerElement(dm, &epe));
533: PetscCall(DMStagGetBoundaryTypes(dm, &boundary_type_x, NULL, NULL));
534: if (stencil_type == DMSTAG_STENCIL_NONE) {
535: /* Couple all DOF at each location to each other */
536: DMStagStencil *row_vertex, *row_element;
538: PetscCall(PetscMalloc1(dof[0], &row_vertex));
539: for (PetscInt c = 0; c < dof[0]; ++c) {
540: row_vertex[c].loc = DMSTAG_LEFT;
541: row_vertex[c].c = c;
542: }
544: PetscCall(PetscMalloc1(dof[1], &row_element));
545: for (PetscInt c = 0; c < dof[1]; ++c) {
546: row_element[c].loc = DMSTAG_ELEMENT;
547: row_element[c].c = c;
548: }
550: for (PetscInt e = start; e < start + n + n_extra; ++e) {
551: {
552: for (PetscInt c = 0; c < dof[0]; ++c) row_vertex[c].i = e;
553: PetscCall(DMStagMatSetValuesStencil(dm, A, dof[0], row_vertex, dof[0], row_vertex, NULL, INSERT_VALUES));
554: }
555: if (e < N) {
556: for (PetscInt c = 0; c < dof[1]; ++c) row_element[c].i = e;
557: PetscCall(DMStagMatSetValuesStencil(dm, A, dof[1], row_element, dof[1], row_element, NULL, INSERT_VALUES));
558: }
559: }
560: PetscCall(PetscFree(row_vertex));
561: PetscCall(PetscFree(row_element));
562: } else if (stencil_type == DMSTAG_STENCIL_STAR || stencil_type == DMSTAG_STENCIL_BOX) {
563: DMStagStencil *col, *row;
565: PetscCall(PetscMalloc1(epe, &row));
566: {
567: PetscInt nrows = 0;
568: for (PetscInt c = 0; c < dof[0]; ++c) {
569: row[nrows].c = c;
570: row[nrows].loc = DMSTAG_LEFT;
571: ++nrows;
572: }
573: for (PetscInt c = 0; c < dof[1]; ++c) {
574: row[nrows].c = c;
575: row[nrows].loc = DMSTAG_ELEMENT;
576: ++nrows;
577: }
578: }
579: PetscCall(PetscMalloc1(epe, &col));
580: {
581: PetscInt ncols = 0;
582: for (PetscInt c = 0; c < dof[0]; ++c) {
583: col[ncols].c = c;
584: col[ncols].loc = DMSTAG_LEFT;
585: ++ncols;
586: }
587: for (PetscInt c = 0; c < dof[1]; ++c) {
588: col[ncols].c = c;
589: col[ncols].loc = DMSTAG_ELEMENT;
590: ++ncols;
591: }
592: }
593: for (PetscInt e = start; e < start + n + n_extra; ++e) {
594: for (PetscInt i = 0; i < epe; ++i) row[i].i = e;
595: for (PetscInt offset = -stencil_width; offset <= stencil_width; ++offset) {
596: const PetscInt e_offset = e + offset;
598: /* Only set values corresponding to elements which can have non-dummy entries,
599: meaning those that map to unknowns in the global representation. In the periodic
600: case, this is the entire stencil, but in all other cases, only includes a single
601: "extra" element which is partially outside the physical domain (those points in the
602: global representation */
603: if (boundary_type_x == DM_BOUNDARY_PERIODIC || (e_offset < N + 1 && e_offset >= 0)) {
604: for (PetscInt i = 0; i < epe; ++i) col[i].i = e_offset;
605: PetscCall(DMStagMatSetValuesStencil(dm, A, epe, row, epe, col, NULL, INSERT_VALUES));
606: }
607: }
608: }
609: PetscCall(PetscFree(row));
610: PetscCall(PetscFree(col));
611: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported stencil type %s", DMStagStencilTypes[stencil_type]);
612: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
613: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }