Actual source code: stag1d.c
1: /*
2: Functions specific to the 1-dimensional implementation of DMStag
3: */
4: #include <petsc/private/dmstagimpl.h>
6: /*@C
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
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: [](chapter_stag), `DMSTAG`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMDestroy()`, `DMView()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMLocalToGlobalBegin()`, `DMDACreate1d()`
38: @*/
39: PETSC_EXTERN 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: MPI_Comm_size(comm, &size);
44: DMCreate(comm, dm);
45: DMSetDimension(*dm, 1);
46: DMStagInitialize(bndx, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, M, 0, 0, size, 0, 0, dof0, dof1, 0, 0, stencilType, stencilWidth, lx, NULL, NULL, *dm);
47: return 0;
48: }
50: PETSC_INTERN PetscErrorCode DMStagRestrictSimple_1d(DM dmf, Vec xf_local, DM dmc, Vec xc_local)
51: {
52: PetscScalar **LA_xf, **LA_xc;
53: PetscInt i, start, n, nextra, N;
54: PetscInt d, dof[2];
55: PetscInt slot_left_coarse, slot_element_coarse, slot_left_fine, slot_element_fine;
57: DMStagGetDOF(dmc, &dof[0], &dof[1], NULL, NULL);
58: DMStagGetCorners(dmc, &start, NULL, NULL, &n, NULL, NULL, &nextra, NULL, NULL);
59: DMStagGetGlobalSizes(dmc, &N, NULL, NULL);
60: if (PetscDefined(USE_DEBUG)) {
61: PetscInt dof_check[2], n_fine, start_fine;
63: DMStagGetDOF(dmf, &dof_check[0], &dof_check[1], NULL, NULL);
64: DMStagGetCorners(dmf, &start_fine, NULL, NULL, &n_fine, NULL, NULL, NULL, NULL, NULL);
65: for (d = 0; d < 2; ++d)
69: {
70: PetscInt size_local, entries_local;
72: DMStagGetEntriesLocal(dmf, &entries_local);
73: VecGetLocalSize(xf_local, &size_local);
75: }
76: {
77: PetscInt size_local, entries_local;
79: DMStagGetEntriesLocal(dmc, &entries_local);
80: VecGetLocalSize(xc_local, &size_local);
82: }
83: }
84: VecZeroEntries(xc_local);
85: DMStagVecGetArray(dmf, xf_local, &LA_xf);
86: DMStagVecGetArray(dmc, xc_local, &LA_xc);
87: DMStagGetLocationSlot(dmf, DMSTAG_LEFT, 0, &slot_left_fine);
88: DMStagGetLocationSlot(dmf, DMSTAG_ELEMENT, 0, &slot_element_fine);
89: DMStagGetLocationSlot(dmc, DMSTAG_LEFT, 0, &slot_left_coarse);
90: DMStagGetLocationSlot(dmc, DMSTAG_ELEMENT, 0, &slot_element_coarse);
91: for (i = start; i < start + n + nextra; ++i) {
92: for (d = 0; d < dof[0]; ++d) LA_xc[i][slot_left_coarse + d] = LA_xf[2 * i][slot_left_fine + d];
93: if (i < N) {
94: for (d = 0; d < dof[1]; ++d) LA_xc[i][slot_element_coarse + d] = 0.5 * (LA_xf[2 * i][slot_element_fine + d] + LA_xf[2 * i + 1][slot_element_fine + d]);
95: }
96: }
97: DMStagVecRestoreArray(dmf, xf_local, &LA_xf);
98: DMStagVecRestoreArray(dmc, xc_local, &LA_xc);
99: return 0;
100: }
102: PETSC_INTERN PetscErrorCode DMStagSetUniformCoordinatesExplicit_1d(DM dm, PetscReal xmin, PetscReal xmax)
103: {
104: DM_Stag *stagCoord;
105: DM dmCoord;
106: Vec coordLocal;
107: PetscReal h, min;
108: PetscScalar **arr;
109: PetscInt start_ghost, n_ghost, s;
110: PetscInt ileft, ielement;
112: DMGetCoordinateDM(dm, &dmCoord);
113: stagCoord = (DM_Stag *)dmCoord->data;
114: for (s = 0; s < 2; ++s) {
116: stagCoord->dof[s]);
117: }
118: DMCreateLocalVector(dmCoord, &coordLocal);
120: DMStagVecGetArray(dmCoord, coordLocal, &arr);
121: if (stagCoord->dof[0]) DMStagGetLocationSlot(dmCoord, DMSTAG_LEFT, 0, &ileft);
122: if (stagCoord->dof[1]) DMStagGetLocationSlot(dmCoord, DMSTAG_ELEMENT, 0, &ielement);
123: DMStagGetGhostCorners(dmCoord, &start_ghost, NULL, NULL, &n_ghost, NULL, NULL);
125: min = xmin;
126: h = (xmax - xmin) / stagCoord->N[0];
128: for (PetscInt ind = start_ghost; ind < start_ghost + n_ghost; ++ind) {
129: if (stagCoord->dof[0]) {
130: const PetscReal off = 0.0;
131: arr[ind][ileft] = min + ((PetscReal)ind + off) * h;
132: }
133: if (stagCoord->dof[1]) {
134: const PetscReal off = 0.5;
135: arr[ind][ielement] = min + ((PetscReal)ind + off) * h;
136: }
137: }
138: DMStagVecRestoreArray(dmCoord, coordLocal, &arr);
139: DMSetCoordinatesLocal(dm, coordLocal);
140: VecDestroy(&coordLocal);
141: return 0;
142: }
144: /* Helper functions used in DMSetUp_Stag() */
145: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM);
147: PETSC_INTERN PetscErrorCode DMSetUp_Stag_1d(DM dm)
148: {
149: DM_Stag *const stag = (DM_Stag *)dm->data;
150: PetscMPIInt size, rank;
151: MPI_Comm comm;
152: PetscInt j;
154: PetscObjectGetComm((PetscObject)dm, &comm);
155: MPI_Comm_size(comm, &size);
156: MPI_Comm_rank(comm, &rank);
158: /* Check Global size */
161: /* Local sizes */
163: if (!stag->l[0]) {
164: /* Divide equally, giving an extra elements to higher ranks */
165: PetscMalloc1(stag->nRanks[0], &stag->l[0]);
166: 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);
167: }
168: {
169: PetscInt Nchk = 0;
170: for (j = 0; j < size; ++j) Nchk += stag->l[0][j];
172: }
173: stag->n[0] = stag->l[0][rank];
175: /* Rank (trivial in 1d) */
176: stag->rank[0] = rank;
177: stag->firstRank[0] = (PetscBool)(rank == 0);
178: stag->lastRank[0] = (PetscBool)(rank == size - 1);
180: /* Local (unghosted) numbers of entries */
181: stag->entriesPerElement = stag->dof[0] + stag->dof[1];
182: switch (stag->boundaryType[0]) {
183: case DM_BOUNDARY_NONE:
184: case DM_BOUNDARY_GHOSTED:
185: stag->entries = stag->n[0] * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0);
186: break;
187: case DM_BOUNDARY_PERIODIC:
188: stag->entries = stag->n[0] * stag->entriesPerElement;
189: break;
190: default:
191: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
192: }
194: /* Starting element */
195: stag->start[0] = 0;
196: for (j = 0; j < stag->rank[0]; ++j) stag->start[0] += stag->l[0][j];
198: /* Local/ghosted size and starting element */
199: switch (stag->boundaryType[0]) {
200: case DM_BOUNDARY_NONE:
201: switch (stag->stencilType) {
202: case DMSTAG_STENCIL_NONE: /* Only dummy cells on the right */
203: stag->startGhost[0] = stag->start[0];
204: stag->nGhost[0] = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
205: break;
206: case DMSTAG_STENCIL_STAR:
207: case DMSTAG_STENCIL_BOX:
208: stag->startGhost[0] = stag->firstRank[0] ? stag->start[0] : stag->start[0] - stag->stencilWidth;
209: stag->nGhost[0] = stag->n[0];
210: stag->nGhost[0] += stag->firstRank[0] ? 0 : stag->stencilWidth;
211: stag->nGhost[0] += stag->lastRank[0] ? 1 : stag->stencilWidth;
212: break;
213: default:
214: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
215: }
216: break;
217: case DM_BOUNDARY_GHOSTED:
218: switch (stag->stencilType) {
219: case DMSTAG_STENCIL_NONE:
220: stag->startGhost[0] = stag->start[0];
221: stag->nGhost[0] = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
222: break;
223: case DMSTAG_STENCIL_STAR:
224: case DMSTAG_STENCIL_BOX:
225: stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
226: stag->nGhost[0] = stag->n[0] + 2 * stag->stencilWidth + (stag->lastRank[0] && stag->stencilWidth == 0 ? 1 : 0);
227: break;
228: default:
229: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
230: }
231: break;
232: case DM_BOUNDARY_PERIODIC:
233: switch (stag->stencilType) {
234: case DMSTAG_STENCIL_NONE:
235: stag->startGhost[0] = stag->start[0];
236: stag->nGhost[0] = stag->n[0];
237: break;
238: case DMSTAG_STENCIL_STAR:
239: case DMSTAG_STENCIL_BOX:
240: stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
241: stag->nGhost[0] = stag->n[0] + 2 * stag->stencilWidth;
242: break;
243: default:
244: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
245: }
246: break;
247: default:
248: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
249: }
251: /* Total size of ghosted/local representation */
252: stag->entriesGhost = stag->nGhost[0] * stag->entriesPerElement;
254: /* Define neighbors */
255: PetscMalloc1(3, &stag->neighbors);
256: if (stag->firstRank[0]) {
257: switch (stag->boundaryType[0]) {
258: case DM_BOUNDARY_GHOSTED:
259: case DM_BOUNDARY_NONE:
260: stag->neighbors[0] = -1;
261: break;
262: case DM_BOUNDARY_PERIODIC:
263: stag->neighbors[0] = stag->nRanks[0] - 1;
264: break;
265: default:
266: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
267: }
268: } else {
269: stag->neighbors[0] = stag->rank[0] - 1;
270: }
271: stag->neighbors[1] = stag->rank[0];
272: if (stag->lastRank[0]) {
273: switch (stag->boundaryType[0]) {
274: case DM_BOUNDARY_GHOSTED:
275: case DM_BOUNDARY_NONE:
276: stag->neighbors[2] = -1;
277: break;
278: case DM_BOUNDARY_PERIODIC:
279: stag->neighbors[2] = 0;
280: break;
281: default:
282: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
283: }
284: } else {
285: stag->neighbors[2] = stag->rank[0] + 1;
286: }
290: /* Create global->local VecScatter and ISLocalToGlobalMapping */
291: {
292: PetscInt *idxLocal, *idxGlobal, *idxGlobalAll;
293: PetscInt i, iLocal, d, entriesToTransferTotal, ghostOffsetStart, ghostOffsetEnd, nNonDummyGhost;
294: IS isLocal, isGlobal;
296: /* The offset on the right (may not be equal to the stencil width, as we
297: always have at least one ghost element, to account for the boundary
298: point, and may with ghosted boundaries), and the number of non-dummy ghost elements */
299: ghostOffsetStart = stag->start[0] - stag->startGhost[0];
300: ghostOffsetEnd = stag->startGhost[0] + stag->nGhost[0] - (stag->start[0] + stag->n[0]);
301: nNonDummyGhost = stag->nGhost[0] - (stag->lastRank[0] ? ghostOffsetEnd : 0) - (stag->firstRank[0] ? ghostOffsetStart : 0);
303: /* Compute the number of non-dummy entries in the local representation
304: This is equal to the number of non-dummy elements in the local (ghosted) representation,
305: plus some extra entries on the right boundary on the last rank*/
306: switch (stag->boundaryType[0]) {
307: case DM_BOUNDARY_GHOSTED:
308: case DM_BOUNDARY_NONE:
309: entriesToTransferTotal = nNonDummyGhost * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0);
310: break;
311: case DM_BOUNDARY_PERIODIC:
312: entriesToTransferTotal = stag->entriesGhost; /* No dummy points */
313: break;
314: default:
315: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
316: }
318: PetscMalloc1(entriesToTransferTotal, &idxLocal);
319: PetscMalloc1(entriesToTransferTotal, &idxGlobal);
320: PetscMalloc1(stag->entriesGhost, &idxGlobalAll);
321: if (stag->boundaryType[0] == DM_BOUNDARY_NONE) {
322: PetscInt count = 0, countAll = 0;
323: /* Left ghost points and native points */
324: for (i = stag->startGhost[0], iLocal = 0; iLocal < nNonDummyGhost; ++i, ++iLocal) {
325: for (d = 0; d < stag->entriesPerElement; ++d, ++count, ++countAll) {
326: idxLocal[count] = iLocal * stag->entriesPerElement + d;
327: idxGlobal[count] = i * stag->entriesPerElement + d;
328: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
329: }
330: }
331: /* Ghost points on the right
332: Special case for last (partial dummy) element on the last rank */
333: if (stag->lastRank[0]) {
334: i = stag->N[0];
335: iLocal = (stag->nGhost[0] - ghostOffsetEnd);
336: /* Only vertex (0-cell) dofs in global representation */
337: for (d = 0; d < stag->dof[0]; ++d, ++count, ++countAll) {
338: idxGlobal[count] = i * stag->entriesPerElement + d;
339: idxLocal[count] = iLocal * stag->entriesPerElement + d;
340: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
341: }
342: for (d = stag->dof[0]; d < stag->entriesPerElement; ++d, ++countAll) { /* Additional dummy entries */
343: idxGlobalAll[countAll] = -1;
344: }
345: }
346: } else if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC) {
347: PetscInt count = 0, iLocal = 0; /* No dummy points, so idxGlobal and idxGlobalAll are identical */
348: const PetscInt iMin = stag->firstRank[0] ? stag->start[0] : stag->startGhost[0];
349: const PetscInt iMax = stag->lastRank[0] ? stag->startGhost[0] + stag->nGhost[0] - stag->stencilWidth : stag->startGhost[0] + stag->nGhost[0];
350: /* Ghost points on the left */
351: if (stag->firstRank[0]) {
352: for (i = stag->N[0] - stag->stencilWidth; iLocal < stag->stencilWidth; ++i, ++iLocal) {
353: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
354: idxGlobal[count] = i * stag->entriesPerElement + d;
355: idxLocal[count] = iLocal * stag->entriesPerElement + d;
356: idxGlobalAll[count] = idxGlobal[count];
357: }
358: }
359: }
360: /* Native points */
361: for (i = iMin; i < iMax; ++i, ++iLocal) {
362: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
363: idxGlobal[count] = i * stag->entriesPerElement + d;
364: idxLocal[count] = iLocal * stag->entriesPerElement + d;
365: idxGlobalAll[count] = idxGlobal[count];
366: }
367: }
368: /* Ghost points on the right */
369: if (stag->lastRank[0]) {
370: for (i = 0; iLocal < stag->nGhost[0]; ++i, ++iLocal) {
371: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
372: idxGlobal[count] = i * stag->entriesPerElement + d;
373: idxLocal[count] = iLocal * stag->entriesPerElement + d;
374: idxGlobalAll[count] = idxGlobal[count];
375: }
376: }
377: }
378: } else if (stag->boundaryType[0] == DM_BOUNDARY_GHOSTED) {
379: PetscInt count = 0, countAll = 0;
380: /* Dummy elements on the left, on the first rank */
381: if (stag->firstRank[0]) {
382: for (iLocal = 0; iLocal < ghostOffsetStart; ++iLocal) {
383: /* Complete elements full of dummy entries */
384: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
385: }
386: i = 0; /* nonDummy entries start with global entry 0 */
387: } else {
388: /* nonDummy entries start as usual */
389: i = stag->startGhost[0];
390: iLocal = 0;
391: }
393: /* non-Dummy entries */
394: {
395: PetscInt iLocalNonDummyMax = stag->firstRank[0] ? nNonDummyGhost + ghostOffsetStart : nNonDummyGhost;
396: for (; iLocal < iLocalNonDummyMax; ++i, ++iLocal) {
397: for (d = 0; d < stag->entriesPerElement; ++d, ++count, ++countAll) {
398: idxLocal[count] = iLocal * stag->entriesPerElement + d;
399: idxGlobal[count] = i * stag->entriesPerElement + d;
400: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
401: }
402: }
403: }
405: /* (partial) dummy elements on the right, on the last rank */
406: if (stag->lastRank[0]) {
407: /* First one is partial dummy */
408: i = stag->N[0];
409: iLocal = (stag->nGhost[0] - ghostOffsetEnd);
410: for (d = 0; d < stag->dof[0]; ++d, ++count, ++countAll) { /* Only vertex (0-cell) dofs in global representation */
411: idxLocal[count] = iLocal * stag->entriesPerElement + d;
412: idxGlobal[count] = i * stag->entriesPerElement + d;
413: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
414: }
415: for (d = stag->dof[0]; d < stag->entriesPerElement; ++d, ++countAll) { /* Additional dummy entries */
416: idxGlobalAll[countAll] = -1;
417: }
418: for (iLocal = stag->nGhost[0] - ghostOffsetEnd + 1; iLocal < stag->nGhost[0]; ++iLocal) {
419: /* Additional dummy elements */
420: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
421: }
422: }
423: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
425: /* Create Local IS (transferring pointer ownership) */
426: ISCreateGeneral(PetscObjectComm((PetscObject)dm), entriesToTransferTotal, idxLocal, PETSC_OWN_POINTER, &isLocal);
428: /* Create Global IS (transferring pointer ownership) */
429: ISCreateGeneral(PetscObjectComm((PetscObject)dm), entriesToTransferTotal, idxGlobal, PETSC_OWN_POINTER, &isGlobal);
431: /* Create stag->gtol, which doesn't include dummy entries */
432: {
433: Vec local, global;
434: VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm), 1, stag->entries, PETSC_DECIDE, NULL, &global);
435: VecCreateSeqWithArray(PETSC_COMM_SELF, stag->entriesPerElement, stag->entriesGhost, NULL, &local);
436: VecScatterCreate(global, isGlobal, local, isLocal, &stag->gtol);
437: VecDestroy(&global);
438: VecDestroy(&local);
439: }
441: /* In special cases, create a dedicated injective local-to-global map */
442: if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) DMStagPopulateLocalToGlobalInjective(dm);
444: /* Destroy ISs */
445: ISDestroy(&isLocal);
446: ISDestroy(&isGlobal);
448: /* Create local-to-global map (transferring pointer ownership) */
449: ISLocalToGlobalMappingCreate(comm, 1, stag->entriesGhost, idxGlobalAll, PETSC_OWN_POINTER, &dm->ltogmap);
450: }
452: /* Precompute location offsets */
453: DMStagComputeLocationOffsets_1d(dm);
455: /* View from Options */
456: DMViewFromOptions(dm, NULL, "-dm_view");
458: return 0;
459: }
461: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM dm)
462: {
463: DM_Stag *const stag = (DM_Stag *)dm->data;
464: const PetscInt epe = stag->entriesPerElement;
466: PetscMalloc1(DMSTAG_NUMBER_LOCATIONS, &stag->locationOffsets);
467: stag->locationOffsets[DMSTAG_LEFT] = 0;
468: stag->locationOffsets[DMSTAG_ELEMENT] = stag->locationOffsets[DMSTAG_LEFT] + stag->dof[0];
469: stag->locationOffsets[DMSTAG_RIGHT] = stag->locationOffsets[DMSTAG_LEFT] + epe;
470: return 0;
471: }
473: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_1d(DM dm)
474: {
475: DM_Stag *const stag = (DM_Stag *)dm->data;
476: PetscInt *idxLocal, *idxGlobal;
477: PetscInt i, iLocal, d, count;
478: IS isLocal, isGlobal;
480: PetscMalloc1(stag->entries, &idxLocal);
481: PetscMalloc1(stag->entries, &idxGlobal);
482: count = 0;
483: iLocal = stag->start[0] - stag->startGhost[0];
484: for (i = stag->start[0]; i < stag->start[0] + stag->n[0]; ++i, ++iLocal) {
485: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
486: idxGlobal[count] = i * stag->entriesPerElement + d;
487: idxLocal[count] = iLocal * stag->entriesPerElement + d;
488: }
489: }
490: if (stag->lastRank[0] && stag->boundaryType[0] != DM_BOUNDARY_PERIODIC) {
491: i = stag->start[0] + stag->n[0];
492: iLocal = stag->start[0] - stag->startGhost[0] + stag->n[0];
493: for (d = 0; d < stag->dof[0]; ++d, ++count) {
494: idxGlobal[count] = i * stag->entriesPerElement + d;
495: idxLocal[count] = iLocal * stag->entriesPerElement + d;
496: }
497: }
498: ISCreateGeneral(PetscObjectComm((PetscObject)dm), stag->entries, idxLocal, PETSC_OWN_POINTER, &isLocal);
499: ISCreateGeneral(PetscObjectComm((PetscObject)dm), stag->entries, idxGlobal, PETSC_OWN_POINTER, &isGlobal);
500: {
501: Vec local, global;
502: VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm), 1, stag->entries, PETSC_DECIDE, NULL, &global);
503: VecCreateSeqWithArray(PETSC_COMM_SELF, stag->entriesPerElement, stag->entriesGhost, NULL, &local);
504: VecScatterCreate(local, isLocal, global, isGlobal, &stag->ltog_injective);
505: VecDestroy(&global);
506: VecDestroy(&local);
507: }
508: ISDestroy(&isLocal);
509: ISDestroy(&isGlobal);
510: return 0;
511: }
513: PETSC_INTERN PetscErrorCode DMCreateMatrix_Stag_1D_AIJ_Assemble(DM dm, Mat A)
514: {
515: DMStagStencilType stencil_type;
516: PetscInt dof[2], start, n, n_extra, stencil_width, N, epe;
517: DMBoundaryType boundary_type_x;
519: DMStagGetDOF(dm, &dof[0], &dof[1], NULL, NULL);
520: DMStagGetStencilType(dm, &stencil_type);
521: DMStagGetStencilWidth(dm, &stencil_width);
522: DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, &n_extra, NULL, NULL);
523: DMStagGetGlobalSizes(dm, &N, NULL, NULL);
524: DMStagGetEntriesPerElement(dm, &epe);
525: DMStagGetBoundaryTypes(dm, &boundary_type_x, NULL, NULL);
526: if (stencil_type == DMSTAG_STENCIL_NONE) {
527: /* Couple all DOF at each location to each other */
528: DMStagStencil *row_vertex, *row_element;
530: PetscMalloc1(dof[0], &row_vertex);
531: for (PetscInt c = 0; c < dof[0]; ++c) {
532: row_vertex[c].loc = DMSTAG_LEFT;
533: row_vertex[c].c = c;
534: }
536: PetscMalloc1(dof[1], &row_element);
537: for (PetscInt c = 0; c < dof[1]; ++c) {
538: row_element[c].loc = DMSTAG_ELEMENT;
539: row_element[c].c = c;
540: }
542: for (PetscInt e = start; e < start + n + n_extra; ++e) {
543: {
544: for (PetscInt c = 0; c < dof[0]; ++c) row_vertex[c].i = e;
545: DMStagMatSetValuesStencil(dm, A, dof[0], row_vertex, dof[0], row_vertex, NULL, INSERT_VALUES);
546: }
547: if (e < N) {
548: for (PetscInt c = 0; c < dof[1]; ++c) row_element[c].i = e;
549: DMStagMatSetValuesStencil(dm, A, dof[1], row_element, dof[1], row_element, NULL, INSERT_VALUES);
550: }
551: }
552: PetscFree(row_vertex);
553: PetscFree(row_element);
554: } else if (stencil_type == DMSTAG_STENCIL_STAR || stencil_type == DMSTAG_STENCIL_BOX) {
555: DMStagStencil *col, *row;
557: PetscMalloc1(epe, &row);
558: {
559: PetscInt nrows = 0;
560: for (PetscInt c = 0; c < dof[0]; ++c) {
561: row[nrows].c = c;
562: row[nrows].loc = DMSTAG_LEFT;
563: ++nrows;
564: }
565: for (PetscInt c = 0; c < dof[1]; ++c) {
566: row[nrows].c = c;
567: row[nrows].loc = DMSTAG_ELEMENT;
568: ++nrows;
569: }
570: }
571: PetscMalloc1(epe, &col);
572: {
573: PetscInt ncols = 0;
574: for (PetscInt c = 0; c < dof[0]; ++c) {
575: col[ncols].c = c;
576: col[ncols].loc = DMSTAG_LEFT;
577: ++ncols;
578: }
579: for (PetscInt c = 0; c < dof[1]; ++c) {
580: col[ncols].c = c;
581: col[ncols].loc = DMSTAG_ELEMENT;
582: ++ncols;
583: }
584: }
585: for (PetscInt e = start; e < start + n + n_extra; ++e) {
586: for (PetscInt i = 0; i < epe; ++i) row[i].i = e;
587: for (PetscInt offset = -stencil_width; offset <= stencil_width; ++offset) {
588: const PetscInt e_offset = e + offset;
590: /* Only set values corresponding to elements which can have non-dummy entries,
591: meaning those that map to unknowns in the global representation. In the periodic
592: case, this is the entire stencil, but in all other cases, only includes a single
593: "extra" element which is partially outside the physical domain (those points in the
594: global representation */
595: if (boundary_type_x == DM_BOUNDARY_PERIODIC || (e_offset < N + 1 && e_offset >= 0)) {
596: for (PetscInt i = 0; i < epe; ++i) col[i].i = e_offset;
597: DMStagMatSetValuesStencil(dm, A, epe, row, epe, col, NULL, INSERT_VALUES);
598: }
599: }
600: }
601: PetscFree(row);
602: PetscFree(col);
603: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported stencil type %s", DMStagStencilTypes[stencil_type]);
604: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
605: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
606: return 0;
607: }