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