Actual source code: dalocal.c
2: /*
3: Code for manipulating distributed regular arrays in parallel.
4: */
6: #include <petsc/private/dmdaimpl.h>
7: #include <petscbt.h>
8: #include <petscsf.h>
9: #include <petscds.h>
10: #include <petscfe.h>
12: /*
13: This allows the DMDA vectors to properly tell MATLAB their dimensions
14: */
15: #if defined(PETSC_HAVE_MATLAB)
16: #include <engine.h> /* MATLAB include file */
17: #include <mex.h> /* MATLAB include file */
18: static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj, void *mengine)
19: {
20: PetscInt n, m;
21: Vec vec = (Vec)obj;
22: PetscScalar *array;
23: mxArray *mat;
24: DM da;
26: PetscFunctionBegin;
27: PetscCall(VecGetDM(vec, &da));
28: PetscCheck(da, PetscObjectComm((PetscObject)vec), PETSC_ERR_ARG_WRONGSTATE, "Vector not associated with a DMDA");
29: PetscCall(DMDAGetGhostCorners(da, 0, 0, 0, &m, &n, 0));
31: PetscCall(VecGetArray(vec, &array));
32: #if !defined(PETSC_USE_COMPLEX)
33: mat = mxCreateDoubleMatrix(m, n, mxREAL);
34: #else
35: mat = mxCreateDoubleMatrix(m, n, mxCOMPLEX);
36: #endif
37: PetscCall(PetscArraycpy(mxGetPr(mat), array, n * m));
38: PetscCall(PetscObjectName(obj));
39: engPutVariable((Engine *)mengine, obj->name, mat);
41: PetscCall(VecRestoreArray(vec, &array));
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
44: #endif
46: PetscErrorCode DMCreateLocalVector_DA(DM da, Vec *g)
47: {
48: DM_DA *dd = (DM_DA *)da->data;
50: PetscFunctionBegin;
53: PetscCall(VecCreate(PETSC_COMM_SELF, g));
54: PetscCall(VecSetSizes(*g, dd->nlocal, PETSC_DETERMINE));
55: PetscCall(VecSetBlockSize(*g, dd->w));
56: PetscCall(VecSetType(*g, da->vectype));
57: if (dd->nlocal < da->bind_below) {
58: PetscCall(VecSetBindingPropagates(*g, PETSC_TRUE));
59: PetscCall(VecBindToCPU(*g, PETSC_TRUE));
60: }
61: PetscCall(VecSetDM(*g, da));
62: #if defined(PETSC_HAVE_MATLAB)
63: if (dd->w == 1 && da->dim == 2) PetscCall(PetscObjectComposeFunction((PetscObject)*g, "PetscMatlabEnginePut_C", VecMatlabEnginePut_DA2d));
64: #endif
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: /*@
69: DMDAGetNumCells - Get the number of cells in the local piece of the `DMDA`. This includes ghost cells.
71: Input Parameter:
72: . dm - The `DMDA` object
74: Output Parameters:
75: + numCellsX - The number of local cells in the x-direction
76: . numCellsY - The number of local cells in the y-direction
77: . numCellsZ - The number of local cells in the z-direction
78: - numCells - The number of local cells
80: Level: developer
82: .seealso: `DM`, `DMDA`, `DMDAGetCellPoint()`
83: @*/
84: PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells)
85: {
86: DM_DA *da = (DM_DA *)dm->data;
87: const PetscInt dim = dm->dim;
88: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
89: const PetscInt nC = (mx) * (dim > 1 ? (my) * (dim > 2 ? (mz) : 1) : 1);
91: PetscFunctionBegin;
93: if (numCellsX) {
95: *numCellsX = mx;
96: }
97: if (numCellsY) {
99: *numCellsY = my;
100: }
101: if (numCellsZ) {
103: *numCellsZ = mz;
104: }
105: if (numCells) {
107: *numCells = nC;
108: }
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: /*@
113: DMDAGetCellPoint - Get the DM point corresponding to the tuple (i, j, k) in the `DMDA`
115: Input Parameters:
116: + dm - The `DMDA` object
117: - i,j,k - The global indices for the cell
119: Output Parameter:
120: . point - The local `DM` point
122: Level: developer
124: .seealso: `DM`, `DMDA`, `DMDAGetNumCells()`
125: @*/
126: PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
127: {
128: const PetscInt dim = dm->dim;
129: DMDALocalInfo info;
131: PetscFunctionBegin;
134: PetscCall(DMDAGetLocalInfo(dm, &info));
135: if (dim > 0) PetscCheck(!(i < info.gxs) && !(i >= info.gxs + info.gxm), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "X index %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", i, info.gxs, info.gxs + info.gxm);
136: if (dim > 1) PetscCheck(!(j < info.gys) && !(j >= info.gys + info.gym), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Y index %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", j, info.gys, info.gys + info.gym);
137: if (dim > 2) PetscCheck(!(k < info.gzs) && !(k >= info.gzs + info.gzm), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Z index %" PetscInt_FMT PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", k, info.gzs, info.gzs + info.gzm);
138: *point = i + (dim > 1 ? (j + (dim > 2 ? k * info.gym : 0)) * info.gxm : 0);
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
143: {
144: DM_DA *da = (DM_DA *)dm->data;
145: const PetscInt dim = dm->dim;
146: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
147: const PetscInt nVx = mx + 1;
148: const PetscInt nVy = dim > 1 ? (my + 1) : 1;
149: const PetscInt nVz = dim > 2 ? (mz + 1) : 1;
150: const PetscInt nV = nVx * nVy * nVz;
152: PetscFunctionBegin;
153: if (numVerticesX) {
155: *numVerticesX = nVx;
156: }
157: if (numVerticesY) {
159: *numVerticesY = nVy;
160: }
161: if (numVerticesZ) {
163: *numVerticesZ = nVz;
164: }
165: if (numVertices) {
167: *numVertices = nV;
168: }
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
173: {
174: DM_DA *da = (DM_DA *)dm->data;
175: const PetscInt dim = dm->dim;
176: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
177: const PetscInt nxF = (dim > 1 ? (my) * (dim > 2 ? (mz) : 1) : 1);
178: const PetscInt nXF = (mx + 1) * nxF;
179: const PetscInt nyF = mx * (dim > 2 ? mz : 1);
180: const PetscInt nYF = dim > 1 ? (my + 1) * nyF : 0;
181: const PetscInt nzF = mx * (dim > 1 ? my : 0);
182: const PetscInt nZF = dim > 2 ? (mz + 1) * nzF : 0;
184: PetscFunctionBegin;
185: if (numXFacesX) {
187: *numXFacesX = nxF;
188: }
189: if (numXFaces) {
191: *numXFaces = nXF;
192: }
193: if (numYFacesY) {
195: *numYFacesY = nyF;
196: }
197: if (numYFaces) {
199: *numYFaces = nYF;
200: }
201: if (numZFacesZ) {
203: *numZFacesZ = nzF;
204: }
205: if (numZFaces) {
207: *numZFaces = nZF;
208: }
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
213: {
214: const PetscInt dim = dm->dim;
215: PetscInt nC, nV, nXF, nYF, nZF;
217: PetscFunctionBegin;
220: PetscCall(DMDAGetNumCells(dm, NULL, NULL, NULL, &nC));
221: PetscCall(DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV));
222: PetscCall(DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF));
223: if (height == 0) {
224: /* Cells */
225: if (pStart) *pStart = 0;
226: if (pEnd) *pEnd = nC;
227: } else if (height == 1) {
228: /* Faces */
229: if (pStart) *pStart = nC + nV;
230: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
231: } else if (height == dim) {
232: /* Vertices */
233: if (pStart) *pStart = nC;
234: if (pEnd) *pEnd = nC + nV;
235: } else if (height < 0) {
236: /* All points */
237: if (pStart) *pStart = 0;
238: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
239: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %" PetscInt_FMT " in the DA", height);
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
244: {
245: const PetscInt dim = dm->dim;
246: PetscInt nC, nV, nXF, nYF, nZF;
248: PetscFunctionBegin;
251: PetscCall(DMDAGetNumCells(dm, NULL, NULL, NULL, &nC));
252: PetscCall(DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV));
253: PetscCall(DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF));
254: if (depth == dim) {
255: /* Cells */
256: if (pStart) *pStart = 0;
257: if (pEnd) *pEnd = nC;
258: } else if (depth == dim - 1) {
259: /* Faces */
260: if (pStart) *pStart = nC + nV;
261: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
262: } else if (depth == 0) {
263: /* Vertices */
264: if (pStart) *pStart = nC;
265: if (pEnd) *pEnd = nC + nV;
266: } else if (depth < 0) {
267: /* All points */
268: if (pStart) *pStart = 0;
269: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
270: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %" PetscInt_FMT " in the DA", depth);
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
275: {
276: const PetscInt dim = dm->dim;
277: PetscInt nC, nV, nXF, nYF, nZF;
279: PetscFunctionBegin;
280: *coneSize = 0;
281: PetscCall(DMDAGetNumCells(dm, NULL, NULL, NULL, &nC));
282: PetscCall(DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV));
283: PetscCall(DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF));
284: switch (dim) {
285: case 2:
286: if (p >= 0) {
287: if (p < nC) {
288: *coneSize = 4;
289: } else if (p < nC + nV) {
290: *coneSize = 0;
291: } else if (p < nC + nV + nXF + nYF + nZF) {
292: *coneSize = 2;
293: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", p, nC + nV + nXF + nYF + nZF);
294: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %" PetscInt_FMT " is invalid", p);
295: break;
296: case 3:
297: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
298: }
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
303: {
304: const PetscInt dim = dm->dim;
305: PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;
307: PetscFunctionBegin;
308: if (!cone) PetscCall(DMGetWorkArray(dm, 6, MPIU_INT, cone));
309: PetscCall(DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC));
310: PetscCall(DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV));
311: PetscCall(DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF));
312: switch (dim) {
313: case 2:
314: if (p >= 0) {
315: if (p < nC) {
316: const PetscInt cy = p / nCx;
317: const PetscInt cx = p % nCx;
319: (*cone)[0] = cy * nxF + cx + nC + nV;
320: (*cone)[1] = cx * nyF + cy + nyF + nC + nV + nXF;
321: (*cone)[2] = cy * nxF + cx + nxF + nC + nV;
322: (*cone)[3] = cx * nyF + cy + nC + nV + nXF;
323: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
324: } else if (p < nC + nV) {
325: } else if (p < nC + nV + nXF) {
326: const PetscInt fy = (p - nC + nV) / nxF;
327: const PetscInt fx = (p - nC + nV) % nxF;
329: (*cone)[0] = fy * nVx + fx + nC;
330: (*cone)[1] = fy * nVx + fx + 1 + nC;
331: } else if (p < nC + nV + nXF + nYF) {
332: const PetscInt fx = (p - nC + nV + nXF) / nyF;
333: const PetscInt fy = (p - nC + nV + nXF) % nyF;
335: (*cone)[0] = fy * nVx + fx + nC;
336: (*cone)[1] = fy * nVx + fx + nVx + nC;
337: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", p, nC + nV + nXF + nYF + nZF);
338: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %" PetscInt_FMT " is invalid", p);
339: break;
340: case 3:
341: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
342: }
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
347: {
348: PetscFunctionBegin;
349: PetscCall(DMGetWorkArray(dm, 6, MPIU_INT, cone));
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
354: {
355: DM_DA *da = (DM_DA *)dm->data;
356: Vec coordinates;
357: PetscSection section;
358: PetscScalar *coords;
359: PetscReal h[3];
360: PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
362: PetscFunctionBegin;
364: PetscCall(DMDAGetInfo(dm, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
365: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "The following code only works for dim <= 3");
366: h[0] = (xu - xl) / M;
367: h[1] = (yu - yl) / N;
368: h[2] = (zu - zl) / P;
369: PetscCall(DMDAGetDepthStratum(dm, 0, &vStart, &vEnd));
370: PetscCall(DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV));
371: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
372: PetscCall(PetscSectionSetNumFields(section, 1));
373: PetscCall(PetscSectionSetFieldComponents(section, 0, dim));
374: PetscCall(PetscSectionSetChart(section, vStart, vEnd));
375: for (v = vStart; v < vEnd; ++v) PetscCall(PetscSectionSetDof(section, v, dim));
376: PetscCall(PetscSectionSetUp(section));
377: PetscCall(PetscSectionGetStorageSize(section, &size));
378: PetscCall(VecCreateSeq(PETSC_COMM_SELF, size, &coordinates));
379: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
380: PetscCall(VecGetArray(coordinates, &coords));
381: for (k = 0; k < nVz; ++k) {
382: PetscInt ind[3], d, off;
384: ind[0] = 0;
385: ind[1] = 0;
386: ind[2] = k + da->zs;
387: for (j = 0; j < nVy; ++j) {
388: ind[1] = j + da->ys;
389: for (i = 0; i < nVx; ++i) {
390: const PetscInt vertex = (k * nVy + j) * nVx + i + vStart;
392: PetscCall(PetscSectionGetOffset(section, vertex, &off));
393: ind[0] = i + da->xs;
394: for (d = 0; d < dim; ++d) coords[off + d] = h[d] * ind[d];
395: }
396: }
397: }
398: PetscCall(VecRestoreArray(coordinates, &coords));
399: PetscCall(DMSetCoordinateSection(dm, PETSC_DETERMINE, section));
400: PetscCall(DMSetCoordinatesLocal(dm, coordinates));
401: PetscCall(PetscSectionDestroy(§ion));
402: PetscCall(VecDestroy(&coordinates));
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: /* ------------------------------------------------------------------- */
408: /*@C
409: DMDAGetArray - Gets a work array for a `DMDA`
411: Input Parameters:
412: + da - information about my local patch
413: - ghosted - do you want arrays for the ghosted or nonghosted patch
415: Output Parameter:
416: . vptr - array data structured
418: Level: advanced
420: Note:
421: The vector values are NOT initialized and may have garbage in them, so you may need
422: to zero them.
424: .seealso: `DM`, `DMDA`, `DMDARestoreArray()`
425: @*/
426: PetscErrorCode DMDAGetArray(DM da, PetscBool ghosted, void *vptr)
427: {
428: PetscInt j, i, xs, ys, xm, ym, zs, zm;
429: char *iarray_start;
430: void **iptr = (void **)vptr;
431: DM_DA *dd = (DM_DA *)da->data;
433: PetscFunctionBegin;
435: if (ghosted) {
436: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
437: if (dd->arrayghostedin[i]) {
438: *iptr = dd->arrayghostedin[i];
439: iarray_start = (char *)dd->startghostedin[i];
440: dd->arrayghostedin[i] = NULL;
441: dd->startghostedin[i] = NULL;
443: goto done;
444: }
445: }
446: xs = dd->Xs;
447: ys = dd->Ys;
448: zs = dd->Zs;
449: xm = dd->Xe - dd->Xs;
450: ym = dd->Ye - dd->Ys;
451: zm = dd->Ze - dd->Zs;
452: } else {
453: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
454: if (dd->arrayin[i]) {
455: *iptr = dd->arrayin[i];
456: iarray_start = (char *)dd->startin[i];
457: dd->arrayin[i] = NULL;
458: dd->startin[i] = NULL;
460: goto done;
461: }
462: }
463: xs = dd->xs;
464: ys = dd->ys;
465: zs = dd->zs;
466: xm = dd->xe - dd->xs;
467: ym = dd->ye - dd->ys;
468: zm = dd->ze - dd->zs;
469: }
471: switch (da->dim) {
472: case 1: {
473: void *ptr;
475: PetscCall(PetscMalloc(xm * sizeof(PetscScalar), &iarray_start));
477: ptr = (void *)(iarray_start - xs * sizeof(PetscScalar));
478: *iptr = (void *)ptr;
479: break;
480: }
481: case 2: {
482: void **ptr;
484: PetscCall(PetscMalloc((ym + 1) * sizeof(void *) + xm * ym * sizeof(PetscScalar), &iarray_start));
486: ptr = (void **)(iarray_start + xm * ym * sizeof(PetscScalar) - ys * sizeof(void *));
487: for (j = ys; j < ys + ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar) * (xm * (j - ys) - xs);
488: *iptr = (void *)ptr;
489: break;
490: }
491: case 3: {
492: void ***ptr, **bptr;
494: PetscCall(PetscMalloc((zm + 1) * sizeof(void **) + (ym * zm + 1) * sizeof(void *) + xm * ym * zm * sizeof(PetscScalar), &iarray_start));
496: ptr = (void ***)(iarray_start + xm * ym * zm * sizeof(PetscScalar) - zs * sizeof(void *));
497: bptr = (void **)(iarray_start + xm * ym * zm * sizeof(PetscScalar) + zm * sizeof(void **));
498: for (i = zs; i < zs + zm; i++) ptr[i] = bptr + ((i - zs) * ym - ys);
499: for (i = zs; i < zs + zm; i++) {
500: for (j = ys; j < ys + ym; j++) ptr[i][j] = iarray_start + sizeof(PetscScalar) * (xm * ym * (i - zs) + xm * (j - ys) - xs);
501: }
502: *iptr = (void *)ptr;
503: break;
504: }
505: default:
506: SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", da->dim);
507: }
509: done:
510: /* add arrays to the checked out list */
511: if (ghosted) {
512: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
513: if (!dd->arrayghostedout[i]) {
514: dd->arrayghostedout[i] = *iptr;
515: dd->startghostedout[i] = iarray_start;
516: break;
517: }
518: }
519: } else {
520: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
521: if (!dd->arrayout[i]) {
522: dd->arrayout[i] = *iptr;
523: dd->startout[i] = iarray_start;
524: break;
525: }
526: }
527: }
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: /*@C
532: DMDARestoreArray - Restores an array of derivative types for a `DMDA`
534: Input Parameters:
535: + da - information about my local patch
536: . ghosted - do you want arrays for the ghosted or nonghosted patch
537: - vptr - array data structured
539: Level: advanced
541: .seealso: `DM`, `DMDA`, `DMDAGetArray()`
542: @*/
543: PetscErrorCode DMDARestoreArray(DM da, PetscBool ghosted, void *vptr)
544: {
545: PetscInt i;
546: void **iptr = (void **)vptr, *iarray_start = NULL;
547: DM_DA *dd = (DM_DA *)da->data;
549: PetscFunctionBegin;
551: if (ghosted) {
552: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
553: if (dd->arrayghostedout[i] == *iptr) {
554: iarray_start = dd->startghostedout[i];
555: dd->arrayghostedout[i] = NULL;
556: dd->startghostedout[i] = NULL;
557: break;
558: }
559: }
560: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
561: if (!dd->arrayghostedin[i]) {
562: dd->arrayghostedin[i] = *iptr;
563: dd->startghostedin[i] = iarray_start;
564: break;
565: }
566: }
567: } else {
568: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
569: if (dd->arrayout[i] == *iptr) {
570: iarray_start = dd->startout[i];
571: dd->arrayout[i] = NULL;
572: dd->startout[i] = NULL;
573: break;
574: }
575: }
576: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
577: if (!dd->arrayin[i]) {
578: dd->arrayin[i] = *iptr;
579: dd->startin[i] = iarray_start;
580: break;
581: }
582: }
583: }
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }