Actual source code: stagda.c
petsc-3.12.5 2020-03-29
1: /* Routines to convert between a (subset of) DMStag and DMDA */
3: #include <petscdmda.h>
4: #include <petsc/private/dmstagimpl.h>
5: #include <petscdmdatypes.h>
7: static PetscErrorCode DMStagCreateCompatibleDMDA(DM dm,DMStagStencilLocation loc,PetscInt c,DM *dmda)
8: {
9: PetscErrorCode ierr;
10: DM_Stag * const stag = (DM_Stag*) dm->data;
11: PetscInt dim,i,j,stencilWidth,dof,N[DMSTAG_MAX_DIM];
12: DMDAStencilType stencilType;
13: PetscInt *l[DMSTAG_MAX_DIM];
17: DMGetDimension(dm,&dim);
19: /* Create grid decomposition (to be adjusted later) */
20: for (i=0; i<dim; ++i) {
21: PetscMalloc1(stag->nRanks[i],&l[i]);
22: for (j=0; j<stag->nRanks[i]; ++j) l[i][j] = stag->l[i][j];
23: N[i] = stag->N[i];
24: }
26: /* dof */
27: dof = c < 0 ? -c : 1;
29: /* Determine/adjust sizes */
30: switch (loc) {
31: case DMSTAG_ELEMENT:
32: break;
33: case DMSTAG_LEFT:
34: case DMSTAG_RIGHT:
35: #if defined(PETSC_USE_DEBUG)
36: if (dim<1) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Incompatible dim (%d) and loc(%s) combination",dim,DMStagStencilLocations[loc]);
37: #endif
38: l[0][stag->nRanks[0]-1] += 1; /* extra vertex in direction 0 on last rank in dimension 0 */
39: N[0] += 1;
40: break;
41: case DMSTAG_UP:
42: case DMSTAG_DOWN:
43: #if defined(PETSC_USE_DEBUG)
44: if (dim < 2) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Incompatible dim (%d) and loc(%s) combination",dim,DMStagStencilLocations[loc]);
45: #endif
46: l[1][stag->nRanks[1]-1] += 1; /* extra vertex in direction 1 on last rank in dimension 1 */
47: N[1] += 1;
48: break;
49: case DMSTAG_BACK:
50: case DMSTAG_FRONT:
51: #if defined(PETSC_USE_DEBUG)
52: if (dim < 3) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Incompatible dim (%d) and loc(%s) combination",dim,DMStagStencilLocations[loc]);
53: #endif
54: l[2][stag->nRanks[2]-1] += 1; /* extra vertex in direction 2 on last rank in dimension 2 */
55: N[2] += 1;
56: break;
57: case DMSTAG_DOWN_LEFT :
58: case DMSTAG_DOWN_RIGHT :
59: case DMSTAG_UP_LEFT :
60: case DMSTAG_UP_RIGHT :
61: #if defined(PETSC_USE_DEBUG)
62: if (dim < 2) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Incompatible dim (%d) and loc(%s) combination",dim,DMStagStencilLocations[loc]);
63: #endif
64: for (i=0; i<2; ++i) { /* extra vertex in direction i on last rank in dimension i = 0,1 */
65: l[i][stag->nRanks[i]-1] += 1;
66: N[i] += 1;
67: }
68: break;
69: case DMSTAG_BACK_LEFT:
70: case DMSTAG_BACK_RIGHT:
71: case DMSTAG_FRONT_LEFT:
72: case DMSTAG_FRONT_RIGHT:
73: #if defined(PETSC_USE_DEBUG)
74: if (dim < 3) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Incompatible dim (%d) and loc(%s) combination",dim,DMStagStencilLocations[loc]);
75: #endif
76: for (i=0; i<3; i+=2) { /* extra vertex in direction i on last rank in dimension i = 0,2 */
77: l[i][stag->nRanks[i]-1] += 1;
78: N[i] += 1;
79: }
80: break;
81: case DMSTAG_BACK_DOWN:
82: case DMSTAG_BACK_UP:
83: case DMSTAG_FRONT_DOWN:
84: case DMSTAG_FRONT_UP:
85: #if defined(PETSC_USE_DEBUG)
86: if (dim < 3) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Incompatible dim (%d) and loc(%s) combination",dim,DMStagStencilLocations[loc]);
87: #endif
88: for (i=1; i<3; ++i) { /* extra vertex in direction i on last rank in dimension i = 1,2 */
89: l[i][stag->nRanks[i]-1] += 1;
90: N[i] += 1;
91: }
92: break;
93: case DMSTAG_BACK_DOWN_LEFT:
94: case DMSTAG_BACK_DOWN_RIGHT:
95: case DMSTAG_BACK_UP_LEFT:
96: case DMSTAG_BACK_UP_RIGHT:
97: case DMSTAG_FRONT_DOWN_LEFT:
98: case DMSTAG_FRONT_DOWN_RIGHT:
99: case DMSTAG_FRONT_UP_LEFT:
100: case DMSTAG_FRONT_UP_RIGHT:
101: #if defined(PETSC_USE_DEBUG)
102: if (dim < 3) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Incompatible dim (%d) and loc(%s) combination",dim,DMStagStencilLocations[loc]);
103: #endif
104: for (i=0; i<3; ++i) { /* extra vertex in direction i on last rank in dimension i = 0,1,2 */
105: l[i][stag->nRanks[i]-1] += 1;
106: N[i] += 1;
107: }
108: break;
109: break;
110: default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for location %s",DMStagStencilLocations[loc]);
111: }
113: /* Use the same stencil type */
114: switch (stag->stencilType) {
115: case DMSTAG_STENCIL_STAR: stencilType = DMDA_STENCIL_STAR; stencilWidth = stag->stencilWidth; break;
116: case DMSTAG_STENCIL_BOX : stencilType = DMDA_STENCIL_BOX ; stencilWidth = stag->stencilWidth; break;
117: default: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported Stencil Type %d",stag->stencilType);
118: }
120: /* Create DMDA, using same boundary type */
121: switch (dim) {
122: case 1:
123: DMDACreate1d(PetscObjectComm((PetscObject)dm),stag->boundaryType[0],N[0],dof,stencilWidth,l[0],dmda);
124: break;
125: case 2:
126: DMDACreate2d(PetscObjectComm((PetscObject)dm),stag->boundaryType[0],stag->boundaryType[1],stencilType,N[0],N[1],stag->nRanks[0],stag->nRanks[1],dof,stencilWidth,l[0],l[1],dmda);
127: break;
128: case 3:
129: DMDACreate3d(PetscObjectComm((PetscObject)dm),stag->boundaryType[0],stag->boundaryType[1],stag->boundaryType[2],stencilType,N[0],N[1],N[2],stag->nRanks[0],stag->nRanks[1],stag->nRanks[2],dof,stencilWidth,l[0],l[1],l[2],dmda);
130: break;
131: default: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not implemented for dim %d",dim);
132: }
133: for (i=0; i<dim; ++i) {
134: PetscFree(l[i]);
135: }
136: return(0);
137: }
139: /*
140: Helper function to get the number of extra points in a DMDA representation for a given canonical location.
141: */
142: static PetscErrorCode DMStagDMDAGetExtraPoints(DM dm,DMStagStencilLocation locCanonical,PetscInt *extraPoint)
143: {
145: PetscInt dim,d,nExtra[DMSTAG_MAX_DIM];
149: DMGetDimension(dm,&dim);
150: if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
151: DMStagGetCorners(dm,NULL,NULL,NULL,NULL,NULL,NULL,&nExtra[0],&nExtra[1],&nExtra[2]);
152: for(d=0; d<dim; ++d) extraPoint[d] = 0;
153: switch (locCanonical) {
154: case DMSTAG_ELEMENT:
155: break; /* no extra points */
156: case DMSTAG_LEFT:
157: extraPoint[0] = nExtra[0]; break; /* only extra point in x */
158: case DMSTAG_DOWN:
159: extraPoint[1] = nExtra[1]; break; /* only extra point in y */
160: case DMSTAG_BACK:
161: extraPoint[2] = nExtra[2]; break; /* only extra point in z */
162: case DMSTAG_DOWN_LEFT:
163: extraPoint[0] = nExtra[0]; extraPoint[1] = nExtra[1]; break; /* extra point in both x and y */
164: case DMSTAG_BACK_LEFT:
165: extraPoint[0] = nExtra[0]; extraPoint[2] = nExtra[2]; break; /* extra point in both x and z */
166: case DMSTAG_BACK_DOWN:
167: extraPoint[1] = nExtra[1]; extraPoint[2] = nExtra[2]; break; /* extra point in both y and z */
168: case DMSTAG_BACK_DOWN_LEFT:
169: extraPoint[0] = nExtra[0]; extraPoint[1] = nExtra[1]; extraPoint[2] = nExtra[2]; break; /* extra points in x,y,z */
170: default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for location (perhaps not canonical) %s",DMStagStencilLocations[locCanonical]);
171: }
172: return(0);
173: }
175: /*
176: Function much like DMStagMigrateVec(), but which accepts an additional position argument to disambiguate which
177: type of DMDA to migrate to.
178: */
180: static PetscErrorCode DMStagMigrateVecDMDA(DM dm,Vec vec,DMStagStencilLocation loc,PetscInt c,DM dmTo,Vec vecTo)
181: {
183: PetscInt i,j,k,d,dim,dof,dofToMax,start[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],extraPoint[DMSTAG_MAX_DIM];
184: Vec vecLocal;
191: DMGetDimension(dm,&dim);
192: DMDAGetDof(dmTo,&dofToMax);
193: #if defined(PETSC_USE_DEBUG)
194: if (-c > dofToMax) SETERRQ1(PetscObjectComm((PetscObject)dmTo),PETSC_ERR_ARG_OUTOFRANGE,"Invalid negative component value. Must be >= -%D",dofToMax);
195: #endif
196: DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],NULL,NULL,NULL);
197: DMStagDMDAGetExtraPoints(dm,loc,extraPoint);
198: DMStagGetLocationDOF(dm,loc,&dof);
199: DMGetLocalVector(dm,&vecLocal);
200: DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);
201: DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);
202: if (dim == 1) {
203: PetscScalar **arrTo;
204: DMDAVecGetArrayDOF(dmTo,vecTo,&arrTo);
205: if (c < 0) {
206: const PetscInt dofTo = -c;
207: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
208: for (d=0; d<PetscMin(dof,dofTo); ++d) {
209: DMStagStencil pos;
210: pos.i = i; pos.loc = loc; pos.c = d;
211: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[i][d]);
212: }
213: for (;d<dofTo; ++d) {
214: arrTo[i][d] = 0.0; /* Pad extra dof with zeros */
215: }
216: }
217: } else {
218: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
219: DMStagStencil pos;
220: pos.i = i; pos.loc = loc; pos.c = c;
221: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[i][0]);
222: }
223: }
224: DMDAVecRestoreArrayDOF(dmTo,vecTo,&arrTo);
225: } else if (dim == 2) {
226: PetscScalar ***arrTo;
227: DMDAVecGetArrayDOF(dmTo,vecTo,&arrTo);
228: if (c < 0) {
229: const PetscInt dofTo = -c;
230: for (j=start[1]; j<start[1] + n[1] + extraPoint[1]; ++j) {
231: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
232: for (d=0; d<PetscMin(dof,dofTo); ++d) {
233: DMStagStencil pos;
234: pos.i = i; pos.j = j; pos.loc = loc; pos.c = d;
235: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[j][i][d]);
236: }
237: for (;d<dofTo; ++d) {
238: arrTo[j][i][d] = 0.0; /* Pad extra dof with zeros */
239: }
240: }
241: }
242: } else {
243: for (j=start[1]; j<start[1] + n[1] + extraPoint[1]; ++j) {
244: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
245: DMStagStencil pos;
246: pos.i = i; pos.j = j; pos.loc = loc; pos.c = c;
247: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[j][i][0]);
248: }
249: }
250: }
251: DMDAVecRestoreArrayDOF(dmTo,vecTo,&arrTo);
252: } else if (dim == 3) {
253: PetscScalar ****arrTo;
254: DMDAVecGetArrayDOF(dmTo,vecTo,&arrTo);
255: if (c < 0) {
256: const PetscInt dofTo = -c;
257: for (k=start[2]; k<start[2] + n[2] + extraPoint[2]; ++k) {
258: for (j=start[1]; j<start[1] + n[1] + extraPoint[1]; ++j) {
259: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
260: for (d=0; d<PetscMin(dof,dofTo); ++d) {
261: DMStagStencil pos;
262: pos.i = i; pos.j = j; pos.k = k; pos.loc = loc; pos.c = d;
263: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[k][j][i][d]);
264: }
265: for (;d<dofTo; ++d) {
266: arrTo[k][j][i][d] = 0.0; /* Pad extra dof with zeros */
267: }
268: }
269: }
270: }
271: } else {
272: for (k=start[2]; k<start[2] + n[2] + extraPoint[2]; ++k) {
273: for (j=start[1]; j<start[1] + n[1] + extraPoint[1]; ++j) {
274: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
275: DMStagStencil pos;
276: pos.i = i; pos.j = j; pos.k = k; pos.loc = loc; pos.c = c;
277: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[k][j][i][0]);
278: }
279: }
280: }
281: }
282: DMDAVecRestoreArrayDOF(dmTo,vecTo,&arrTo);
283: } else SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported dimension %d",dim);
284: DMRestoreLocalVector(dm,&vecLocal);
285: return(0);
286: }
288: /* Transfer coordinates from a DMStag to a DMDA, specifying which location */
289: static PetscErrorCode DMStagTransferCoordinatesToDMDA(DM dmstag,DMStagStencilLocation loc,DM dmda)
290: {
292: PetscInt dim,start[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],extraPoint[DMSTAG_MAX_DIM],d;
293: DM dmstagCoord,dmdaCoord;
294: DMType dmstagCoordType;
295: Vec stagCoord,daCoord;
296: PetscBool daCoordIsStag,daCoordIsProduct;
301: DMGetDimension(dmstag,&dim);
302: DMGetCoordinateDM(dmstag,&dmstagCoord);
303: DMGetCoordinatesLocal(dmstag,&stagCoord); /* Note local */
304: DMGetCoordinateDM(dmda,&dmdaCoord);
305: daCoord = NULL;
306: DMGetCoordinates(dmda,&daCoord);
307: if (!daCoord) {
308: DMCreateGlobalVector(dmdaCoord,&daCoord);
309: DMSetCoordinates(dmda,daCoord);
310: VecDestroy(&daCoord);
311: DMGetCoordinates(dmda,&daCoord);
312: }
313: DMGetType(dmstagCoord,&dmstagCoordType);
314: PetscStrcmp(dmstagCoordType,DMSTAG,&daCoordIsStag);
315: PetscStrcmp(dmstagCoordType,DMPRODUCT,&daCoordIsProduct);
316: DMStagGetCorners(dmstag,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],NULL,NULL,NULL);
317: DMStagDMDAGetExtraPoints(dmstag,loc,extraPoint);
318: if (dim == 1) {
319: PetscInt ex;
320: PetscScalar **cArrDa;
321: DMDAVecGetArrayDOF(dmdaCoord,daCoord,&cArrDa);
322: if (daCoordIsStag) {
323: PetscInt slot;
324: PetscScalar **cArrStag;
325: DMStagGetLocationSlot(dmstagCoord,loc,0,&slot);
326: DMStagVecGetArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
327: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
328: cArrDa[ex][0] = cArrStag[ex][slot];
329: }
330: DMStagVecRestoreArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
331: } else if (daCoordIsProduct) {
332: PetscScalar **cArrX;
333: DMStagGet1dCoordinateArraysDOFRead(dmstag,&cArrX,NULL,NULL);
334: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
335: cArrDa[ex][0] = cArrX[ex][0];
336: }
337: DMStagRestore1dCoordinateArraysDOFRead(dmstag,&cArrX,NULL,NULL);
338: } else SETERRQ(PetscObjectComm((PetscObject)dmstag),PETSC_ERR_SUP,"Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
339: DMDAVecRestoreArrayDOF(dmdaCoord,daCoord,&cArrDa);
340: } else if (dim == 2) {
341: PetscInt ex,ey;
342: PetscScalar ***cArrDa;
343: DMDAVecGetArrayDOF(dmdaCoord,daCoord,&cArrDa);
344: if (daCoordIsStag) {
345: PetscInt slot;
346: PetscScalar ***cArrStag;
347: DMStagGetLocationSlot(dmstagCoord,loc,0,&slot);
348: DMStagVecGetArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
349: for (ey=start[1]; ey<start[1] + n[1] + extraPoint[1]; ++ey) {
350: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
351: for (d=0; d<2; ++d) {
352: cArrDa[ey][ex][d] = cArrStag[ey][ex][slot+d];
353: }
354: }
355: }
356: DMStagVecRestoreArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
357: } else if (daCoordIsProduct) {
358: PetscScalar **cArrX,**cArrY;
359: DMStagGet1dCoordinateArraysDOFRead(dmstag,&cArrX,&cArrY,NULL);
360: for (ey=start[1]; ey<start[1] + n[1] + extraPoint[1]; ++ey) {
361: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
362: cArrDa[ey][ex][0] = cArrX[ex][0];
363: cArrDa[ey][ex][1] = cArrY[ey][0];
364: }
365: }
366: DMStagRestore1dCoordinateArraysDOFRead(dmstag,&cArrX,&cArrY,NULL);
367: } else SETERRQ(PetscObjectComm((PetscObject)dmstag),PETSC_ERR_SUP,"Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
368: DMDAVecRestoreArrayDOF(dmdaCoord,daCoord,&cArrDa);
369: } else if (dim == 3) {
370: PetscInt ex,ey,ez;
371: PetscScalar ****cArrDa;
372: DMDAVecGetArrayDOF(dmdaCoord,daCoord,&cArrDa);
373: if (daCoordIsStag) {
374: PetscInt slot;
375: PetscScalar ****cArrStag;
376: DMStagGetLocationSlot(dmstagCoord,loc,0,&slot);
377: DMStagVecGetArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
378: for (ez=start[2]; ez<start[2] + n[2] + extraPoint[2]; ++ez) {
379: for (ey=start[1]; ey<start[1] + n[1] + extraPoint[1]; ++ey) {
380: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
381: for (d=0; d<3; ++d) {
382: cArrDa[ez][ey][ex][d] = cArrStag[ez][ey][ex][slot+d];
383: }
384: }
385: }
386: }
387: DMStagVecRestoreArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
388: } else if (daCoordIsProduct) {
389: PetscScalar **cArrX,**cArrY,**cArrZ;
390: DMStagGet1dCoordinateArraysDOFRead(dmstag,&cArrX,&cArrY,&cArrZ);
391: for (ez=start[2]; ez<start[2] + n[2] + extraPoint[2]; ++ez) {
392: for (ey=start[1]; ey<start[1] + n[1] + extraPoint[1]; ++ey) {
393: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
394: cArrDa[ez][ey][ex][0] = cArrX[ex][0];
395: cArrDa[ez][ey][ex][1] = cArrY[ey][0];
396: cArrDa[ez][ey][ex][2] = cArrZ[ez][0];
397: }
398: }
399: }
400: DMStagRestore1dCoordinateArraysDOFRead(dmstag,&cArrX,&cArrY,&cArrZ);
401: } else SETERRQ(PetscObjectComm((PetscObject)dmstag),PETSC_ERR_SUP,"Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
402: DMDAVecRestoreArrayDOF(dmdaCoord,daCoord,&cArrDa);
403: } else SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported dimension %d",dim);
404: return(0);
405: }
407: /*
408: Convert to a location value with only BACK, DOWN, LEFT, and ELEMENT involved (makes looping easier)
409: */
410: static PetscErrorCode DMStagStencilLocationCanonicalize(DMStagStencilLocation loc,DMStagStencilLocation *locCanonical)
411: {
413: switch (loc) {
414: case DMSTAG_ELEMENT:
415: *locCanonical = DMSTAG_ELEMENT;
416: break;
417: case DMSTAG_LEFT:
418: case DMSTAG_RIGHT:
419: *locCanonical = DMSTAG_LEFT;
420: break;
421: case DMSTAG_DOWN:
422: case DMSTAG_UP:
423: *locCanonical = DMSTAG_DOWN;
424: break;
425: case DMSTAG_BACK:
426: case DMSTAG_FRONT:
427: *locCanonical = DMSTAG_BACK;
428: break;
429: case DMSTAG_DOWN_LEFT :
430: case DMSTAG_DOWN_RIGHT :
431: case DMSTAG_UP_LEFT :
432: case DMSTAG_UP_RIGHT :
433: *locCanonical = DMSTAG_DOWN_LEFT;
434: break;
435: case DMSTAG_BACK_LEFT:
436: case DMSTAG_BACK_RIGHT:
437: case DMSTAG_FRONT_LEFT:
438: case DMSTAG_FRONT_RIGHT:
439: *locCanonical = DMSTAG_BACK_LEFT;
440: break;
441: case DMSTAG_BACK_DOWN:
442: case DMSTAG_BACK_UP:
443: case DMSTAG_FRONT_DOWN:
444: case DMSTAG_FRONT_UP:
445: *locCanonical = DMSTAG_BACK_DOWN;
446: break;
447: case DMSTAG_BACK_DOWN_LEFT:
448: case DMSTAG_BACK_DOWN_RIGHT:
449: case DMSTAG_BACK_UP_LEFT:
450: case DMSTAG_BACK_UP_RIGHT:
451: case DMSTAG_FRONT_DOWN_LEFT:
452: case DMSTAG_FRONT_DOWN_RIGHT:
453: case DMSTAG_FRONT_UP_LEFT:
454: case DMSTAG_FRONT_UP_RIGHT:
455: *locCanonical = DMSTAG_BACK_DOWN_LEFT;
456: break;
457: default :
458: *locCanonical = DMSTAG_NULL_LOCATION;
459: break;
460: }
461: return(0);
462: }
464: /*@C
465: DMStagVecSplitToDMDA - create a DMDA and Vec from a DMStag and Vec
467: Logically Collective
469: High-level helper function which accepts a DMStag, a global vector, and location/dof,
470: and generates a corresponding DMDA and Vec.
472: Input Parameters:
473: + dm - the DMStag object
474: . vec- Vec object associated with dm
475: . loc - which subgrid to extract (see DMStagStencilLocation)
476: - c - which component to extract (see note below)
478: Output Parameters:
479: + pda - the new DMDA
480: - pdavec - the new Vec
482: Notes:
483: If a c value of -k is provided, the first k dof for that position are extracted,
484: padding with zero values if needbe. If a non-negative value is provided, a single
485: dof is extracted.
487: The caller is responsible for destroying the created DMDA and Vec.
489: Level: advanced
491: .seealso: DMSTAG, DMDA, DMStagMigrateVec(), DMStagCreateCompatibleDMStag()
492: @*/
493: PetscErrorCode DMStagVecSplitToDMDA(DM dm,Vec vec,DMStagStencilLocation loc,PetscInt c,DM *pda,Vec *pdavec)
494: {
495: PetscErrorCode ierr;
496: PetscInt dim,locdof;
497: DM da,coordDM;
498: Vec davec;
499: DMStagStencilLocation locCanonical;
504: DMGetDimension(dm,&dim);
505: DMStagGetLocationDOF(dm,loc,&locdof);
506: if (c >= locdof) SETERRQ3(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has %D dof, but component %D requested\n",DMStagStencilLocations[loc],locdof,c);
507: DMStagStencilLocationCanonicalize(loc,&locCanonical);
508: DMStagCreateCompatibleDMDA(dm,locCanonical,c,pda);
509: da = *pda;
510: DMSetUp(*pda);
511: DMGetCoordinateDM(dm,&coordDM);
512: if (coordDM) {
513: DMStagTransferCoordinatesToDMDA(dm,locCanonical,da);
514: }
515: DMCreateGlobalVector(da,pdavec);
516: davec = *pdavec;
517: DMStagMigrateVecDMDA(dm,vec,locCanonical,c,da,davec);
518: return(0);
519: }