Actual source code: stag1d.c

petsc-master 2019-10-13
Report Typos and Errors
  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 faces, edges, 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 grid points
 15: . dof0 - number of degrees of freedom per vertex/point/node/0-cell
 16: . dof1 - number of degrees of freedom per element/edge/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() a 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:   Notes:
 31:   You must call DMSetUp() after this call before using the DM.
 32:   If you wish to use the options database (see the keys above) to change values in the DMStag, you must call
 33:   DMSetFromOptions() after this function but before DMSetUp().

 35:   Level: beginner

 37: .seealso: 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: {
 42:   DM_Stag        *stag;
 43:   PetscMPIInt    size;

 46:   MPI_Comm_size(comm,&size);
 47:   DMCreate(comm,dm);
 48:   DMSetDimension(*dm,1); /* Must precede DMSetType */
 49:   DMSetType(*dm,DMSTAG);
 50:   stag = (DM_Stag*)(*dm)->data;

 52:   /* Global sizes and flags (derived quantities set in DMSetUp_Stag) */
 53:   stag->boundaryType[0] = bndx;
 54:   stag->N[0]            = M;
 55:   stag->nRanks[0]       = size;
 56:   stag->stencilType     = stencilType;
 57:   stag->stencilWidth    = stencilWidth;
 58:   DMStagSetDOF(*dm,dof0,dof1,0,0);
 59:   DMStagSetOwnershipRanges(*dm,lx,NULL,NULL);

 61:   return(0);
 62: }

 64: PETSC_INTERN PetscErrorCode DMStagSetUniformCoordinatesExplicit_1d(DM dm,PetscReal xmin,PetscReal xmax)
 65: {
 67:   DM_Stag        *stagCoord;
 68:   DM             dmCoord;
 69:   Vec            coordLocal,coord;
 70:   PetscReal      h,min;
 71:   PetscScalar    **arr;
 72:   PetscInt       ind,start,n,nExtra,s;
 73:   PetscInt       ileft,ielement;

 76:   DMGetCoordinateDM(dm, &dmCoord);
 77:   stagCoord = (DM_Stag*) dmCoord->data;
 78:   for (s=0; s<2; ++s) {
 79:     if (stagCoord->dof[s] !=0 && stagCoord->dof[s] != 1) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Coordinate DM in 1 dimensions must have 0 or 1 dof on each stratum, but stratum %d has %d dof",s,stagCoord->dof[s]);
 80:   }
 81:   DMGetLocalVector(dmCoord,&coordLocal);

 83:   DMStagVecGetArrayDOF(dmCoord,coordLocal,&arr);
 84:   if (stagCoord->dof[0]) {
 85:     DMStagGetLocationSlot(dmCoord,DMSTAG_LEFT,0,&ileft);
 86:   }
 87:   if (stagCoord->dof[1]) {
 88:     DMStagGetLocationSlot(dmCoord,DMSTAG_ELEMENT,0,&ielement);
 89:   }
 90:   DMStagGetCorners(dmCoord,&start,NULL,NULL,&n,NULL,NULL,&nExtra,NULL,NULL);

 92:   min = xmin;
 93:   h = (xmax-xmin)/stagCoord->N[0];

 95:   for(ind=start; ind<start + n + nExtra; ++ind) {
 96:     if (stagCoord->dof[0]) {
 97:       const PetscReal off = 0.0;
 98:         arr[ind][ileft] = min + ((PetscReal)ind + off) * h;
 99:     }
100:     if (stagCoord->dof[1]) {
101:       const PetscReal off = 0.5;
102:         arr[ind][ielement] = min + ((PetscReal)ind + off) * h;
103:     }
104:   }
105:   DMStagVecRestoreArrayDOF(dmCoord,coordLocal,&arr);
106:   DMCreateGlobalVector(dmCoord,&coord);
107:   DMLocalToGlobalBegin(dmCoord,coordLocal,INSERT_VALUES,coord);
108:   DMLocalToGlobalEnd(dmCoord,coordLocal,INSERT_VALUES,coord);
109:   DMSetCoordinates(dm,coord);
110:   PetscLogObjectParent((PetscObject)dm,(PetscObject)coord);
111:   VecDestroy(&coord);
112:   DMRestoreLocalVector(dmCoord,&coordLocal);
113:   return(0);
114: }

116: /* Helper functions used in DMSetUp_Stag() */
117: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM);

