Actual source code: stagda.c

  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:   DM_Stag *const  stag = (DM_Stag *)dm->data;
 10:   PetscInt        dim, i, j, stencilWidth, dof, N[DMSTAG_MAX_DIM];
 11:   DMDAStencilType stencilType;
 12:   PetscInt       *l[DMSTAG_MAX_DIM];

 14:   PetscFunctionBegin;
 16:   PetscCall(DMGetDimension(dm, &dim));

 18:   /* Create grid decomposition (to be adjusted later) */
 19:   for (i = 0; i < dim; ++i) {
 20:     PetscCall(PetscMalloc1(stag->nRanks[i], &l[i]));
 21:     for (j = 0; j < stag->nRanks[i]; ++j) l[i][j] = stag->l[i][j];
 22:     N[i] = stag->N[i];
 23:   }

 25:   /* dof */
 26:   dof = c < 0 ? -c : 1;

 28:   /* Determine/adjust sizes */
 29:   switch (loc) {
 30:   case DMSTAG_ELEMENT:
 31:     break;
 32:   case DMSTAG_LEFT:
 33:   case DMSTAG_RIGHT:
 34:     PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
 35:     l[0][stag->nRanks[0] - 1] += 1; /* extra vertex in direction 0 on last rank in dimension 0 */
 36:     N[0] += 1;
 37:     break;
 38:   case DMSTAG_UP:
 39:   case DMSTAG_DOWN:
 40:     PetscCheck(dim >= 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
 41:     l[1][stag->nRanks[1] - 1] += 1; /* extra vertex in direction 1 on last rank in dimension 1 */
 42:     N[1] += 1;
 43:     break;
 44:   case DMSTAG_BACK:
 45:   case DMSTAG_FRONT:
 46:     PetscCheck(dim >= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
 47:     l[2][stag->nRanks[2] - 1] += 1; /* extra vertex in direction 2 on last rank in dimension 2 */
 48:     N[2] += 1;
 49:     break;
 50:   case DMSTAG_DOWN_LEFT:
 51:   case DMSTAG_DOWN_RIGHT:
 52:   case DMSTAG_UP_LEFT:
 53:   case DMSTAG_UP_RIGHT:
 54:     PetscCheck(dim >= 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
 55:     for (i = 0; i < 2; ++i) { /* extra vertex in direction i on last rank in dimension i = 0,1 */
 56:       l[i][stag->nRanks[i] - 1] += 1;
 57:       N[i] += 1;
 58:     }
 59:     break;
 60:   case DMSTAG_BACK_LEFT:
 61:   case DMSTAG_BACK_RIGHT:
 62:   case DMSTAG_FRONT_LEFT:
 63:   case DMSTAG_FRONT_RIGHT:
 64:     PetscCheck(dim >= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
 65:     for (i = 0; i < 3; i += 2) { /* extra vertex in direction i on last rank in dimension i = 0,2 */
 66:       l[i][stag->nRanks[i] - 1] += 1;
 67:       N[i] += 1;
 68:     }
 69:     break;
 70:   case DMSTAG_BACK_DOWN:
 71:   case DMSTAG_BACK_UP:
 72:   case DMSTAG_FRONT_DOWN:
 73:   case DMSTAG_FRONT_UP:
 74:     PetscCheck(dim >= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
 75:     for (i = 1; i < 3; ++i) { /* extra vertex in direction i on last rank in dimension i = 1,2 */
 76:       l[i][stag->nRanks[i] - 1] += 1;
 77:       N[i] += 1;
 78:     }
 79:     break;
 80:   case DMSTAG_BACK_DOWN_LEFT:
 81:   case DMSTAG_BACK_DOWN_RIGHT:
 82:   case DMSTAG_BACK_UP_LEFT:
 83:   case DMSTAG_BACK_UP_RIGHT:
 84:   case DMSTAG_FRONT_DOWN_LEFT:
 85:   case DMSTAG_FRONT_DOWN_RIGHT:
 86:   case DMSTAG_FRONT_UP_LEFT:
 87:   case DMSTAG_FRONT_UP_RIGHT:
 88:     PetscCheck(dim >= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Incompatible dim (%" PetscInt_FMT ") and loc(%s) combination", dim, DMStagStencilLocations[loc]);
 89:     for (i = 0; i < 3; ++i) { /* extra vertex in direction i on last rank in dimension i = 0,1,2 */
 90:       l[i][stag->nRanks[i] - 1] += 1;
 91:       N[i] += 1;
 92:     }
 93:     break;
 94:   default:
 95:     SETERRQ(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:
101:     stencilType  = DMDA_STENCIL_STAR;
102:     stencilWidth = stag->stencilWidth;
103:     break;
104:   case DMSTAG_STENCIL_BOX:
105:     stencilType  = DMDA_STENCIL_BOX;
106:     stencilWidth = stag->stencilWidth;
107:     break;
108:   default:
109:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported Stencil Type %d", stag->stencilType);
110:   }

112:   /* Create DMDA, using same boundary type */
113:   switch (dim) {
114:   case 1:
115:     PetscCall(DMDACreate1d(PetscObjectComm((PetscObject)dm), stag->boundaryType[0], N[0], dof, stencilWidth, l[0], dmda));
116:     break;
117:   case 2:
118:     PetscCall(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));
119:     break;
120:   case 3:
121:     PetscCall(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));
122:     break;
123:   default:
124:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "not implemented for dim %" PetscInt_FMT, dim);
125:   }
126:   for (i = 0; i < dim; ++i) PetscCall(PetscFree(l[i]));
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

130: /*
131: Helper function to get the number of extra points in a DMDA representation for a given canonical location.
132: */
133: static PetscErrorCode DMStagDMDAGetExtraPoints(DM dm, DMStagStencilLocation locCanonical, PetscInt *extraPoint)
134: {
135:   PetscInt dim, d, nExtra[DMSTAG_MAX_DIM];

137:   PetscFunctionBegin;
139:   PetscCall(DMGetDimension(dm, &dim));
140:   PetscCheck(dim <= DMSTAG_MAX_DIM, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " dimensions", dim);
141:   PetscCall(DMStagGetCorners(dm, NULL, NULL, NULL, NULL, NULL, NULL, &nExtra[0], &nExtra[1], &nExtra[2]));
142:   for (d = 0; d < dim; ++d) extraPoint[d] = 0;
143:   switch (locCanonical) {
144:   case DMSTAG_ELEMENT:
145:     break; /* no extra points */
146:   case DMSTAG_LEFT:
147:     extraPoint[0] = nExtra[0];
148:     break; /* only extra point in x */
149:   case DMSTAG_DOWN:
150:     extraPoint[1] = nExtra[1];
151:     break; /* only extra point in y */
152:   case DMSTAG_BACK:
153:     extraPoint[2] = nExtra[2];
154:     break; /* only extra point in z */
155:   case DMSTAG_DOWN_LEFT:
156:     extraPoint[0] = nExtra[0];
157:     extraPoint[1] = nExtra[1];
158:     break; /* extra point in both x and y  */
159:   case DMSTAG_BACK_LEFT:
160:     extraPoint[0] = nExtra[0];
161:     extraPoint[2] = nExtra[2];
162:     break; /* extra point in both x and z  */
163:   case DMSTAG_BACK_DOWN:
164:     extraPoint[1] = nExtra[1];
165:     extraPoint[2] = nExtra[2];
166:     break; /* extra point in both y and z  */
167:   case DMSTAG_BACK_DOWN_LEFT:
168:     extraPoint[0] = nExtra[0];
169:     extraPoint[1] = nExtra[1];
170:     extraPoint[2] = nExtra[2];
171:     break; /* extra points in x,y,z */
172:   default:
173:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location (perhaps not canonical) %s", DMStagStencilLocations[locCanonical]);
174:   }
175:   PetscFunctionReturn(PETSC_SUCCESS);
176: }

178: /*
179: Function much like DMStagMigrateVec(), but which accepts an additional position argument to disambiguate which
180: type of DMDA to migrate to.
181: */

183: static PetscErrorCode DMStagMigrateVecDMDA(DM dm, Vec vec, DMStagStencilLocation loc, PetscInt c, DM dmTo, Vec vecTo)
184: {
185:   PetscInt i, j, k, d, dim, dof, dofToMax, start[DMSTAG_MAX_DIM], n[DMSTAG_MAX_DIM], extraPoint[DMSTAG_MAX_DIM];
186:   Vec      vecLocal;

188:   PetscFunctionBegin;
193:   PetscCall(DMGetDimension(dm, &dim));
194:   PetscCall(DMDAGetDof(dmTo, &dofToMax));
195:   PetscCheck(-c <= dofToMax, PetscObjectComm((PetscObject)dmTo), PETSC_ERR_ARG_OUTOFRANGE, "Invalid negative component value. Must be >= -%" PetscInt_FMT, dofToMax);
196:   PetscCall(DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], NULL, NULL, NULL));
197:   PetscCall(DMStagDMDAGetExtraPoints(dm, loc, extraPoint));
198:   PetscCall(DMStagGetLocationDOF(dm, loc, &dof));
199:   PetscCall(DMGetLocalVector(dm, &vecLocal));
200:   PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal));
201:   PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal));
202:   if (dim == 1) {
203:     PetscScalar **arrTo;
204:     PetscCall(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;
211:           pos.loc = loc;
212:           pos.c   = d;
213:           PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[i][d]));
214:         }
215:         for (; d < dofTo; ++d) { arrTo[i][d] = 0.0; /* Pad extra dof with zeros */ }
216:       }
217:     } else {
218:       for (i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
219:         DMStagStencil pos;
220:         pos.i   = i;
221:         pos.loc = loc;
222:         pos.c   = c;
223:         PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[i][0]));
224:       }
225:     }
226:     PetscCall(DMDAVecRestoreArrayDOF(dmTo, vecTo, &arrTo));
227:   } else if (dim == 2) {
228:     PetscScalar ***arrTo;
229:     PetscCall(DMDAVecGetArrayDOF(dmTo, vecTo, &arrTo));
230:     if (c < 0) {
231:       const PetscInt dofTo = -c;
232:       for (j = start[1]; j < start[1] + n[1] + extraPoint[1]; ++j) {
233:         for (i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
234:           for (d = 0; d < PetscMin(dof, dofTo); ++d) {
235:             DMStagStencil pos;
236:             pos.i   = i;
237:             pos.j   = j;
238:             pos.loc = loc;
239:             pos.c   = d;
240:             PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[j][i][d]));
241:           }
242:           for (; d < dofTo; ++d) { arrTo[j][i][d] = 0.0; /* Pad extra dof with zeros */ }
243:         }
244:       }
245:     } else {
246:       for (j = start[1]; j < start[1] + n[1] + extraPoint[1]; ++j) {
247:         for (i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
248:           DMStagStencil pos;
249:           pos.i   = i;
250:           pos.j   = j;
251:           pos.loc = loc;
252:           pos.c   = c;
253:           PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[j][i][0]));
254:         }
255:       }
256:     }
257:     PetscCall(DMDAVecRestoreArrayDOF(dmTo, vecTo, &arrTo));
258:   } else if (dim == 3) {
259:     PetscScalar ****arrTo;
260:     PetscCall(DMDAVecGetArrayDOF(dmTo, vecTo, &arrTo));
261:     if (c < 0) {
262:       const PetscInt dofTo = -c;
263:       for (k = start[2]; k < start[2] + n[2] + extraPoint[2]; ++k) {
264:         for (j = start[1]; j < start[1] + n[1] + extraPoint[1]; ++j) {
265:           for (i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
266:             for (d = 0; d < PetscMin(dof, dofTo); ++d) {
267:               DMStagStencil pos;
268:               pos.i   = i;
269:               pos.j   = j;
270:               pos.k   = k;
271:               pos.loc = loc;
272:               pos.c   = d;
273:               PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[k][j][i][d]));
274:             }
275:             for (; d < dofTo; ++d) { arrTo[k][j][i][d] = 0.0; /* Pad extra dof with zeros */ }
276:           }
277:         }
278:       }
279:     } else {
280:       for (k = start[2]; k < start[2] + n[2] + extraPoint[2]; ++k) {
281:         for (j = start[1]; j < start[1] + n[1] + extraPoint[1]; ++j) {
282:           for (i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
283:             DMStagStencil pos;
284:             pos.i   = i;
285:             pos.j   = j;
286:             pos.k   = k;
287:             pos.loc = loc;
288:             pos.c   = c;
289:             PetscCall(DMStagVecGetValuesStencil(dm, vecLocal, 1, &pos, &arrTo[k][j][i][0]));
290:           }
291:         }
292:       }
293:     }
294:     PetscCall(DMDAVecRestoreArrayDOF(dmTo, vecTo, &arrTo));
295:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
296:   PetscCall(DMRestoreLocalVector(dm, &vecLocal));
297:   PetscFunctionReturn(PETSC_SUCCESS);
298: }

