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() 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: {
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 DMStagSetUniformCoordinatesExplicit_1d(DM dm,PetscReal xmin,PetscReal xmax)
51: {
52: DM_Stag *stagCoord;
53: DM dmCoord;
54: Vec coordLocal;
55: PetscReal h,min;
56: PetscScalar **arr;
57: PetscInt start_ghost,n_ghost,s;
58: PetscInt ileft,ielement;
60: DMGetCoordinateDM(dm, &dmCoord);
61: stagCoord = (DM_Stag*) dmCoord->data;
62: for (s=0; s<2; ++s) {
64: }
65: DMCreateLocalVector(dmCoord,&coordLocal);
67: DMStagVecGetArray(dmCoord,coordLocal,&arr);
68: if (stagCoord->dof[0]) {
69: DMStagGetLocationSlot(dmCoord,DMSTAG_LEFT,0,&ileft);
70: }
71: if (stagCoord->dof[1]) {
72: DMStagGetLocationSlot(dmCoord,DMSTAG_ELEMENT,0,&ielement);
73: }
74: DMStagGetGhostCorners(dmCoord,&start_ghost,NULL,NULL,&n_ghost,NULL,NULL);
76: min = xmin;
77: h = (xmax-xmin)/stagCoord->N[0];
79: for (PetscInt ind=start_ghost; ind<start_ghost + n_ghost; ++ind) {
80: if (stagCoord->dof[0]) {
81: const PetscReal off = 0.0;
82: arr[ind][ileft] = min + ((PetscReal)ind + off) * h;
83: }
84: if (stagCoord->dof[1]) {
85: const PetscReal off = 0.5;
86: arr[ind][ielement] = min + ((PetscReal)ind + off) * h;
87: }
88: }
89: DMStagVecRestoreArray(dmCoord,coordLocal,&arr);
90: DMSetCoordinatesLocal(dm,coordLocal);
91: PetscLogObjectParent((PetscObject)dm,(PetscObject)coordLocal);
92: VecDestroy(&coordLocal);
93: return 0;
94: }
96: /* Helper functions used in DMSetUp_Stag() */
97: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM);
99: PETSC_INTERN PetscErrorCode DMSetUp_Stag_1d(DM dm)
100: {
101: DM_Stag * const stag = (DM_Stag*)dm->data;
102: PetscMPIInt size,rank;
103: MPI_Comm comm;
104: PetscInt j;
106: PetscObjectGetComm((PetscObject)dm,&comm);
107: MPI_Comm_size(comm,&size);
108: MPI_Comm_rank(comm,&rank);
110: /* Check Global size */
113: /* Local sizes */
115: if (!stag->l[0]) {
116: /* Divide equally, giving an extra elements to higher ranks */
117: PetscMalloc1(stag->nRanks[0],&stag->l[0]);
118: 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);
119: }
120: {
121: PetscInt Nchk = 0;
122: for (j=0; j<size; ++j) Nchk += stag->l[0][j];
124: }
125: stag->n[0] = stag->l[0][rank];
127: /* Rank (trivial in 1d) */
128: stag->rank[0] = rank;
129: stag->firstRank[0] = (PetscBool)(rank == 0);
130: stag->lastRank[0] = (PetscBool)(rank == size-1);
132: /* Local (unghosted) numbers of entries */
133: stag->entriesPerElement = stag->dof[0] + stag->dof[1];
134: switch (stag->boundaryType[0]) {
135: case DM_BOUNDARY_NONE:
136: case DM_BOUNDARY_GHOSTED: stag->entries = stag->n[0] * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0); break;
137: case DM_BOUNDARY_PERIODIC: stag->entries = stag->n[0] * stag->entriesPerElement; break;
138: default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
139: }
141: /* Starting element */
142: stag->start[0] = 0;
143: for (j=0; j<stag->rank[0]; ++j) stag->start[0] += stag->l[0][j];
145: /* Local/ghosted size and starting element */
146: switch (stag->boundaryType[0]) {
147: case DM_BOUNDARY_NONE :
148: switch (stag->stencilType) {
149: case DMSTAG_STENCIL_NONE : /* Only dummy cells on the right */
150: stag->startGhost[0] = stag->start[0];
151: stag->nGhost[0] = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
152: break;
153: case DMSTAG_STENCIL_STAR :
154: case DMSTAG_STENCIL_BOX :
155: stag->startGhost[0] = stag->firstRank[0] ? stag->start[0]: stag->start[0] - stag->stencilWidth;
156: stag->nGhost[0] = stag->n[0];
157: stag->nGhost[0] += stag->firstRank[0] ? 0 : stag->stencilWidth;
158: stag->nGhost[0] += stag->lastRank[0] ? 1 : stag->stencilWidth;
159: break;
160: default :
161: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
162: }
163: break;
164: case DM_BOUNDARY_GHOSTED:
165: switch (stag->stencilType) {
166: case DMSTAG_STENCIL_NONE :
167: stag->startGhost[0] = stag->start[0];
168: stag->nGhost[0] = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
169: break;
170: case DMSTAG_STENCIL_STAR :
171: case DMSTAG_STENCIL_BOX :
172: stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
173: stag->nGhost[0] = stag->n[0] + 2*stag->stencilWidth + (stag->lastRank[0] && stag->stencilWidth == 0 ? 1 : 0);
174: break;
175: default :
176: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
177: }
178: break;
179: case DM_BOUNDARY_PERIODIC:
180: switch (stag->stencilType) {
181: case DMSTAG_STENCIL_NONE :
182: stag->startGhost[0] = stag->start[0];
183: stag->nGhost[0] = stag->n[0];
184: break;
185: case DMSTAG_STENCIL_STAR :
186: case DMSTAG_STENCIL_BOX :
187: stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
188: stag->nGhost[0] = stag->n[0] + 2*stag->stencilWidth;
189: break;
190: default :
191: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unrecognized ghost stencil type %d",stag->stencilType);
192: }
193: break;
194: default :
195: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
196: }
198: /* Total size of ghosted/local representation */
199: stag->entriesGhost = stag->nGhost[0]*stag->entriesPerElement;
201: /* Define neighbors */
202: PetscMalloc1(3,&stag->neighbors);
203: if (stag->firstRank[0]) {
204: switch (stag->boundaryType[0]) {
205: case DM_BOUNDARY_GHOSTED:
206: case DM_BOUNDARY_NONE: stag->neighbors[0] = -1; break;
207: case DM_BOUNDARY_PERIODIC: stag->neighbors[0] = stag->nRanks[0]-1; break;
208: default : SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
209: }
210: } else {
211: stag->neighbors[0] = stag->rank[0]-1;
212: }
213: stag->neighbors[1] = stag->rank[0];
214: if (stag->lastRank[0]) {
215: switch (stag->boundaryType[0]) {
216: case DM_BOUNDARY_GHOSTED:
217: case DM_BOUNDARY_NONE: stag->neighbors[2] = -1; break;
218: case DM_BOUNDARY_PERIODIC: stag->neighbors[2] = 0; break;
219: default : SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
220: }
221: } else {
222: stag->neighbors[2] = stag->rank[0]+1;
223: }
225: if (stag->n[0] < stag->stencilWidth) {
226: SETERRQ(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);
227: }
229: /* Create global->local VecScatter and ISLocalToGlobalMapping */
230: {
231: PetscInt *idxLocal,*idxGlobal,*idxGlobalAll;
232: PetscInt i,iLocal,d,entriesToTransferTotal,ghostOffsetStart,ghostOffsetEnd,nNonDummyGhost;
233: IS isLocal,isGlobal;
235: /* The offset on the right (may not be equal to the stencil width, as we
236: always have at least one ghost element, to account for the boundary
237: point, and may with ghosted boundaries), and the number of non-dummy ghost elements */
238: ghostOffsetStart = stag->start[0] - stag->startGhost[0];
239: ghostOffsetEnd = stag->startGhost[0]+stag->nGhost[0] - (stag->start[0]+stag->n[0]);
240: nNonDummyGhost = stag->nGhost[0] - (stag->lastRank[0] ? ghostOffsetEnd : 0) - (stag->firstRank[0] ? ghostOffsetStart : 0);
242: /* Compute the number of non-dummy entries in the local representation
243: This is equal to the number of non-dummy elements in the local (ghosted) representation,
244: plus some extra entries on the right boundary on the last rank*/
245: switch (stag->boundaryType[0]) {
246: case DM_BOUNDARY_GHOSTED:
247: case DM_BOUNDARY_NONE:
248: entriesToTransferTotal = nNonDummyGhost * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0);
249: break;
250: case DM_BOUNDARY_PERIODIC:
251: entriesToTransferTotal = stag->entriesGhost; /* No dummy points */
252: break;
253: default :
254: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
255: }
257: PetscMalloc1(entriesToTransferTotal,&idxLocal);
258: PetscMalloc1(entriesToTransferTotal,&idxGlobal);
259: PetscMalloc1(stag->entriesGhost,&idxGlobalAll);
260: if (stag->boundaryType[0] == DM_BOUNDARY_NONE) {
261: PetscInt count = 0,countAll = 0;
262: /* Left ghost points and native points */
263: for (i=stag->startGhost[0], iLocal=0; iLocal<nNonDummyGhost; ++i,++iLocal) {
264: for (d=0; d<stag->entriesPerElement; ++d,++count,++countAll) {
265: idxLocal [count] = iLocal * stag->entriesPerElement + d;
266: idxGlobal[count] = i * stag->entriesPerElement + d;
267: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
268: }
269: }
270: /* Ghost points on the right
271: Special case for last (partial dummy) element on the last rank */
272: if (stag->lastRank[0]) {
273: i = stag->N[0];
274: iLocal = (stag->nGhost[0]-ghostOffsetEnd);
275: /* Only vertex (0-cell) dofs in global representation */
276: for (d=0; d<stag->dof[0]; ++d,++count,++countAll) {
277: idxGlobal[count] = i * stag->entriesPerElement + d;
278: idxLocal [count] = iLocal * stag->entriesPerElement + d;
279: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
280: }
281: for (d=stag->dof[0]; d<stag->entriesPerElement; ++d,++countAll) { /* Additional dummy entries */
282: idxGlobalAll[countAll] = -1;
283: }
284: }
285: } else if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC) {
286: PetscInt count = 0,iLocal = 0; /* No dummy points, so idxGlobal and idxGlobalAll are identical */
287: const PetscInt iMin = stag->firstRank[0] ? stag->start[0] : stag->startGhost[0];
288: const PetscInt iMax = stag->lastRank[0] ? stag->startGhost[0] + stag->nGhost[0] - stag->stencilWidth : stag->startGhost[0] + stag->nGhost[0];
289: /* Ghost points on the left */
290: if (stag->firstRank[0]) {
291: for (i=stag->N[0]-stag->stencilWidth; iLocal<stag->stencilWidth; ++i,++iLocal) {
292: for (d=0; d<stag->entriesPerElement; ++d,++count) {
293: idxGlobal[count] = i * stag->entriesPerElement + d;
294: idxLocal [count] = iLocal * stag->entriesPerElement + d;
295: idxGlobalAll[count] = idxGlobal[count];
296: }
297: }
298: }
299: /* Native points */
300: for (i=iMin; i<iMax; ++i,++iLocal) {
301: for (d=0; d<stag->entriesPerElement; ++d,++count) {
302: idxGlobal[count] = i * stag->entriesPerElement + d;
303: idxLocal [count] = iLocal * stag->entriesPerElement + d;
304: idxGlobalAll[count] = idxGlobal[count];
305: }
306: }
307: /* Ghost points on the right */
308: if (stag->lastRank[0]) {
309: for (i=0; iLocal<stag->nGhost[0]; ++i,++iLocal) {
310: for (d=0; d<stag->entriesPerElement; ++d,++count) {
311: idxGlobal[count] = i * stag->entriesPerElement + d;
312: idxLocal [count] = iLocal * stag->entriesPerElement + d;
313: idxGlobalAll[count] = idxGlobal[count];
314: }
315: }
316: }
317: } else if (stag->boundaryType[0] == DM_BOUNDARY_GHOSTED) {
318: PetscInt count = 0,countAll = 0;
319: /* Dummy elements on the left, on the first rank */
320: if (stag->firstRank[0]) {
321: for (iLocal=0; iLocal<ghostOffsetStart; ++iLocal) {
322: /* Complete elements full of dummy entries */
323: for (d=0; d<stag->entriesPerElement; ++d,++countAll) {
324: idxGlobalAll[countAll] = -1;
325: }
326: }
327: i = 0; /* nonDummy entries start with global entry 0 */
328: } else {
329: /* nonDummy entries start as usual */
330: i = stag->startGhost[0];
331: iLocal = 0;
332: }
334: /* non-Dummy entries */
335: {
336: PetscInt iLocalNonDummyMax = stag->firstRank[0] ? nNonDummyGhost + ghostOffsetStart : nNonDummyGhost;
337: for (; iLocal<iLocalNonDummyMax; ++i,++iLocal) {
338: for (d=0; d<stag->entriesPerElement; ++d,++count,++countAll) {
339: idxLocal [count] = iLocal * stag->entriesPerElement + d;
340: idxGlobal[count] = i * stag->entriesPerElement + d;
341: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
342: }
343: }
344: }
346: /* (partial) dummy elements on the right, on the last rank */
347: if (stag->lastRank[0]) {
348: /* First one is partial dummy */
349: i = stag->N[0];
350: iLocal = (stag->nGhost[0]-ghostOffsetEnd);
351: for (d=0; d<stag->dof[0]; ++d,++count,++countAll) { /* Only vertex (0-cell) dofs in global representation */
352: idxLocal [count] = iLocal * stag->entriesPerElement + d;
353: idxGlobal[count] = i * stag->entriesPerElement + d;
354: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
355: }
356: for (d=stag->dof[0]; d<stag->entriesPerElement; ++d,++countAll) { /* Additional dummy entries */
357: idxGlobalAll[countAll] = -1;
358: }
359: for (iLocal = stag->nGhost[0] - ghostOffsetEnd + 1; iLocal < stag->nGhost[0]; ++iLocal) {
360: /* Additional dummy elements */
361: for (d=0; d<stag->entriesPerElement; ++d,++countAll) {
362: idxGlobalAll[countAll] = -1;
363: }
364: }
365: }
366: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported x boundary type %s",DMBoundaryTypes[stag->boundaryType[0]]);
368: /* Create Local IS (transferring pointer ownership) */
369: ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxLocal,PETSC_OWN_POINTER,&isLocal);
371: /* Create Global IS (transferring pointer ownership) */
372: ISCreateGeneral(PetscObjectComm((PetscObject)dm),entriesToTransferTotal,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
374: /* Create stag->gtol, which doesn't include dummy entries */
375: {
376: Vec local,global;
377: VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
378: VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
379: VecScatterCreate(global,isGlobal,local,isLocal,&stag->gtol);
380: VecDestroy(&global);
381: VecDestroy(&local);
382: }
384: /* In special cases, create a dedicated injective local-to-global map */
385: if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) {
386: DMStagPopulateLocalToGlobalInjective(dm);
387: }
389: /* Destroy ISs */
390: ISDestroy(&isLocal);
391: ISDestroy(&isGlobal);
393: /* Create local-to-global map (transferring pointer ownership) */
394: ISLocalToGlobalMappingCreate(comm,1,stag->entriesGhost,idxGlobalAll,PETSC_OWN_POINTER,&dm->ltogmap);
395: PetscLogObjectParent((PetscObject)dm,(PetscObject)dm->ltogmap);
396: }
398: /* Precompute location offsets */
399: DMStagComputeLocationOffsets_1d(dm);
401: /* View from Options */
402: DMViewFromOptions(dm,NULL,"-dm_view");
404: return 0;
405: }
407: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM dm)
408: {
409: DM_Stag * const stag = (DM_Stag*)dm->data;
410: const PetscInt epe = stag->entriesPerElement;
412: PetscMalloc1(DMSTAG_NUMBER_LOCATIONS,&stag->locationOffsets);
413: stag->locationOffsets[DMSTAG_LEFT] = 0;
414: stag->locationOffsets[DMSTAG_ELEMENT] = stag->locationOffsets[DMSTAG_LEFT] + stag->dof[0];
415: stag->locationOffsets[DMSTAG_RIGHT] = stag->locationOffsets[DMSTAG_LEFT] + epe;
416: return 0;
417: }
419: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_1d(DM dm)
420: {
421: DM_Stag * const stag = (DM_Stag*)dm->data;
422: PetscInt *idxLocal,*idxGlobal;
423: PetscInt i,iLocal,d,count;
424: IS isLocal,isGlobal;
426: PetscMalloc1(stag->entries,&idxLocal);
427: PetscMalloc1(stag->entries,&idxGlobal);
428: count = 0;
429: iLocal = stag->start[0]-stag->startGhost[0];
430: for (i=stag->start[0]; i<stag->start[0]+stag->n[0]; ++i,++iLocal) {
431: for (d=0; d<stag->entriesPerElement; ++d,++count) {
432: idxGlobal[count] = i * stag->entriesPerElement + d;
433: idxLocal [count] = iLocal * stag->entriesPerElement + d;
434: }
435: }
436: if (stag->lastRank[0] && stag->boundaryType[0] != DM_BOUNDARY_PERIODIC) {
437: i = stag->start[0]+stag->n[0];
438: iLocal = stag->start[0]-stag->startGhost[0] + stag->n[0];
439: for (d=0; d<stag->dof[0]; ++d,++count) {
440: idxGlobal[count] = i * stag->entriesPerElement + d;
441: idxLocal [count] = iLocal * stag->entriesPerElement + d;
442: }
443: }
444: ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxLocal,PETSC_OWN_POINTER,&isLocal);
445: ISCreateGeneral(PetscObjectComm((PetscObject)dm),stag->entries,idxGlobal,PETSC_OWN_POINTER,&isGlobal);
446: {
447: Vec local,global;
448: VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),1,stag->entries,PETSC_DECIDE,NULL,&global);
449: VecCreateSeqWithArray(PETSC_COMM_SELF,stag->entriesPerElement,stag->entriesGhost,NULL,&local);
450: VecScatterCreate(local,isLocal,global,isGlobal,&stag->ltog_injective);
451: VecDestroy(&global);
452: VecDestroy(&local);
453: }
454: ISDestroy(&isLocal);
455: ISDestroy(&isGlobal);
456: return 0;
457: }
459: PETSC_INTERN PetscErrorCode DMCreateMatrix_Stag_1D_AIJ_Assemble(DM dm,Mat A)
460: {
461: DMStagStencilType stencil_type;
462: PetscInt dof[2],start,n,n_extra,stencil_width,N,epe;
463: DMBoundaryType boundary_type_x;
465: DMStagGetDOF(dm,&dof[0],&dof[1],NULL,NULL);
466: DMStagGetStencilType(dm,&stencil_type);
467: DMStagGetStencilWidth(dm,&stencil_width);
468: DMStagGetCorners(dm,&start,NULL,NULL,&n,NULL,NULL,&n_extra,NULL,NULL);
469: DMStagGetGlobalSizes(dm,&N,NULL,NULL);
470: DMStagGetEntriesPerElement(dm,&epe);
471: DMStagGetBoundaryTypes(dm,&boundary_type_x,NULL,NULL);
472: if (stencil_type == DMSTAG_STENCIL_NONE) {
473: /* Couple all DOF at each location to each other */
474: DMStagStencil *row_vertex,*row_element;
476: PetscMalloc1(dof[0],&row_vertex);
477: for (PetscInt c=0; c<dof[0]; ++c) {
478: row_vertex[c].loc = DMSTAG_LEFT;
479: row_vertex[c].c = c;
480: }
482: PetscMalloc1(dof[1],&row_element);
483: for (PetscInt c=0; c<dof[1]; ++c) {
484: row_element[c].loc = DMSTAG_ELEMENT;
485: row_element[c].c = c;
486: }
488: for (PetscInt e=start; e<start+n+n_extra; ++e) {
489: {
490: for (PetscInt c=0; c<dof[0]; ++c){
491: row_vertex[c].i = e;
492: }
493: DMStagMatSetValuesStencil(dm,A,dof[0],row_vertex,dof[0],row_vertex,NULL,INSERT_VALUES);
494: }
495: if (e < N) {
496: for (PetscInt c=0; c<dof[1]; ++c) {
497: row_element[c].i = e;
498: }
499: DMStagMatSetValuesStencil(dm,A,dof[1],row_element,dof[1],row_element,NULL,INSERT_VALUES);
500: }
501: }
502: PetscFree(row_vertex);
503: PetscFree(row_element);
504: } else if (stencil_type == DMSTAG_STENCIL_STAR || stencil_type == DMSTAG_STENCIL_BOX) {
505: DMStagStencil *col,*row;
507: PetscMalloc1(epe,&row);
508: {
509: PetscInt nrows = 0;
510: for (PetscInt c=0; c<dof[0]; ++c) {
511: row[nrows].c = c;
512: row[nrows].loc = DMSTAG_LEFT;
513: ++nrows;
514: }
515: for (PetscInt c=0; c<dof[1]; ++c) {
516: row[nrows].c = c;
517: row[nrows].loc = DMSTAG_ELEMENT;
518: ++nrows;
519: }
520: }
521: PetscMalloc1(epe,&col);
522: {
523: PetscInt ncols = 0;
524: for (PetscInt c=0; c<dof[0]; ++c) {
525: col[ncols].c = c;
526: col[ncols].loc = DMSTAG_LEFT;
527: ++ncols;
528: }
529: for (PetscInt c=0; c<dof[1]; ++c) {
530: col[ncols].c = c;
531: col[ncols].loc = DMSTAG_ELEMENT;
532: ++ncols;
533: }
534: }
535: for (PetscInt e=start; e<start+n+n_extra; ++e) {
536: for (PetscInt i=0; i<epe; ++i) {
537: row[i].i = e;
538: }
539: for (PetscInt offset = -stencil_width; offset<=stencil_width; ++offset) {
540: const PetscInt e_offset = e + offset;
542: /* Only set values corresponding to elements which can have non-dummy entries,
543: meaning those that map to unknowns in the global representation. In the periodic
544: case, this is the entire stencil, but in all other cases, only includes a single
545: "extra" element which is partially outside the physical domain (those points in the
546: global representation */
547: if (boundary_type_x == DM_BOUNDARY_PERIODIC || (e_offset < N+1 && e_offset >= 0)) {
548: for (PetscInt i=0; i<epe; ++i) {
549: col[i].i = e_offset;
550: }
551: DMStagMatSetValuesStencil(dm,A,epe,row,epe,col,NULL,INSERT_VALUES);
552: }
553: }
554: }
555: PetscFree(row);
556: PetscFree(col);
557: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported stencil type %s",DMStagStencilTypes[stencil_type]);
558: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
559: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
560: return 0;
561: }