Actual source code: stag1d.c
petsc-3.13.6 2020-09-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: PetscMPIInt size;
45: MPI_Comm_size(comm,&size);
46: DMCreate(comm,dm);
47: DMSetDimension(*dm,1);
48: DMStagInitialize(bndx,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,M,0,0,size,0,0,dof0,dof1,0,0,stencilType,stencilWidth,lx,NULL,NULL,*dm);
49: return(0);
50: }
52: PETSC_INTERN PetscErrorCode DMStagSetUniformCoordinatesExplicit_1d(DM dm,PetscReal xmin,PetscReal xmax)
53: {
55: DM_Stag *stagCoord;
56: DM dmCoord;
57: Vec coordLocal,coord;
58: PetscReal h,min;
59: PetscScalar **arr;
60: PetscInt ind,start,n,nExtra,s;
61: PetscInt ileft,ielement;
64: DMGetCoordinateDM(dm, &dmCoord);
65: stagCoord = (DM_Stag*) dmCoord->data;
66: for (s=0; s<2; ++s) {
67: 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]);
68: }
69: DMGetLocalVector(dmCoord,&coordLocal);
71: DMStagVecGetArray(dmCoord,coordLocal,&arr);
72: if (stagCoord->dof[0]) {
73: DMStagGetLocationSlot(dmCoord,DMSTAG_LEFT,0,&ileft);
74: }
75: if (stagCoord->dof[1]) {
76: DMStagGetLocationSlot(dmCoord,DMSTAG_ELEMENT,0,&ielement);
77: }
78: DMStagGetCorners(dmCoord,&start,NULL,NULL,&n,NULL,NULL,&nExtra,NULL,NULL);
80: min = xmin;
81: h = (xmax-xmin)/stagCoord->N[0];
83: for(ind=start; ind<start + n + nExtra; ++ind) {
84: if (stagCoord->dof[0]) {
85: const PetscReal off = 0.0;
86: arr[ind][ileft] = min + ((PetscReal)ind + off) * h;
87: }
88: if (stagCoord->dof[1]) {
89: const PetscReal off = 0.5;
90: arr[ind][ielement] = min + ((PetscReal)ind + off) * h;
91: }
92: }
93: DMStagVecRestoreArray(dmCoord,coordLocal,&arr);
94: DMCreateGlobalVector(dmCoord,&coord);
95: DMLocalToGlobalBegin(dmCoord,coordLocal,INSERT_VALUES,coord);
96: DMLocalToGlobalEnd(dmCoord,coordLocal,INSERT_VALUES,coord);
97: DMSetCoordinates(dm,coord);
98: PetscLogObjectParent((PetscObject)dm,(PetscObject)coord);
99: VecDestroy(&coord);
100: DMRestoreLocalVector(dmCoord,&coordLocal);
101: return(0);
102: }
104: /* Helper functions used in DMSetUp_Stag() */
105: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM);
107: PETSC_INTERN PetscErrorCode DMSetUp_Stag_1d(DM dm)
108: {
109: PetscErrorCode ierr;
110: DM_Stag * const stag = (DM_Stag*)dm->data;
111: PetscMPIInt size,rank;
112: MPI_Comm comm;
113: PetscInt j;
116: PetscObjectGetComm((PetscObject)dm,&comm);
117: MPI_Comm_size(comm,&size);
118: MPI_Comm_rank(comm,&rank);
120: /* Check Global size */
121: if (stag->N[0] < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Global grid size of %D < 1 specified",stag->N[0]);
123: /* Local sizes */
124: if (stag->N[0] < size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"More ranks (%d) than elements (%D) specified",size,stag->N[0]);
125: if (!stag->l[0]) {
126: /* Divide equally, giving an extra elements to higher ranks */
127: PetscMalloc1(stag->nRanks[0],&stag->l[0]);
128: 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);
129: }
130: {
131: PetscInt Nchk = 0;
132: for (j=0; j<size; ++j) Nchk += stag->l[0][j];
133: 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]);
134: }
135: stag->n[0] = stag->l[0][rank];
137: /* Rank (trivial in 1d) */
138: stag->rank[0] = rank;
139: stag->firstRank[0] = (PetscBool)(rank == 0);
140: stag->lastRank[0] = (PetscBool)(rank == size-1);
142: /* Local (unghosted) numbers of entries */
143: stag->entriesPerElement = stag->dof[0] + stag->dof[1];
144: switch (stag->boundaryType[0]) {
145: case DM_BOUNDARY_NONE:
146: case DM_BOUNDARY_GHOSTED: stag->entries = stag->n[0] * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0); break;
147: case DM_BOUNDARY_PERIODIC: stag->entries = stag->n[0] * stag->entriesPerElement; break;
148: default: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
149: }
151: /* Starting element */
152: stag->start[0] = 0;
153: for(j=0; j<stag->rank[0]; ++j) stag->start[0] += stag->l[0][j];
155: /* Local/ghosted size and starting element */
156: switch (stag->boundaryType[0]) {
157: case DM_BOUNDARY_NONE :
158: switch (stag->stencilType) {
159: case DMSTAG_STENCIL_NONE : /* Only dummy cells on the right */
160: stag->startGhost[0] = stag->start[0];
161: stag->nGhost[0] = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
162: break;
163: case DMSTAG_STENCIL_STAR :
164: case DMSTAG_STENCIL_BOX :
165: stag->startGhost[0] = stag->firstRank[0] ? stag->start[0]: stag->start[0] - stag->stencilWidth;
166: stag->nGhost[0] = stag->n[0];
167: stag->nGhost[0] += stag->firstRank[0] ? 0 : stag->stencilWidth;
168: stag->nGhost[0] += stag->lastRank[0] ? 1 : stag->stencilWidth;
169: break;
170: default :
171: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
172: }
173: break;
174: case DM_BOUNDARY_GHOSTED:
175: switch (stag->stencilType) {
176: case DMSTAG_STENCIL_NONE :
177: stag->startGhost[0] = stag->start[0];
178: stag->nGhost[0] = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
179: break;
180: case DMSTAG_STENCIL_STAR :
181: case DMSTAG_STENCIL_BOX :
182: stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
183: stag->nGhost[0] = stag->n[0] + 2*stag->stencilWidth + (stag->lastRank[0] && stag->stencilWidth == 0 ? 1 : 0);
184: break;
185: default :
186: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
187: }
188: break;
189: case DM_BOUNDARY_PERIODIC:
190: switch (stag->stencilType) {
191: case DMSTAG_STENCIL_NONE :
192: stag->startGhost[0] = stag->start[0];
193: stag->nGhost[0] = stag->n[0];
194: break;
195: case DMSTAG_STENCIL_STAR :
196: case DMSTAG_STENCIL_BOX :
197: stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
198: stag->nGhost[0] = stag->n[0] + 2*stag->stencilWidth;
199: break;
200: default :
201: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
202: }
203: break;
204: default :
205: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
206: }
208: /* Total size of ghosted/local represention */
209: stag->entriesGhost = stag->nGhost[0]*stag->entriesPerElement;
211: /* Define neighbors */
212: PetscMalloc1(3,&stag->neighbors);
213: if (stag->firstRank[0]) {
214: switch (stag->boundaryType[0]) {
215: case DM_BOUNDARY_GHOSTED:
216: case DM_BOUNDARY_NONE: stag->neighbors[0] = -1; break;
217: case DM_BOUNDARY_PERIODIC: stag->neighbors[0] = stag->nRanks[0]-1; break;
218: default : SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
219: }
220: } else {
221: stag->neighbors[0] = stag->rank[0]-1;
222: }
223: stag->neighbors[1] = stag->rank[0];
224: if (stag->lastRank[0]) {
225: switch (stag->boundaryType[0]) {
226: case DM_BOUNDARY_GHOSTED:
227: case DM_BOUNDARY_NONE: stag->neighbors[2] = -1; break;
228: case DM_BOUNDARY_PERIODIC: stag->neighbors[2] = 0; break;
229: default : SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
230: }
231: } else {
232: stag->neighbors[2] = stag->rank[0]+1;
233: }
235: if (stag->n[0] < stag->stencilWidth) {
236: 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);
237: }
239: /* Create global->local VecScatter and ISLocalToGlobalMapping */
240: {
241: PetscInt *idxLocal,*idxGlobal,*idxGlobalAll;
242: PetscInt i,iLocal,d,entriesToTransferTotal,ghostOffsetStart,ghostOffsetEnd,nNonDummyGhost;
243: IS isLocal,isGlobal;
245: /* The offset on the right (may not be equal to the stencil width, as we
246: always have at least one ghost element, to account for the boundary
247: point, and may with ghosted boundaries), and the number of non-dummy ghost elements */
248: ghostOffsetStart = stag->start[0] - stag->startGhost[0];
249: ghostOffsetEnd = stag->startGhost[0]+stag->nGhost[0] - (stag->start[0]+stag->n[0]);
250: nNonDummyGhost = stag->nGhost[0] - (stag->lastRank[0] ? ghostOffsetEnd : 0) - (stag->firstRank[0] ? ghostOffsetStart : 0);
252: /* Compute the number of non-dummy entries in the local representation
253: This is equal to the number of non-dummy elements in the local (ghosted) representation,
254: plus some extra entries on the right boundary on the last rank*/
255: switch (stag->boundaryType[0]) {
256: case DM_BOUNDARY_GHOSTED:
257: case DM_BOUNDARY_NONE:
258: entriesToTransferTotal = nNonDummyGhost * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0);
259: break;
260: case DM_BOUNDARY_PERIODIC:
261: entriesToTransferTotal = stag->entriesGhost; /* No dummy points */
262: break;
263: default :
264: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
265: }
267: PetscMalloc1(entriesToTransferTotal,&idxLocal);
268: PetscMalloc1(entriesToTransferTotal,&idxGlobal);
269: PetscMalloc1(stag->entriesGhost,&idxGlobalAll);
270: if (stag->boundaryType[0] == DM_BOUNDARY_NONE) {
271: PetscInt count = 0,countAll = 0;
272: /* Left ghost points and native points */
273: for (i=stag->startGhost[0], iLocal=0; iLocal<nNonDummyGhost; ++i,++iLocal) {
274: for (d=0; d<stag->entriesPerElement; ++d,++count,++countAll) {
275: idxLocal [count] = iLocal * stag->entriesPerElement + d;
276: idxGlobal[count] = i * stag->entriesPerElement + d;
277: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
278: }
279: }
280: /* Ghost points on the right
281: Special case for last (partial dummy) element on the last rank */
282: if (stag->lastRank[0] ) {
283: i = stag->N[0];
284: iLocal = (stag->nGhost[0]-ghostOffsetEnd);
285: /* Only vertex (0-cell) dofs in global representation */
286: for (d=0; d<stag->dof[0]; ++d,++count,++countAll) {
287: idxGlobal[count] = i * stag->entriesPerElement + d;
288: idxLocal [count] = iLocal * stag->entriesPerElement + d;
289: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
290: }
291: for (d=stag->dof[0]; d<stag->entriesPerElement; ++d,++countAll) { /* Additional dummy entries */
292: idxGlobalAll[countAll] = -1;
293: }
294: }
295: } else if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC) {
296: PetscInt count = 0,iLocal = 0; /* No dummy points, so idxGlobal and idxGlobalAll are identical */
297: const PetscInt iMin = stag->firstRank[0] ? stag->start[0] : stag->startGhost[0];
298: const PetscInt iMax = stag->lastRank[0] ? stag->startGhost[0] + stag->nGhost[0] - stag->stencilWidth : stag->startGhost[0] + stag->nGhost[0];
299: /* Ghost points on the left */
300: if (stag->firstRank[0]) {
301: for (i=stag->N[0]-stag->stencilWidth; iLocal<stag->stencilWidth; ++i,++iLocal) {
302: for (d=0; d<stag->entriesPerElement; ++d,++count) {
303: idxGlobal[count] = i * stag->entriesPerElement + d;
304: idxLocal [count] = iLocal * stag->entriesPerElement + d;
305: idxGlobalAll[count] = idxGlobal[count];
306: }
307: }
308: }
309: /* Native points */
310: for (i=iMin; i<iMax; ++i,++iLocal) {
311: for (d=0; d<stag->entriesPerElement; ++d,++count) {
312: idxGlobal[count] = i * stag->entriesPerElement + d;
313: idxLocal [count] = iLocal * stag->entriesPerElement + d;
314: idxGlobalAll[count] = idxGlobal[count];
315: }
316: }
317: /* Ghost points on the right */
318: if (stag->lastRank[0]) {
319: for (i=0; iLocal<stag->nGhost[0]; ++i,++iLocal) {
320: for (d=0; d<stag->entriesPerElement; ++d,++count) {
321: idxGlobal[count] = i * stag->entriesPerElement + d;
322: idxLocal [count] = iLocal * stag->entriesPerElement + d;
323: idxGlobalAll[count] = idxGlobal[count];
324: }
325: }
326: }
327: } else if (stag->boundaryType[0] == DM_BOUNDARY_GHOSTED) {
328: PetscInt count = 0,countAll = 0;
329: /* Dummy elements on the left, on the first rank */
330: if (stag->firstRank[0]) {
331: for(iLocal=0; iLocal<ghostOffsetStart; ++iLocal) {
332: /* Complete elements full of dummy entries */
333: for (d=0; d<stag->entriesPerElement; ++d,++countAll) {
334: idxGlobalAll[countAll] = -1;
335: }
336: }
337: i = 0; /* nonDummy entries start with global entry 0 */
338: } else {
339: /* nonDummy entries start as usual */
340: i = stag->startGhost[0];
341: iLocal = 0;
342: }
344: /* non-Dummy entries */
345: {
346: PetscInt iLocalNonDummyMax = stag->firstRank[0] ? nNonDummyGhost + ghostOffsetStart : nNonDummyGhost;
347: for (; iLocal<iLocalNonDummyMax; ++i,++iLocal) {
348: for (d=0; d<stag->entriesPerElement; ++d,++count,++countAll) {
349: idxLocal [count] = iLocal * stag->entriesPerElement + d;
350: idxGlobal[count] = i * stag->entriesPerElement + d;
351: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
352: }
353: }
354: }
356: /* (partial) dummy elements on the right, on the last rank */
357: if (stag->lastRank[0]) {
358: /* First one is partial dummy */
359: i = stag->N[0];
360: iLocal = (stag->nGhost[0]-ghostOffsetEnd);
361: for (d=0; d<stag->dof[0]; ++d,++count,++countAll) { /* Only vertex (0-cell) dofs in global representation */
362: idxLocal [count] = iLocal * stag->entriesPerElement + d;
363: idxGlobal[count] = i * stag->entriesPerElement + d;
364: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
365: }
366: for (d=stag->dof[0]; d<stag->entriesPerElement; ++d,++countAll) { /* Additional dummy entries */
367: idxGlobalAll[countAll] = -1;
368: }
369: for (iLocal = stag->nGhost[0] - ghostOffsetEnd + 1; iLocal < stag->nGhost[0]; ++iLocal) {
370: /* Additional dummy elements */
371: for (d=0; d<stag->entriesPerElement; ++d,++countAll) {
372: idxGlobalAll[countAll] = -1;
373: }
374: }
375: }
376: } else SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
378: /* Create Local IS (transferring pointer ownership) */
379: ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxLocal,PETSC_OWN_POINTER,&isLocal);
381: /* Create Global IS (transferring pointer ownership) */
382: ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
384: /* Create stag->gtol, which doesn't include dummy entries */
385: {
386: Vec local,global;
387: VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
388: VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
389: VecScatterCreate(global,isGlobal,local,isLocal,&stag->gtol);
390: VecDestroy(&global);
391: VecDestroy(&local);
392: }
394: /* In special cases, create a dedicated injective local-to-global map */
395: if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) {
396: DMStagPopulateLocalToGlobalInjective(dm);
397: }
399: /* Destroy ISs */
400: ISDestroy(&isLocal);
401: ISDestroy(&isGlobal);
403: /* Create local-to-global map (transferring pointer ownership) */
404: ISLocalToGlobalMappingCreate(comm,1,stag->entriesGhost,idxGlobalAll,PETSC_OWN_POINTER,&dm->ltogmap);
405: PetscLogObjectParent((PetscObject)dm,(PetscObject)dm->ltogmap);
406: }
408: /* Precompute location offsets */
409: DMStagComputeLocationOffsets_1d(dm);
411: /* View from Options */
412: DMViewFromOptions(dm,NULL,"-dm_view");
414: return(0);
415: }
417: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM dm)
418: {
419: PetscErrorCode ierr;
420: DM_Stag * const stag = (DM_Stag*)dm->data;
421: const PetscInt epe = stag->entriesPerElement;
424: PetscMalloc1(DMSTAG_NUMBER_LOCATIONS,&stag->locationOffsets);
425: stag->locationOffsets[DMSTAG_LEFT] = 0;
426: stag->locationOffsets[DMSTAG_ELEMENT] = stag->locationOffsets[DMSTAG_LEFT] + stag->dof[0];
427: stag->locationOffsets[DMSTAG_RIGHT] = stag->locationOffsets[DMSTAG_LEFT] + epe;
428: return(0);
429: }
431: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_1d(DM dm)
432: {
433: PetscErrorCode ierr;
434: DM_Stag * const stag = (DM_Stag*)dm->data;
435: PetscInt *idxLocal,*idxGlobal;
436: PetscInt i,iLocal,d,count;
437: IS isLocal,isGlobal;
440: PetscMalloc1(stag->entries,&idxLocal);
441: PetscMalloc1(stag->entries,&idxGlobal);
442: count = 0;
443: iLocal = stag->start[0]-stag->startGhost[0];
444: for (i=stag->start[0]; i<stag->start[0]+stag->n[0]; ++i,++iLocal) {
445: for (d=0; d<stag->entriesPerElement; ++d,++count) {
446: idxGlobal[count] = i * stag->entriesPerElement + d;
447: idxLocal [count] = iLocal * stag->entriesPerElement + d;
448: }
449: }
450: if (stag->lastRank[0] && stag->boundaryType[0] != DM_BOUNDARY_PERIODIC) {
451: i = stag->start[0]+stag->n[0];
452: iLocal = stag->start[0]-stag->startGhost[0] + stag->n[0];
453: for (d=0; d<stag->dof[0]; ++d,++count) {
454: idxGlobal[count] = i * stag->entriesPerElement + d;
455: idxLocal [count] = iLocal * stag->entriesPerElement + d;
456: }
457: }
458: ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxLocal,PETSC_OWN_POINTER,&isLocal);
459: ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
460: {
461: Vec local,global;
462: VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
463: VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
464: VecScatterCreate(local,isLocal,global,isGlobal,&stag->ltog_injective);
465: VecDestroy(&global);
466: VecDestroy(&local);
467: }
468: ISDestroy(&isLocal);
469: ISDestroy(&isGlobal);
470: return(0);
471: }