300: /* Transfer coordinates from a DMStag to a DMDA, specifying which location */
301: static PetscErrorCode DMStagTransferCoordinatesToDMDA(DM dmstag, DMStagStencilLocation loc, DM dmda)
302: {
303:   PetscInt  dim, start[DMSTAG_MAX_DIM], n[DMSTAG_MAX_DIM], extraPoint[DMSTAG_MAX_DIM], d;
304:   DM        dmstagCoord, dmdaCoord;
305:   DMType    dmstagCoordType;
306:   Vec       stagCoord, daCoord;
307:   PetscBool daCoordIsStag, daCoordIsProduct;

309:   PetscFunctionBegin;
312:   PetscCall(DMGetDimension(dmstag, &dim));
313:   PetscCall(DMGetCoordinateDM(dmstag, &dmstagCoord));
314:   PetscCall(DMGetCoordinatesLocal(dmstag, &stagCoord)); /* Note local */
315:   PetscCall(DMGetCoordinateDM(dmda, &dmdaCoord));
316:   daCoord = NULL;
317:   PetscCall(DMGetCoordinates(dmda, &daCoord));
318:   if (!daCoord) {
319:     PetscCall(DMCreateGlobalVector(dmdaCoord, &daCoord));
320:     PetscCall(DMSetCoordinates(dmda, daCoord));
321:     PetscCall(VecDestroy(&daCoord));
322:     PetscCall(DMGetCoordinates(dmda, &daCoord));
323:   }
324:   PetscCall(DMGetType(dmstagCoord, &dmstagCoordType));
325:   PetscCall(PetscStrcmp(dmstagCoordType, DMSTAG, &daCoordIsStag));
326:   PetscCall(PetscStrcmp(dmstagCoordType, DMPRODUCT, &daCoordIsProduct));
327:   PetscCall(DMStagGetCorners(dmstag, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], NULL, NULL, NULL));
328:   PetscCall(DMStagDMDAGetExtraPoints(dmstag, loc, extraPoint));
329:   if (dim == 1) {
330:     PetscInt      ex;
331:     PetscScalar **cArrDa;
332:     PetscCall(DMDAVecGetArrayDOF(dmdaCoord, daCoord, &cArrDa));
333:     if (daCoordIsStag) {
334:       PetscInt      slot;
335:       PetscScalar **cArrStag;
336:       PetscCall(DMStagGetLocationSlot(dmstagCoord, loc, 0, &slot));
337:       PetscCall(DMStagVecGetArrayRead(dmstagCoord, stagCoord, &cArrStag));
338:       for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) cArrDa[ex][0] = cArrStag[ex][slot];
339:       PetscCall(DMStagVecRestoreArrayRead(dmstagCoord, stagCoord, &cArrStag));
340:     } else if (daCoordIsProduct) {
341:       PetscScalar **cArrX;
342:       PetscCall(DMStagGetProductCoordinateArraysRead(dmstag, &cArrX, NULL, NULL));
343:       for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) cArrDa[ex][0] = cArrX[ex][0];
344:       PetscCall(DMStagRestoreProductCoordinateArraysRead(dmstag, &cArrX, NULL, NULL));
345:     } else SETERRQ(PetscObjectComm((PetscObject)dmstag), PETSC_ERR_SUP, "Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
346:     PetscCall(DMDAVecRestoreArrayDOF(dmdaCoord, daCoord, &cArrDa));
347:   } else if (dim == 2) {
348:     PetscInt       ex, ey;
349:     PetscScalar ***cArrDa;
350:     PetscCall(DMDAVecGetArrayDOF(dmdaCoord, daCoord, &cArrDa));
351:     if (daCoordIsStag) {
352:       PetscInt       slot;
353:       PetscScalar ***cArrStag;
354:       PetscCall(DMStagGetLocationSlot(dmstagCoord, loc, 0, &slot));
355:       PetscCall(DMStagVecGetArrayRead(dmstagCoord, stagCoord, &cArrStag));
356:       for (ey = start[1]; ey < start[1] + n[1] + extraPoint[1]; ++ey) {
357:         for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) {
358:           for (d = 0; d < 2; ++d) cArrDa[ey][ex][d] = cArrStag[ey][ex][slot + d];
359:         }
360:       }
361:       PetscCall(DMStagVecRestoreArrayRead(dmstagCoord, stagCoord, &cArrStag));
362:     } else if (daCoordIsProduct) {
363:       PetscScalar **cArrX, **cArrY;
364:       PetscCall(DMStagGetProductCoordinateArraysRead(dmstag, &cArrX, &cArrY, NULL));
365:       for (ey = start[1]; ey < start[1] + n[1] + extraPoint[1]; ++ey) {
366:         for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) {
367:           cArrDa[ey][ex][0] = cArrX[ex][0];
368:           cArrDa[ey][ex][1] = cArrY[ey][0];
369:         }
370:       }
371:       PetscCall(DMStagRestoreProductCoordinateArraysRead(dmstag, &cArrX, &cArrY, NULL));
372:     } else SETERRQ(PetscObjectComm((PetscObject)dmstag), PETSC_ERR_SUP, "Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
373:     PetscCall(DMDAVecRestoreArrayDOF(dmdaCoord, daCoord, &cArrDa));
374:   } else if (dim == 3) {
375:     PetscInt        ex, ey, ez;
376:     PetscScalar ****cArrDa;
377:     PetscCall(DMDAVecGetArrayDOF(dmdaCoord, daCoord, &cArrDa));
378:     if (daCoordIsStag) {
379:       PetscInt        slot;
380:       PetscScalar ****cArrStag;
381:       PetscCall(DMStagGetLocationSlot(dmstagCoord, loc, 0, &slot));
382:       PetscCall(DMStagVecGetArrayRead(dmstagCoord, stagCoord, &cArrStag));
383:       for (ez = start[2]; ez < start[2] + n[2] + extraPoint[2]; ++ez) {
384:         for (ey = start[1]; ey < start[1] + n[1] + extraPoint[1]; ++ey) {
385:           for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) {
386:             for (d = 0; d < 3; ++d) cArrDa[ez][ey][ex][d] = cArrStag[ez][ey][ex][slot + d];
387:           }
388:         }
389:       }
390:       PetscCall(DMStagVecRestoreArrayRead(dmstagCoord, stagCoord, &cArrStag));
391:     } else if (daCoordIsProduct) {
392:       PetscScalar **cArrX, **cArrY, **cArrZ;
393:       PetscCall(DMStagGetProductCoordinateArraysRead(dmstag, &cArrX, &cArrY, &cArrZ));
394:       for (ez = start[2]; ez < start[2] + n[2] + extraPoint[2]; ++ez) {
395:         for (ey = start[1]; ey < start[1] + n[1] + extraPoint[1]; ++ey) {
396:           for (ex = start[0]; ex < start[0] + n[0] + extraPoint[0]; ++ex) {
397:             cArrDa[ez][ey][ex][0] = cArrX[ex][0];
398:             cArrDa[ez][ey][ex][1] = cArrY[ey][0];
399:             cArrDa[ez][ey][ex][2] = cArrZ[ez][0];
400:           }
401:         }
402:       }
403:       PetscCall(DMStagRestoreProductCoordinateArraysRead(dmstag, &cArrX, &cArrY, &cArrZ));
404:     } else SETERRQ(PetscObjectComm((PetscObject)dmstag), PETSC_ERR_SUP, "Stag to DA coordinate transfer only supported for DMStag coordinate DM of type DMstag or DMProduct");
405:     PetscCall(DMDAVecRestoreArrayDOF(dmdaCoord, daCoord, &cArrDa));
406:   } else SETERRQ(PetscObjectComm((PetscObject)dmstag), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
407:   PetscFunctionReturn(PETSC_SUCCESS);
408: }

