Actual source code: stag1d.c
petsc-3.12.5 2020-03-29
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: }