xref: /petsc/src/dm/impls/stag/stagda.c (revision 834855d6effb0d027771461c8e947ee1ce5a1e17) 
1 /* Routines to convert between a (subset of) DMStag and DMDA */
2 
3 #include <petscdmda.h>                /*I  "petscdmda.h"   I*/
4 #include <petsc/private/dmstagimpl.h> /*I  "petscdmstag.h"   I*/
5 #include <petscdmdatypes.h>
6 
DMStagCreateCompatibleDMDA(DM dm,DMStagStencilLocation loc,PetscInt c,DM * dmda)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];
13 
14   PetscFunctionBegin;
15   PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMSTAG);
16   PetscCall(DMGetDimension(dm, &dim));
17 
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   }
24 
25   /* dof */
26   dof = c < 0 ? -c : 1;
27 
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   }
97 
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   }
111 
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 }
129 
130 /*
131 Helper function to get the number of extra points in a DMDA representation for a given canonical location.
132 */
DMStagDMDAGetExtraPoints(DM dm,DMStagStencilLocation locCanonical,PetscInt * extraPoint)133 static PetscErrorCode DMStagDMDAGetExtraPoints(DM dm, DMStagStencilLocation locCanonical, PetscInt *extraPoint)
134 {
135   PetscInt dim, d, nExtra[DMSTAG_MAX_DIM];
136 
137   PetscFunctionBegin;
138   PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMSTAG);
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 }
177 
178 /*
179 Function much like DMStagMigrateVec(), but which accepts an additional position argument to disambiguate which
180 type of DMDA to migrate to.
181 */
182 
DMStagMigrateVecDMDA(DM dm,Vec vec,DMStagStencilLocation loc,PetscInt c,DM dmTo,Vec vecTo)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;
187 
188   PetscFunctionBegin;
189   PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMSTAG);
190   PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
191   PetscValidHeaderSpecificType(dmTo, DM_CLASSID, 5, DMDA);
192   PetscValidHeaderSpecific(vecTo, VEC_CLASSID, 6);
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 }
299 
300 /* Transfer coordinates from a DMStag to a DMDA, specifying which location */
DMStagTransferCoordinatesToDMDA(DM dmstag,DMStagStencilLocation loc,DM dmda)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;
308 
309   PetscFunctionBegin;
310   PetscValidHeaderSpecificType(dmstag, DM_CLASSID, 1, DMSTAG);
311   PetscValidHeaderSpecificType(dmda, DM_CLASSID, 3, DMDA);
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 }
409 
410 /*@
411   DMStagVecSplitToDMDA - create a `DMDA` and `Vec` from a subgrid of a `DMSTAG` and its `Vec`
412 
413   Collective
414 
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)
420 
421   Output Parameters:
422 + pda    - the `DMDA`
423 - pdavec - the new `Vec`
424 
425   Level: advanced
426 
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.
431 
432   The caller is responsible for destroying the created `DMDA` and `Vec`.
433 
434 .seealso: [](ch_stag), `DMSTAG`, `DMDA`, `DMStagStencilLocation`, `DM`, `Vec`, `DMStagMigrateVec()`, `DMStagCreateCompatibleDMStag()`
435 @*/
DMStagVecSplitToDMDA(DM dm,Vec vec,DMStagStencilLocation loc,PetscInt c,DM * pda,Vec * pdavec)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;
442 
443   PetscFunctionBegin;
444   PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMSTAG);
445   PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
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 }
462