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