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), &section));
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(&section));
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: }