Actual source code: dalocal.c
1: /*
2: Code for manipulating distributed regular arrays in parallel.
3: */
5: #include <petsc/private/dmdaimpl.h>
6: #include <petscbt.h>
7: #include <petscsf.h>
8: #include <petscds.h>
9: #include <petscfe.h>
11: /*
12: This allows the DMDA vectors to properly tell MATLAB their dimensions
13: */
14: #if defined(PETSC_HAVE_MATLAB)
15: #include <engine.h> /* MATLAB include file */
16: #include <mex.h> /* MATLAB include file */
17: static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj, void *mengine)
18: {
19: PetscInt n, m;
20: Vec vec = (Vec)obj;
21: PetscScalar *array;
22: mxArray *mat;
23: DM da;
25: PetscFunctionBegin;
26: PetscCall(VecGetDM(vec, &da));
27: PetscCheck(da, PetscObjectComm((PetscObject)vec), PETSC_ERR_ARG_WRONGSTATE, "Vector not associated with a DMDA");
28: PetscCall(DMDAGetGhostCorners(da, 0, 0, 0, &m, &n, 0));
30: PetscCall(VecGetArray(vec, &array));
31: #if !defined(PETSC_USE_COMPLEX)
32: mat = mxCreateDoubleMatrix(m, n, mxREAL);
33: #else
34: mat = mxCreateDoubleMatrix(m, n, mxCOMPLEX);
35: #endif
36: PetscCall(PetscArraycpy(mxGetPr(mat), array, n * m));
37: PetscCall(PetscObjectName(obj));
38: engPutVariable((Engine *)mengine, obj->name, mat);
40: PetscCall(VecRestoreArray(vec, &array));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
43: #endif
45: PetscErrorCode DMCreateLocalVector_DA(DM da, Vec *g)
46: {
47: DM_DA *dd = (DM_DA *)da->data;
49: PetscFunctionBegin;
51: PetscAssertPointer(g, 2);
52: PetscCall(VecCreate(PETSC_COMM_SELF, g));
53: PetscCall(VecSetSizes(*g, dd->nlocal, PETSC_DETERMINE));
54: PetscCall(VecSetBlockSize(*g, dd->w));
55: PetscCall(VecSetType(*g, da->vectype));
56: if (dd->nlocal < da->bind_below) {
57: PetscCall(VecSetBindingPropagates(*g, PETSC_TRUE));
58: PetscCall(VecBindToCPU(*g, PETSC_TRUE));
59: }
60: PetscCall(VecSetDM(*g, da));
61: #if defined(PETSC_HAVE_MATLAB)
62: if (dd->w == 1 && da->dim == 2) PetscCall(PetscObjectComposeFunction((PetscObject)*g, "PetscMatlabEnginePut_C", VecMatlabEnginePut_DA2d));
63: #endif
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: /*@
68: DMDAGetNumCells - Get the number of cells (or vertices) in the local piece of the `DMDA`. This includes ghost cells.
70: Input Parameter:
71: . dm - The `DMDA` object
73: Output Parameters:
74: + numCellsX - The number of local cells in the x-direction
75: . numCellsY - The number of local cells in the y-direction
76: . numCellsZ - The number of local cells in the z-direction
77: - numCells - The number of local cells
79: Level: developer
81: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetCellPoint()`
82: @*/
83: PetscErrorCode DMDAGetNumCells(DM dm, PeOp PetscInt *numCellsX, PeOp PetscInt *numCellsY, PeOp PetscInt *numCellsZ, PeOp PetscInt *numCells)
84: {
85: DM_DA *da = (DM_DA *)dm->data;
86: const PetscInt dim = dm->dim;
87: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
88: const PetscInt nC = (mx) * (dim > 1 ? (my) * (dim > 2 ? (mz) : 1) : 1);
90: PetscFunctionBegin;
92: if (numCellsX) {
93: PetscAssertPointer(numCellsX, 2);
94: *numCellsX = mx;
95: }
96: if (numCellsY) {
97: PetscAssertPointer(numCellsY, 3);
98: *numCellsY = my;
99: }
100: if (numCellsZ) {
101: PetscAssertPointer(numCellsZ, 4);
102: *numCellsZ = mz;
103: }
104: if (numCells) {
105: PetscAssertPointer(numCells, 5);
106: *numCells = nC;
107: }
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: /*@
112: DMDAGetCellPoint - Get the `DM` point corresponding to the tuple (i, j, k) in the `DMDA`
114: Input Parameters:
115: + dm - The `DMDA` object
116: . i - The global x index for the cell
117: . j - The global y index for the cell
118: - k - The global z index for the cell
120: Output Parameter:
121: . point - The local `DM` point
123: Level: developer
125: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetNumCells()`
126: @*/
127: PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
128: {
129: const PetscInt dim = dm->dim;
130: DMDALocalInfo info;
132: PetscFunctionBegin;
134: PetscAssertPointer(point, 5);
135: PetscCall(DMDAGetLocalInfo(dm, &info));
136: 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);
137: 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);
138: if (dim > 2) PetscCheck(!(k < info.gzs) && !(k >= info.gzs + info.gzm), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Z index %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", k, info.gzs, info.gzs + info.gzm);
139: *point = i + (dim > 1 ? (j + (dim > 2 ? k * info.gym : 0)) * info.gxm : 0);
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
144: {
145: DM_DA *da = (DM_DA *)dm->data;
146: const PetscInt dim = dm->dim;
147: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
148: const PetscInt nVx = mx + 1;
149: const PetscInt nVy = dim > 1 ? (my + 1) : 1;
150: const PetscInt nVz = dim > 2 ? (mz + 1) : 1;
151: const PetscInt nV = nVx * nVy * nVz;
153: PetscFunctionBegin;
154: if (numVerticesX) {
155: PetscAssertPointer(numVerticesX, 2);
156: *numVerticesX = nVx;
157: }
158: if (numVerticesY) {
159: PetscAssertPointer(numVerticesY, 3);
160: *numVerticesY = nVy;
161: }
162: if (numVerticesZ) {
163: PetscAssertPointer(numVerticesZ, 4);
164: *numVerticesZ = nVz;
165: }
166: if (numVertices) {
167: PetscAssertPointer(numVertices, 5);
168: *numVertices = nV;
169: }
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
174: {
175: DM_DA *da = (DM_DA *)dm->data;
176: const PetscInt dim = dm->dim;
177: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
178: const PetscInt nxF = (dim > 1 ? (my) * (dim > 2 ? (mz) : 1) : 1);
179: const PetscInt nXF = (mx + 1) * nxF;
180: const PetscInt nyF = mx * (dim > 2 ? mz : 1);
181: const PetscInt nYF = dim > 1 ? (my + 1) * nyF : 0;
182: const PetscInt nzF = mx * (dim > 1 ? my : 0);
183: const PetscInt nZF = dim > 2 ? (mz + 1) * nzF : 0;
185: PetscFunctionBegin;
186: if (numXFacesX) {
187: PetscAssertPointer(numXFacesX, 2);
188: *numXFacesX = nxF;
189: }
190: if (numXFaces) {
191: PetscAssertPointer(numXFaces, 3);
192: *numXFaces = nXF;
193: }
194: if (numYFacesY) {
195: PetscAssertPointer(numYFacesY, 4);
196: *numYFacesY = nyF;
197: }
198: if (numYFaces) {
199: PetscAssertPointer(numYFaces, 5);
200: *numYFaces = nYF;
201: }
202: if (numZFacesZ) {
203: PetscAssertPointer(numZFacesZ, 6);
204: *numZFacesZ = nzF;
205: }
206: if (numZFaces) {
207: PetscAssertPointer(numZFaces, 7);
208: *numZFaces = nZF;
209: }
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*@
214: DMDAGetHeightStratum - Get the bounds [`start`, `end`) for all points at a certain height.
216: Not Collective
218: Input Parameters:
219: + dm - The `DMDA` object
220: - height - The requested height
222: Output Parameters:
223: + pStart - The first point at this `height`
224: - pEnd - One beyond the last point at this `height`
226: Level: developer
228: Note:
229: See `DMPlexGetHeightStratum()` for the meaning of these values
231: .seealso: [](ch_unstructured), `DM`, `DMDA`, `DMPlexGetDepthStratum()`, `DMPlexGetHeightStratum()`, `DMPlexGetCellTypeStratum()`, `DMPlexGetDepth()`,
232: `DMPlexGetDepthLabel()`, `DMPlexGetPointDepth()`, `DMPlexSymmetrize()`, `DMPlexInterpolate()`, `DMDAGetDepthStratum()`
233: @*/
234: PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PeOp PetscInt *pStart, PeOp PetscInt *pEnd)
235: {
236: const PetscInt dim = dm->dim;
237: PetscInt nC, nV, nXF, nYF, nZF;
239: PetscFunctionBegin;
240: if (pStart) PetscAssertPointer(pStart, 3);
241: if (pEnd) PetscAssertPointer(pEnd, 4);
242: PetscCall(DMDAGetNumCells(dm, NULL, NULL, NULL, &nC));
243: PetscCall(DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV));
244: PetscCall(DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF));
245: if (height == 0) {
246: /* Cells */
247: if (pStart) *pStart = 0;
248: if (pEnd) *pEnd = nC;
249: } else if (height == 1) {
250: /* Faces */
251: if (pStart) *pStart = nC + nV;
252: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
253: } else if (height == dim) {
254: /* Vertices */
255: if (pStart) *pStart = nC;
256: if (pEnd) *pEnd = nC + nV;
257: } else if (height < 0) {
258: /* All points */
259: if (pStart) *pStart = 0;
260: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
261: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %" PetscInt_FMT " in the DA", height);
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /*@
266: DMDAGetDepthStratum - Get the bounds [`start`, `end`) for all points at a certain depth.
268: Not Collective
270: Input Parameters:
271: + dm - The `DMDA` object
272: - depth - The requested depth
274: Output Parameters:
275: + pStart - The first point at this `depth`
276: - pEnd - One beyond the last point at this `depth`
278: Level: developer
280: Note:
281: See `DMPlexGetDepthStratum()` for the meaning of these values
283: .seealso: [](ch_unstructured), `DM`, `DMDA`, `DMPlexGetDepthStratum()`, `DMPlexGetHeightStratum()`, `DMPlexGetCellTypeStratum()`, `DMPlexGetDepth()`,
284: `DMPlexGetDepthLabel()`, `DMPlexGetPointDepth()`, `DMPlexSymmetrize()`, `DMPlexInterpolate()`, `DMDAGetHeightStratum()`
285: @*/
286: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PeOp PetscInt *pStart, PeOp PetscInt *pEnd)
287: {
288: const PetscInt dim = dm->dim;
289: PetscInt nC, nV, nXF, nYF, nZF;
291: PetscFunctionBegin;
292: if (pStart) PetscAssertPointer(pStart, 3);
293: if (pEnd) PetscAssertPointer(pEnd, 4);
294: PetscCall(DMDAGetNumCells(dm, NULL, NULL, NULL, &nC));
295: PetscCall(DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV));
296: PetscCall(DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF));
297: if (depth == dim) {
298: /* Cells */
299: if (pStart) *pStart = 0;
300: if (pEnd) *pEnd = nC;
301: } else if (depth == dim - 1) {
302: /* Faces */
303: if (pStart) *pStart = nC + nV;
304: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
305: } else if (depth == 0) {
306: /* Vertices */
307: if (pStart) *pStart = nC;
308: if (pEnd) *pEnd = nC + nV;
309: } else if (depth < 0) {
310: /* All points */
311: if (pStart) *pStart = 0;
312: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
313: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %" PetscInt_FMT " in the DA", depth);
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: /*@
318: DMDASetVertexCoordinates - Sets the lower and upper coordinates for a `DMDA`
320: Logically Collective
322: Input Parameters:
323: + dm - The `DMDA` object
324: . xl - the lower x coordinate
325: . xu - the upper x coordinate
326: . yl - the lower y coordinate
327: . yu - the upper y coordinate
328: . zl - the lower z coordinate
329: - zu - the upper z coordinate
331: Level: intermediate
333: .seealso: [](ch_unstructured), `DM`, `DMDA`
334: @*/
335: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
336: {
337: DM_DA *da = (DM_DA *)dm->data;
338: Vec coordinates;
339: PetscSection section;
340: PetscScalar *coords;
341: PetscReal h[3];
342: PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
344: PetscFunctionBegin;
346: PetscCall(DMDAGetInfo(dm, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
347: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "The following code only works for dim <= 3");
348: h[0] = (xu - xl) / M;
349: h[1] = (yu - yl) / N;
350: h[2] = (zu - zl) / P;
351: PetscCall(DMDAGetDepthStratum(dm, 0, &vStart, &vEnd));
352: PetscCall(DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV));
353: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
354: PetscCall(PetscSectionSetNumFields(section, 1));
355: PetscCall(PetscSectionSetFieldComponents(section, 0, dim));
356: PetscCall(PetscSectionSetChart(section, vStart, vEnd));
357: for (v = vStart; v < vEnd; ++v) PetscCall(PetscSectionSetDof(section, v, dim));
358: PetscCall(PetscSectionSetUp(section));
359: PetscCall(PetscSectionGetStorageSize(section, &size));
360: PetscCall(VecCreateSeq(PETSC_COMM_SELF, size, &coordinates));
361: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
362: PetscCall(VecGetArray(coordinates, &coords));
363: for (k = 0; k < nVz; ++k) {
364: PetscInt ind[3], d, off;
366: ind[0] = 0;
367: ind[1] = 0;
368: ind[2] = k + da->zs;
369: for (j = 0; j < nVy; ++j) {
370: ind[1] = j + da->ys;
371: for (i = 0; i < nVx; ++i) {
372: const PetscInt vertex = (k * nVy + j) * nVx + i + vStart;
374: PetscCall(PetscSectionGetOffset(section, vertex, &off));
375: ind[0] = i + da->xs;
376: for (d = 0; d < dim; ++d) coords[off + d] = h[d] * ind[d];
377: }
378: }
379: }
380: PetscCall(VecRestoreArray(coordinates, &coords));
381: PetscCall(DMSetCoordinateSection(dm, PETSC_DETERMINE, section));
382: PetscCall(DMSetCoordinatesLocal(dm, coordinates));
383: PetscCall(PetscSectionDestroy(§ion));
384: PetscCall(VecDestroy(&coordinates));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: /*@C
389: DMDAGetArray - Gets a work array for a `DMDA`
391: Input Parameters:
392: + da - a `DMDA`
393: - ghosted - do you want arrays for the ghosted or nonghosted patch
395: Output Parameter:
396: . vptr - array data structured
398: Level: advanced
400: Notes:
401: The vector values are NOT initialized and may have garbage in them, so you may need
402: to zero them.
404: Use `DMDARestoreArray()` to return the array
406: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDARestoreArray()`
407: @*/
408: PetscErrorCode DMDAGetArray(DM da, PetscBool ghosted, void *vptr)
409: {
410: PetscInt j, i, xs, ys, xm, ym, zs, zm;
411: char *iarray_start;
412: void **iptr = (void **)vptr;
413: DM_DA *dd = (DM_DA *)da->data;
415: PetscFunctionBegin;
417: if (ghosted) {
418: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
419: if (dd->arrayghostedin[i]) {
420: *iptr = dd->arrayghostedin[i];
421: iarray_start = (char *)dd->startghostedin[i];
422: dd->arrayghostedin[i] = NULL;
423: dd->startghostedin[i] = NULL;
425: goto done;
426: }
427: }
428: xs = dd->Xs;
429: ys = dd->Ys;
430: zs = dd->Zs;
431: xm = dd->Xe - dd->Xs;
432: ym = dd->Ye - dd->Ys;
433: zm = dd->Ze - dd->Zs;
434: } else {
435: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
436: if (dd->arrayin[i]) {
437: *iptr = dd->arrayin[i];
438: iarray_start = (char *)dd->startin[i];
439: dd->arrayin[i] = NULL;
440: dd->startin[i] = NULL;
442: goto done;
443: }
444: }
445: xs = dd->xs;
446: ys = dd->ys;
447: zs = dd->zs;
448: xm = dd->xe - dd->xs;
449: ym = dd->ye - dd->ys;
450: zm = dd->ze - dd->zs;
451: }
453: switch (da->dim) {
454: case 1: {
455: void *ptr;
457: PetscCall(PetscMalloc(xm * sizeof(PetscScalar), &iarray_start));
459: ptr = (void *)((PetscScalar *)iarray_start - xs);
460: *iptr = ptr;
461: break;
462: }
463: case 2: {
464: void **ptr;
466: PetscCall(PetscMalloc((ym + 1) * sizeof(void *) + xm * ym * sizeof(PetscScalar), &iarray_start));
468: ptr = (void **)(iarray_start + xm * ym * sizeof(PetscScalar) - ys * sizeof(void *));
469: for (j = ys; j < ys + ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar) * (xm * (j - ys) - xs);
470: *iptr = (void *)ptr;
471: break;
472: }
473: case 3: {
474: void ***ptr, **bptr;
476: PetscCall(PetscMalloc((zm + 1) * sizeof(void **) + (ym * zm + 1) * sizeof(void *) + xm * ym * zm * sizeof(PetscScalar), &iarray_start));
478: ptr = (void ***)(iarray_start + xm * ym * zm * sizeof(PetscScalar) - zs * sizeof(void *));
479: bptr = (void **)(iarray_start + xm * ym * zm * sizeof(PetscScalar) + zm * sizeof(void **));
480: for (i = zs; i < zs + zm; i++) ptr[i] = bptr + ((i - zs) * ym - ys);
481: for (i = zs; i < zs + zm; i++) {
482: for (j = ys; j < ys + ym; j++) ptr[i][j] = iarray_start + sizeof(PetscScalar) * (xm * ym * (i - zs) + xm * (j - ys) - xs);
483: }
484: *iptr = (void *)ptr;
485: break;
486: }
487: default:
488: SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", da->dim);
489: }
491: done:
492: /* add arrays to the checked out list */
493: if (ghosted) {
494: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
495: if (!dd->arrayghostedout[i]) {
496: dd->arrayghostedout[i] = *iptr;
497: dd->startghostedout[i] = iarray_start;
498: break;
499: }
500: }
501: } else {
502: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
503: if (!dd->arrayout[i]) {
504: dd->arrayout[i] = *iptr;
505: dd->startout[i] = iarray_start;
506: break;
507: }
508: }
509: }
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
513: /*@C
514: DMDARestoreArray - Restores an array for a `DMDA` obtained with `DMDAGetArray()`
516: Input Parameters:
517: + da - information about my local patch
518: . ghosted - do you want arrays for the ghosted or nonghosted patch
519: - vptr - array data structured
521: Level: advanced
523: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetArray()`
524: @*/
525: PetscErrorCode DMDARestoreArray(DM da, PetscBool ghosted, void *vptr)
526: {
527: PetscInt i;
528: void **iptr = (void **)vptr, *iarray_start = NULL;
529: DM_DA *dd = (DM_DA *)da->data;
531: PetscFunctionBegin;
533: if (ghosted) {
534: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
535: if (dd->arrayghostedout[i] == *iptr) {
536: iarray_start = dd->startghostedout[i];
537: dd->arrayghostedout[i] = NULL;
538: dd->startghostedout[i] = NULL;
539: break;
540: }
541: }
542: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
543: if (!dd->arrayghostedin[i]) {
544: dd->arrayghostedin[i] = *iptr;
545: dd->startghostedin[i] = iarray_start;
546: break;
547: }
548: }
549: } else {
550: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
551: if (dd->arrayout[i] == *iptr) {
552: iarray_start = dd->startout[i];
553: dd->arrayout[i] = NULL;
554: dd->startout[i] = NULL;
555: break;
556: }
557: }
558: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
559: if (!dd->arrayin[i]) {
560: dd->arrayin[i] = *iptr;
561: dd->startin[i] = iarray_start;
562: break;
563: }
564: }
565: }
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }