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: }