410: /*@
411:   DMStagVecSplitToDMDA - create a `DMDA` and `Vec` from a subgrid of a `DMSTAG` and its `Vec`

413:   Collective

415:   Input Parameters:
416: + dm  - the `DMSTAG` object
417: . vec - `Vec` object associated with `dm`
418: . loc - which subgrid to extract (see `DMStagStencilLocation`)
419: - c   - which component to extract (see note below)

421:   Output Parameters:
422: + pda    - the `DMDA`
423: - pdavec - the new `Vec`

425:   Level: advanced

427:   Notes:
428:   If a `c` value of `-k` is provided, the first `k` DOF for that position are extracted,
429:   padding with zero values if needed. If a non-negative value is provided, a single
430:   DOF is extracted.

432:   The caller is responsible for destroying the created `DMDA` and `Vec`.

434: .seealso: [](ch_stag), `DMSTAG`, `DMDA`, `DMStagStencilLocation`, `DM`, `Vec`, `DMStagMigrateVec()`, `DMStagCreateCompatibleDMStag()`
435: @*/
436: PetscErrorCode DMStagVecSplitToDMDA(DM dm, Vec vec, DMStagStencilLocation loc, PetscInt c, DM *pda, Vec *pdavec)
437: {
438:   PetscInt              dim, locdof;
439:   DM                    da, coordDM;
440:   Vec                   davec;
441:   DMStagStencilLocation locCanonical;

443:   PetscFunctionBegin;
446:   PetscCall(DMGetDimension(dm, &dim));
447:   PetscCall(DMStagGetLocationDOF(dm, loc, &locdof));
448:   PetscCheck(c < locdof, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Location %s has %" PetscInt_FMT " dof, but component %" PetscInt_FMT " requested", DMStagStencilLocations[loc], locdof, c);
449:   PetscCall(DMStagStencilLocationCanonicalize(loc, &locCanonical));
450:   PetscCall(DMStagCreateCompatibleDMDA(dm, locCanonical, c, pda));
451:   da = *pda;
452:   PetscCall(DMSetUp(*pda));
453:   if (dm->coordinates[0].dm != NULL) {
454:     PetscCall(DMGetCoordinateDM(dm, &coordDM));
455:     PetscCall(DMStagTransferCoordinatesToDMDA(dm, locCanonical, da));
456:   }
457:   PetscCall(DMCreateGlobalVector(da, pdavec));
458:   davec = *pdavec;
459:   PetscCall(DMStagMigrateVecDMDA(dm, vec, locCanonical, c, da, davec));
460:   PetscFunctionReturn(PETSC_SUCCESS);
461: }