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