Actual source code: stagda.c
petsc-3.11.4 2019-09-28
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: DMStagGetCorners(dm,NULL,NULL,NULL,NULL,NULL,NULL,&nExtra[0],&nExtra[1],&nExtra[2]);
151: for(d=0; d<dim; ++d) extraPoint[d] = 0;
152: switch (locCanonical) {
153: case DMSTAG_ELEMENT:
154: break; /* no extra points */
155: case DMSTAG_LEFT:
156: extraPoint[0] = nExtra[0]; break; /* only extra point in x */
157: case DMSTAG_DOWN:
158: extraPoint[1] = nExtra[1];break; /* only extra point in y */
159: case DMSTAG_BACK:
160: extraPoint[2] = nExtra[2]; break; /* only extra point in z */
161: case DMSTAG_DOWN_LEFT:
162: extraPoint[0] = nExtra[0]; extraPoint[1] = nExtra[1]; break; /* extra point in both x and y */
163: case DMSTAG_BACK_LEFT:
164: extraPoint[0] = nExtra[0]; extraPoint[2] = nExtra[2]; break; /* extra point in both x and z */
165: case DMSTAG_BACK_DOWN:
166: extraPoint[1] = nExtra[1]; extraPoint[2] = nExtra[2]; break; /* extra point in both y and z */
167: case DMSTAG_BACK_DOWN_LEFT:
168: extraPoint[0] = nExtra[0]; extraPoint[1] = nExtra[1]; extraPoint[2] = nExtra[2]; break; /* extra points in x,y,z */
169: default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for location (perhaps not canonical) %s",DMStagStencilLocations[locCanonical]);
170: }
171: return(0);
172: }
174: /*
175: Function much like DMStagMigrateVec(), but which accepts an additional position argument to disambiguate which
176: type of DMDA to migrate to.
177: */
179: static PetscErrorCode DMStagMigrateVecDMDA(DM dm,Vec vec,DMStagStencilLocation loc,PetscInt c,DM dmTo,Vec vecTo)
180: {
182: PetscInt i,j,k,d,dim,dof,dofToMax,start[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],extraPoint[DMSTAG_MAX_DIM];
183: Vec vecLocal;
190: DMGetDimension(dm,&dim);
191: DMDAGetDof(dmTo,&dofToMax);
192: #if defined(PETSC_USE_DEBUG)
193: if (-c > dofToMax) SETERRQ1(PetscObjectComm((PetscObject)dmTo),PETSC_ERR_ARG_OUTOFRANGE,"Invalid negative component value. Must be >= -%D",dofToMax);
194: #endif
195: DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],NULL,NULL,NULL);
196: DMStagDMDAGetExtraPoints(dm,loc,extraPoint);
197: DMStagGetLocationDOF(dm,loc,&dof);
198: DMGetLocalVector(dm,&vecLocal);
199: DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);
200: DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);
201: if (dim == 1) {
202: PetscScalar **arrTo;
203: DMDAVecGetArrayDOF(dmTo,vecTo,&arrTo);
204: if (c < 0) {
205: const PetscInt dofTo = -c;
206: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
207: for (d=0; d<PetscMin(dof,dofTo); ++d) {
208: DMStagStencil pos;
209: pos.i = i; pos.loc = loc; pos.c = d;
210: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[i][d]);
211: }
212: for (;d<dofTo; ++d) {
213: arrTo[i][d] = 0.0; /* Pad extra dof with zeros */
214: }
215: }
216: } else {
217: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
218: DMStagStencil pos;
219: pos.i = i; pos.loc = loc; pos.c = c;
220: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[i][0]);
221: }
222: }
223: DMDAVecRestoreArrayDOF(dmTo,vecTo,&arrTo);
224: } else if (dim == 2) {
225: PetscScalar ***arrTo;
226: DMDAVecGetArrayDOF(dmTo,vecTo,&arrTo);
227: if (c < 0) {
228: const PetscInt dofTo = -c;
229: for (j=start[1]; j<start[1] + n[1] + extraPoint[1]; ++j) {
230: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
231: for (d=0; d<PetscMin(dof,dofTo); ++d) {
232: DMStagStencil pos;
233: pos.i = i; pos.j = j; pos.loc = loc; pos.c = d;
234: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[j][i][d]);
235: }
236: for (;d<dofTo; ++d) {
237: arrTo[j][i][d] = 0.0; /* Pad extra dof with zeros */
238: }
239: }
240: }
241: } else {
242: for (j=start[1]; j<start[1] + n[1] + extraPoint[1]; ++j) {
243: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
244: DMStagStencil pos;
245: pos.i = i; pos.j = j; pos.loc = loc; pos.c = c;
246: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[j][i][0]);
247: }
248: }
249: }
250: DMDAVecRestoreArrayDOF(dmTo,vecTo,&arrTo);
251: } else if (dim == 3) {
252: PetscScalar ****arrTo;
253: DMDAVecGetArrayDOF(dmTo,vecTo,&arrTo);
254: if (c < 0) {
255: const PetscInt dofTo = -c;
256: for (k=start[2]; k<start[2] + n[2] + extraPoint[2]; ++k) {
257: for (j=start[1]; j<start[1] + n[1] + extraPoint[1]; ++j) {
258: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
259: for (d=0; d<PetscMin(dof,dofTo); ++d) {
260: DMStagStencil pos;
261: pos.i = i; pos.j = j; pos.k = k; pos.loc = loc; pos.c = d;
262: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[k][j][i][d]);
263: }
264: for (;d<dofTo; ++d) {
265: arrTo[k][j][i][d] = 0.0; /* Pad extra dof with zeros */
266: }
267: }
268: }
269: }
270: } else {
271: for (k=start[2]; k<start[2] + n[2] + extraPoint[2]; ++k) {
272: for (j=start[1]; j<start[1] + n[1] + extraPoint[1]; ++j) {
273: for (i=start[0]; i<start[0] + n[0] + extraPoint[0]; ++i) {
274: DMStagStencil pos;
275: pos.i = i; pos.j = j; pos.k = k; pos.loc = loc; pos.c = c;
276: DMStagVecGetValuesStencil(dm,vecLocal,1,&pos,&arrTo[k][j][i][0]);
277: }
278: }
279: }
280: }
281: DMDAVecRestoreArrayDOF(dmTo,vecTo,&arrTo);
282: } else SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported dimension %d",dim);
283: DMRestoreLocalVector(dm,&vecLocal);
284: return(0);
285: }
287: /* Transfer coordinates from a DMStag to a DMDA, specifying which location */
288: static PetscErrorCode DMStagTransferCoordinatesToDMDA(DM dmstag,DMStagStencilLocation loc,DM dmda)
289: {
291: PetscInt dim,start[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],extraPoint[DMSTAG_MAX_DIM],d;
292: DM dmstagCoord,dmdaCoord;
293: DMType dmstagCoordType;
294: Vec stagCoord,daCoord;
295: PetscBool daCoordIsStag,daCoordIsProduct;
300: DMGetDimension(dmstag,&dim);
301: DMGetCoordinateDM(dmstag,&dmstagCoord);
302: DMGetCoordinatesLocal(dmstag,&stagCoord); /* Note local */
303: DMGetCoordinateDM(dmda,&dmdaCoord);
304: daCoord = NULL;
305: DMGetCoordinates(dmda,&daCoord);
306: if (!daCoord) {
307: DMCreateGlobalVector(dmdaCoord,&daCoord);
308: DMSetCoordinates(dmda,daCoord);
309: VecDestroy(&daCoord);
310: DMGetCoordinates(dmda,&daCoord);
311: }
312: DMGetType(dmstagCoord,&dmstagCoordType);
313: PetscStrcmp(dmstagCoordType,DMSTAG,&daCoordIsStag);
314: PetscStrcmp(dmstagCoordType,DMPRODUCT,&daCoordIsProduct);
315: DMStagGetCorners(dmstag,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],NULL,NULL,NULL);
316: DMStagDMDAGetExtraPoints(dmstag,loc,extraPoint);
317: if (dim == 1) {
318: PetscInt ex;
319: PetscScalar **cArrDa;
320: DMDAVecGetArrayDOF(dmdaCoord,daCoord,&cArrDa);
321: if (daCoordIsStag) {
322: PetscInt slot;
323: PetscScalar **cArrStag;
324: DMStagGetLocationSlot(dmstagCoord,loc,0,&slot);
325: DMStagVecGetArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
326: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
327: cArrDa[ex][0] = cArrStag[ex][slot];
328: }
329: DMStagVecRestoreArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
330: } else if (daCoordIsProduct) {
331: PetscScalar **cArrX;
332: DMStagGet1dCoordinateArraysDOFRead(dmstag,&cArrX,NULL,NULL);
333: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
334: cArrDa[ex][0] = cArrX[ex][0];
335: }
336: DMStagRestore1dCoordinateArraysDOFRead(dmstag,&cArrX,NULL,NULL);
337: } else SETERRQ(PetscObjectComm((PetscObject)dmstag),PETSC_ERR_SUP,"Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
338: DMDAVecRestoreArrayDOF(dmdaCoord,daCoord,&cArrDa);
339: } else if (dim == 2) {
340: PetscInt ex,ey;
341: PetscScalar ***cArrDa;
342: DMDAVecGetArrayDOF(dmdaCoord,daCoord,&cArrDa);
343: if (daCoordIsStag) {
344: PetscInt slot;
345: PetscScalar ***cArrStag;
346: DMStagGetLocationSlot(dmstagCoord,loc,0,&slot);
347: DMStagVecGetArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
348: for (ey=start[1]; ey<start[1] + n[1] + extraPoint[1]; ++ey) {
349: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
350: for (d=0; d<2; ++d) {
351: cArrDa[ey][ex][d] = cArrStag[ey][ex][slot+d];
352: }
353: }
354: }
355: DMStagVecRestoreArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
356: } else if (daCoordIsProduct) {
357: PetscScalar **cArrX,**cArrY;
358: DMStagGet1dCoordinateArraysDOFRead(dmstag,&cArrX,&cArrY,NULL);
359: for (ey=start[1]; ey<start[1] + n[1] + extraPoint[1]; ++ey) {
360: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
361: cArrDa[ey][ex][0] = cArrX[ex][0];
362: cArrDa[ey][ex][1] = cArrY[ey][0];
363: }
364: }
365: DMStagRestore1dCoordinateArraysDOFRead(dmstag,&cArrX,&cArrY,NULL);
366: } else SETERRQ(PetscObjectComm((PetscObject)dmstag),PETSC_ERR_SUP,"Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
367: DMDAVecRestoreArrayDOF(dmdaCoord,daCoord,&cArrDa);
368: } else if (dim == 3) {
369: PetscInt ex,ey,ez;
370: PetscScalar ****cArrDa;
371: DMDAVecGetArrayDOF(dmdaCoord,daCoord,&cArrDa);
372: if (daCoordIsStag) {
373: PetscInt slot;
374: PetscScalar ****cArrStag;
375: DMStagGetLocationSlot(dmstagCoord,loc,0,&slot);
376: DMStagVecGetArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
377: for (ez=start[2]; ez<start[2] + n[2] + extraPoint[2]; ++ez) {
378: for (ey=start[1]; ey<start[1] + n[1] + extraPoint[1]; ++ey) {
379: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
380: for (d=0; d<3; ++d) {
381: cArrDa[ez][ey][ex][d] = cArrStag[ez][ey][ex][slot+d];
382: }
383: }
384: }
385: }
386: DMStagVecRestoreArrayDOFRead(dmstagCoord,stagCoord,&cArrStag);
387: } else if (daCoordIsProduct) {
388: PetscScalar **cArrX,**cArrY,**cArrZ;
389: DMStagGet1dCoordinateArraysDOFRead(dmstag,&cArrX,&cArrY,&cArrZ);
390: for (ez=start[2]; ez<start[2] + n[2] + extraPoint[2]; ++ez) {
391: for (ey=start[1]; ey<start[1] + n[1] + extraPoint[1]; ++ey) {
392: for (ex=start[0]; ex<start[0] + n[0] + extraPoint[0]; ++ex) {
393: cArrDa[ez][ey][ex][0] = cArrX[ex][0];
394: cArrDa[ez][ey][ex][1] = cArrY[ey][0];
395: cArrDa[ez][ey][ex][2] = cArrZ[ez][0];
396: }
397: }
398: }
399: DMStagRestore1dCoordinateArraysDOFRead(dmstag,&cArrX,&cArrY,&cArrZ);
400: } else SETERRQ(PetscObjectComm((PetscObject)dmstag),PETSC_ERR_SUP,"Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
401: DMDAVecRestoreArrayDOF(dmdaCoord,daCoord,&cArrDa);
402: } else SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported dimension %d",dim);
403: return(0);
404: }
406: /*
407: Convert to a location value with only BACK, DOWN, LEFT, and ELEMENT involved (makes looping easier)
408: */
409: static PetscErrorCode DMStagStencilLocationCanonicalize(DMStagStencilLocation loc,DMStagStencilLocation *locCanonical)
410: {
412: switch (loc) {
413: case DMSTAG_ELEMENT:
414: *locCanonical = DMSTAG_ELEMENT;
415: break;
416: case DMSTAG_LEFT:
417: case DMSTAG_RIGHT:
418: *locCanonical = DMSTAG_LEFT;
419: break;
420: case DMSTAG_DOWN:
421: case DMSTAG_UP:
422: *locCanonical = DMSTAG_DOWN;
423: break;
424: case DMSTAG_BACK:
425: case DMSTAG_FRONT:
426: *locCanonical = DMSTAG_BACK;
427: break;
428: case DMSTAG_DOWN_LEFT :
429: case DMSTAG_DOWN_RIGHT :
430: case DMSTAG_UP_LEFT :
431: case DMSTAG_UP_RIGHT :
432: *locCanonical = DMSTAG_DOWN_LEFT;
433: break;
434: case DMSTAG_BACK_LEFT:
435: case DMSTAG_BACK_RIGHT:
436: case DMSTAG_FRONT_LEFT:
437: case DMSTAG_FRONT_RIGHT:
438: *locCanonical = DMSTAG_BACK_LEFT;
439: break;
440: case DMSTAG_BACK_DOWN:
441: case DMSTAG_BACK_UP:
442: case DMSTAG_FRONT_DOWN:
443: case DMSTAG_FRONT_UP:
444: *locCanonical = DMSTAG_BACK_DOWN;
445: break;
446: case DMSTAG_BACK_DOWN_LEFT:
447: case DMSTAG_BACK_DOWN_RIGHT:
448: case DMSTAG_BACK_UP_LEFT:
449: case DMSTAG_BACK_UP_RIGHT:
450: case DMSTAG_FRONT_DOWN_LEFT:
451: case DMSTAG_FRONT_DOWN_RIGHT:
452: case DMSTAG_FRONT_UP_LEFT:
453: case DMSTAG_FRONT_UP_RIGHT:
454: *locCanonical = DMSTAG_BACK_DOWN_LEFT;
455: break;
456: default :
457: *locCanonical = DMSTAG_NULL_LOCATION;
458: break;
459: }
460: return(0);
461: }
463: /*@C
464: DMStagVecSplitToDMDA - create a DMDA and Vec from a DMStag and Vec
466: Logically Collective
468: High-level helper function which accepts a DMStag, a global vector, and location/dof,
469: and generates a corresponding DMDA and Vec.
471: Input Parameters:
472: + dm - the DMStag object
473: . vec- Vec object associated with dm
474: . loc - which subgrid to extract (see DMStagStencilLocation)
475: - c - which component to extract (see note below)
477: Output Parameters:
478: + pda - the new DMDA
479: - pdavec - the new Vec
481: Notes:
482: If a c value of -k is provided, the first k dof for that position are extracted,
483: padding with zero values if needbe. If a non-negative value is provided, a single
484: dof is extracted.
486: The caller is responsible for destroying the created DMDA and Vec.
488: Level: advanced
490: .seealso: DMSTAG, DMStagMigrateVec(), DMStagCreateCompatibleDMStag()
491: @*/
492: PetscErrorCode DMStagVecSplitToDMDA(DM dm,Vec vec,DMStagStencilLocation loc,PetscInt c,DM *pda,Vec *pdavec)
493: {
494: PetscErrorCode ierr;
495: PetscInt dim,locdof;
496: DM da,coordDM;
497: Vec davec;
498: DMStagStencilLocation locCanonical;
503: DMGetDimension(dm,&dim);
504: DMStagGetLocationDOF(dm,loc,&locdof);
505: 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);
506: DMStagStencilLocationCanonicalize(loc,&locCanonical);
507: DMStagCreateCompatibleDMDA(dm,locCanonical,c,pda);
508: da = *pda;
509: DMSetUp(*pda);
510: DMGetCoordinateDM(dm,&coordDM);
511: if (coordDM) {
512: DMStagTransferCoordinatesToDMDA(dm,locCanonical,da);
513: }
514: DMCreateGlobalVector(da,pdavec);
515: davec = *pdavec;
516: DMStagMigrateVecDMDA(dm,vec,locCanonical,c,da,davec);
517: return(0);
518: }