119: PETSC_INTERN PetscErrorCode DMSetUp_Stag_1d(DM dm)
120: {
121:   PetscErrorCode  ierr;
122:   DM_Stag * const stag = (DM_Stag*)dm->data;
123:   PetscMPIInt     size,rank;
124:   MPI_Comm        comm;
125:   PetscInt        j;

128:   PetscObjectGetComm((PetscObject)dm,&comm);
129:   MPI_Comm_size(comm,&size);
130:   MPI_Comm_rank(comm,&rank);

132:   /* Check Global size */
133:   if (stag->N[0] < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Global grid size of %D < 1 specified",stag->N[0]);

135:   /* Local sizes */
136:   if (stag->N[0] < size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"More ranks (%d) than elements (%D) specified",size,stag->N[0]);
137:   if (!stag->l[0]) {
138:     /* Divide equally, giving an extra elements to higher ranks */
139:     PetscMalloc1(stag->nRanks[0],&stag->l[0]);
140:     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);
141:   }
142:   {
143:     PetscInt Nchk = 0;
144:     for (j=0; j<size; ++j) Nchk += stag->l[0][j];
145:     if (Nchk != stag->N[0]) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Sum of specified local sizes (%D) is not equal to global size (%D)",Nchk,stag->N[0]);
146:   }
147:   stag->n[0] = stag->l[0][rank];

149:   /* Rank (trivial in 1d) */
150:   stag->rank[0]      = rank;
151:   stag->firstRank[0] = (PetscBool)(rank == 0);
152:   stag->lastRank[0]  = (PetscBool)(rank == size-1);

154:   /* Local (unghosted) numbers of entries */
155:   stag->entriesPerElement = stag->dof[0] + stag->dof[1];
156:   switch (stag->boundaryType[0]) {
157:     case DM_BOUNDARY_NONE:
158:     case DM_BOUNDARY_GHOSTED:  stag->entries = stag->n[0] * stag->entriesPerElement + (stag->lastRank[0] ?  stag->dof[0] : 0); break;
159:     case DM_BOUNDARY_PERIODIC: stag->entries = stag->n[0] * stag->entriesPerElement;                                           break;
160:     default: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
161:   }

163:   /* Starting element */
164:   stag->start[0] = 0;
165:   for(j=0; j<stag->rank[0]; ++j) stag->start[0] += stag->l[0][j];

167:   /* Local/ghosted size and starting element */
168:   switch (stag->boundaryType[0]) {
169:     case DM_BOUNDARY_NONE :
170:       switch (stag->stencilType) {
171:         case DMSTAG_STENCIL_NONE : /* Only dummy cells on the right */
172:           stag->startGhost[0] = stag->start[0];
173:           stag->nGhost[0]     = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
174:           break;
175:         case DMSTAG_STENCIL_STAR :
176:         case DMSTAG_STENCIL_BOX :
177:           stag->startGhost[0] = stag->firstRank[0] ? stag->start[0]: stag->start[0] - stag->stencilWidth;
178:           stag->nGhost[0] = stag->n[0];
179:           stag->nGhost[0] += stag->firstRank[0] ? 0 : stag->stencilWidth;
180:           stag->nGhost[0] += stag->lastRank[0]  ? 1 : stag->stencilWidth;
181:           break;
182:         default :
183:           SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
184:       }
185:       break;
186:     case DM_BOUNDARY_GHOSTED:
187:       switch (stag->stencilType) {
188:         case DMSTAG_STENCIL_NONE :
189:           stag->startGhost[0] = stag->start[0];
190:           stag->nGhost[0]     = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
191:           break;
192:         case DMSTAG_STENCIL_STAR :
193:         case DMSTAG_STENCIL_BOX :
194:           stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
195:           stag->nGhost[0]     = stag->n[0] + 2*stag->stencilWidth + (stag->lastRank[0] && stag->stencilWidth == 0 ? 1 : 0);
196:           break;
197:         default :
198:           SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
199:       }
200:       break;
201:     case DM_BOUNDARY_PERIODIC:
202:       switch (stag->stencilType) {
203:         case DMSTAG_STENCIL_NONE :
204:           stag->startGhost[0] = stag->start[0];
205:           stag->nGhost[0]     = stag->n[0];
206:           break;
207:         case DMSTAG_STENCIL_STAR :
208:         case DMSTAG_STENCIL_BOX :
209:           stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
210:           stag->nGhost[0]     = stag->n[0] + 2*stag->stencilWidth;
211:           break;
212:         default :
213:           SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
214:       }
215:       break;
216:     default :
217:       SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
218:   }

220:   /* Total size of ghosted/local represention */
221:   stag->entriesGhost = stag->nGhost[0]*stag->entriesPerElement;

223:   /* Define neighbors */
224:   PetscMalloc1(3,&stag->neighbors);
225:   if (stag->firstRank[0]) {
226:     switch (stag->boundaryType[0]) {
227:       case DM_BOUNDARY_GHOSTED:
228:       case DM_BOUNDARY_NONE:     stag->neighbors[0] = -1;                break;
229:       case DM_BOUNDARY_PERIODIC: stag->neighbors[0] = stag->nRanks[0]-1; break;
230:       default : SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
231:     }
232:   } else {
233:     stag->neighbors[0] = stag->rank[0]-1;
234:   }
235:   stag->neighbors[1] = stag->rank[0];
236:   if (stag->lastRank[0]) {
237:     switch (stag->boundaryType[0]) {
238:       case DM_BOUNDARY_GHOSTED:
239:       case DM_BOUNDARY_NONE:     stag->neighbors[2] = -1;                break;
240:       case DM_BOUNDARY_PERIODIC: stag->neighbors[2] = 0;                 break;
241:       default : SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
242:     }
243:   } else {
244:     stag->neighbors[2] = stag->rank[0]+1;
245:   }

247:   if (stag->n[0] < stag->stencilWidth) {
248:     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMStag 1d setup does not support local sizes (%d) smaller than the elementwise stencil width (%d)",stag->n[0],stag->stencilWidth);
249:   }

251:   /* Create global->local VecScatter and ISLocalToGlobalMapping */
252:   {
253:     PetscInt *idxLocal,*idxGlobal,*idxGlobalAll;
254:     PetscInt i,iLocal,d,entriesToTransferTotal,ghostOffsetStart,ghostOffsetEnd,nNonDummyGhost;
255:     IS       isLocal,isGlobal;

257:     /* The offset on the right (may not be equal to the stencil width, as we
258:        always have at least one ghost element, to account for the boundary
259:        point, and may with ghosted boundaries), and the number of non-dummy ghost elements */
260:     ghostOffsetStart = stag->start[0] - stag->startGhost[0];
261:     ghostOffsetEnd   = stag->startGhost[0]+stag->nGhost[0] - (stag->start[0]+stag->n[0]);
262:     nNonDummyGhost   = stag->nGhost[0] - (stag->lastRank[0] ? ghostOffsetEnd : 0) - (stag->firstRank[0] ? ghostOffsetStart : 0);

264:     /* Compute the number of non-dummy entries in the local representation
265:        This is equal to the number of non-dummy elements in the local (ghosted) representation,
266:        plus some extra entries on the right boundary on the last rank*/
267:     switch (stag->boundaryType[0]) {
268:       case DM_BOUNDARY_GHOSTED:
269:       case DM_BOUNDARY_NONE:
270:         entriesToTransferTotal = nNonDummyGhost * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0);
271:         break;
272:       case DM_BOUNDARY_PERIODIC:
273:         entriesToTransferTotal = stag->entriesGhost; /* No dummy points */
274:         break;
275:       default :
276:         SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
277:     }

279:     PetscMalloc1(entriesToTransferTotal,&idxLocal);
280:     PetscMalloc1(entriesToTransferTotal,&idxGlobal);
281:     PetscMalloc1(stag->entriesGhost,&idxGlobalAll);
282:     if (stag->boundaryType[0] == DM_BOUNDARY_NONE) {
283:       PetscInt count = 0,countAll = 0;
284:       /* Left ghost points and native points */
285:       for (i=stag->startGhost[0], iLocal=0; iLocal<nNonDummyGhost; ++i,++iLocal) {
286:         for (d=0; d<stag->entriesPerElement; ++d,++count,++countAll) {
287:           idxLocal [count]       = iLocal * stag->entriesPerElement + d;
288:           idxGlobal[count]       = i      * stag->entriesPerElement + d;
289:           idxGlobalAll[countAll] = i      * stag->entriesPerElement + d;
290:         }
291:       }
292:       /* Ghost points on the right
293:          Special case for last (partial dummy) element on the last rank */
294:       if (stag->lastRank[0] ) {
295:         i      = stag->N[0];
296:         iLocal = (stag->nGhost[0]-ghostOffsetEnd);
297:         /* Only vertex (0-cell) dofs in global representation */
298:         for (d=0; d<stag->dof[0]; ++d,++count,++countAll) {
299:           idxGlobal[count]       = i      * stag->entriesPerElement + d;
300:           idxLocal [count]       = iLocal * stag->entriesPerElement + d;
301:           idxGlobalAll[countAll] = i      * stag->entriesPerElement + d;
302:         }
303:         for (d=stag->dof[0]; d<stag->entriesPerElement; ++d,++countAll) { /* Additional dummy entries */
304:           idxGlobalAll[countAll] = -1;
305:         }
306:       }
307:     } else if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC) {
308:       PetscInt count = 0,iLocal = 0; /* No dummy points, so idxGlobal and idxGlobalAll are identical */
309:       const PetscInt iMin = stag->firstRank[0] ? stag->start[0] : stag->startGhost[0];
310:       const PetscInt iMax = stag->lastRank[0] ? stag->startGhost[0] + stag->nGhost[0] - stag->stencilWidth : stag->startGhost[0] + stag->nGhost[0];
311:       /* Ghost points on the left */
312:       if (stag->firstRank[0]) {
313:         for (i=stag->N[0]-stag->stencilWidth; iLocal<stag->stencilWidth; ++i,++iLocal) {
314:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
315:             idxGlobal[count] = i      * stag->entriesPerElement + d;
316:             idxLocal [count] = iLocal * stag->entriesPerElement + d;
317:             idxGlobalAll[count] = idxGlobal[count];
318:           }
319:         }
320:       }
321:       /* Native points */
322:       for (i=iMin; i<iMax; ++i,++iLocal) {
323:         for (d=0; d<stag->entriesPerElement; ++d,++count) {
324:           idxGlobal[count] = i      * stag->entriesPerElement + d;
325:           idxLocal [count] = iLocal * stag->entriesPerElement + d;
326:           idxGlobalAll[count] = idxGlobal[count];
327:         }
328:       }
329:       /* Ghost points on the right */
330:       if (stag->lastRank[0]) {
331:         for (i=0; iLocal<stag->nGhost[0]; ++i,++iLocal) {
332:           for (d=0; d<stag->entriesPerElement; ++d,++count) {
333:             idxGlobal[count] = i      * stag->entriesPerElement + d;
334:             idxLocal [count] = iLocal * stag->entriesPerElement + d;
335:             idxGlobalAll[count] = idxGlobal[count];
336:           }
337:         }
338:       }
339:     } else if (stag->boundaryType[0] == DM_BOUNDARY_GHOSTED) {
340:       PetscInt count = 0,countAll = 0;
341:       /* Dummy elements on the left, on the first rank */
342:       if (stag->firstRank[0]) {
343:         for(iLocal=0; iLocal<ghostOffsetStart; ++iLocal) {
344:           /* Complete elements full of dummy entries */
345:           for (d=0; d<stag->entriesPerElement; ++d,++countAll) {
346:             idxGlobalAll[countAll] = -1;
347:           }
348:         }
349:         i = 0; /* nonDummy entries start with global entry 0 */
350:       } else {
351:         /* nonDummy entries start as usual */
352:         i = stag->startGhost[0];
353:         iLocal = 0;
354:       }

356:       /* non-Dummy entries */
357:       {
358:         PetscInt iLocalNonDummyMax = stag->firstRank[0] ? nNonDummyGhost + ghostOffsetStart : nNonDummyGhost;
359:         for (; iLocal<iLocalNonDummyMax; ++i,++iLocal) {
360:           for (d=0; d<stag->entriesPerElement; ++d,++count,++countAll) {
361:             idxLocal [count]       = iLocal * stag->entriesPerElement + d;
362:             idxGlobal[count]       = i      * stag->entriesPerElement + d;
363:             idxGlobalAll[countAll] = i      * stag->entriesPerElement + d;
364:           }
365:         }
366:       }

368:       /* (partial) dummy elements on the right, on the last rank */
369:       if (stag->lastRank[0]) {
370:         /* First one is partial dummy */
371:         i      = stag->N[0];
372:         iLocal = (stag->nGhost[0]-ghostOffsetEnd);
373:         for (d=0; d<stag->dof[0]; ++d,++count,++countAll) { /* Only vertex (0-cell) dofs in global representation */
374:           idxLocal [count]       = iLocal * stag->entriesPerElement + d;
375:           idxGlobal[count]       = i      * stag->entriesPerElement + d;
376:           idxGlobalAll[countAll] = i      * stag->entriesPerElement + d;
377:         }
378:         for (d=stag->dof[0]; d<stag->entriesPerElement; ++d,++countAll) { /* Additional dummy entries */
379:           idxGlobalAll[countAll] = -1;
380:         }
381:         for (iLocal = stag->nGhost[0] - ghostOffsetEnd + 1; iLocal < stag->nGhost[0]; ++iLocal) {
382:           /* Additional dummy elements */
383:           for (d=0; d<stag->entriesPerElement; ++d,++countAll) {
384:             idxGlobalAll[countAll] = -1;
385:           }
386:         }
387:       }
388:     } else SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);

390:     /* Create Local IS (transferring pointer ownership) */
391:     ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxLocal,PETSC_OWN_POINTER,&isLocal);

393:     /* Create Global IS (transferring pointer ownership) */
394:     ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxGlobal,PETSC_OWN_POINTER,&isGlobal);

396:     /* Create stag->gtol, which doesn't include dummy entries */
397:     {
398:       Vec local,global;
399:       VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
400:       VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
401:       VecScatterCreate(global,isGlobal,local,isLocal,&stag->gtol);
402:       VecDestroy(&global);
403:       VecDestroy(&local);
404:     }

406:     /* In special cases, create a dedicated injective local-to-global map */
407:     if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) {
408:       DMStagPopulateLocalToGlobalInjective(dm);
409:     }

411:     /* Destroy ISs */
412:     ISDestroy(&isLocal);
413:     ISDestroy(&isGlobal);

415:     /* Create local-to-global map (transferring pointer ownership) */
416:     ISLocalToGlobalMappingCreate(comm,1,stag->entriesGhost,idxGlobalAll,PETSC_OWN_POINTER,&dm->ltogmap);
417:     PetscLogObjectParent((PetscObject)dm,(PetscObject)dm->ltogmap);
418:   }

420:   /* Precompute location offsets */
421:   DMStagComputeLocationOffsets_1d(dm);

423:   /* View from Options */
424:   DMViewFromOptions(dm,NULL,"-dm_view");

426:  return(0);
427: }

429: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM dm)
430: {
431:   PetscErrorCode  ierr;
432:   DM_Stag * const stag = (DM_Stag*)dm->data;
433:   const PetscInt  epe = stag->entriesPerElement;

436:   PetscMalloc1(DMSTAG_NUMBER_LOCATIONS,&stag->locationOffsets);
437:   stag->locationOffsets[DMSTAG_LEFT]    = 0;
438:   stag->locationOffsets[DMSTAG_ELEMENT] = stag->locationOffsets[DMSTAG_LEFT] + stag->dof[0];
439:   stag->locationOffsets[DMSTAG_RIGHT]   = stag->locationOffsets[DMSTAG_LEFT] + epe;
440:   return(0);
441: }

443: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_1d(DM dm)
444: {
445:   PetscErrorCode  ierr;
446:   DM_Stag * const stag = (DM_Stag*)dm->data;
447:   PetscInt        *idxLocal,*idxGlobal;
448:   PetscInt        i,iLocal,d,count;
449:   IS              isLocal,isGlobal;

452:   PetscMalloc1(stag->entries,&idxLocal);
453:   PetscMalloc1(stag->entries,&idxGlobal);
454:   count = 0;
455:   iLocal = stag->start[0]-stag->startGhost[0];
456:   for (i=stag->start[0]; i<stag->start[0]+stag->n[0]; ++i,++iLocal) {
457:     for (d=0; d<stag->entriesPerElement; ++d,++count) {
458:       idxGlobal[count] = i      * stag->entriesPerElement + d;
459:       idxLocal [count] = iLocal * stag->entriesPerElement + d;
460:     }
461:   }
462:   if (stag->lastRank[0] && stag->boundaryType[0] != DM_BOUNDARY_PERIODIC) {
463:     i = stag->start[0]+stag->n[0];
464:     iLocal = stag->start[0]-stag->startGhost[0] + stag->n[0];
465:     for (d=0; d<stag->dof[0]; ++d,++count) {
466:       idxGlobal[count] = i      * stag->entriesPerElement + d;
467:       idxLocal [count] = iLocal * stag->entriesPerElement + d;
468:     }
469:   }
470:   ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxLocal,PETSC_OWN_POINTER,&isLocal);
471:   ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
472:   {
473:     Vec local,global;
474:     VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
475:     VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
476:     VecScatterCreate(local,isLocal,global,isGlobal,&stag->ltog_injective);
477:     VecDestroy(&global);
478:     VecDestroy(&local);
479:   }
480:   ISDestroy(&isLocal);
481:   ISDestroy(&isGlobal);
482:   return(0);
483: }