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:   VecGetDM(vec, &da);
 28:   DMDAGetGhostCorners(da, 0, 0, 0, &m, &n, 0);

 30:   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:   PetscArraycpy(mxGetPr(mat), array, n * m);
 37:   PetscObjectName(obj);
 38:   engPutVariable((Engine *)mengine, obj->name, mat);

 40:   VecRestoreArray(vec, &array);
 41:   return 0;
 42: }
 43: #endif

 45: PetscErrorCode DMCreateLocalVector_DA(DM da, Vec *g)
 46: {
 47:   DM_DA *dd = (DM_DA *)da->data;

 51:   VecCreate(PETSC_COMM_SELF, g);
 52:   VecSetSizes(*g, dd->nlocal, PETSC_DETERMINE);
 53:   VecSetBlockSize(*g, dd->w);
 54:   VecSetType(*g, da->vectype);
 55:   if (dd->nlocal < da->bind_below) {
 56:     VecSetBindingPropagates(*g, PETSC_TRUE);
 57:     VecBindToCPU(*g, PETSC_TRUE);
 58:   }
 59:   VecSetDM(*g, da);
 60: #if defined(PETSC_HAVE_MATLAB)
 61:   if (dd->w == 1 && da->dim == 2) PetscObjectComposeFunction((PetscObject)*g, "PetscMatlabEnginePut_C", VecMatlabEnginePut_DA2d);
 62: #endif
 63:   return 0;
 64: }

 66: /*@
 67:   DMDAGetNumCells - Get the number of cells in the local piece of the `DMDA`. This includes ghost cells.

 69:   Input Parameter:
 70: . dm - The `DMDA` object

 72:   Output Parameters:
 73: + numCellsX - The number of local cells in the x-direction
 74: . numCellsY - The number of local cells in the y-direction
 75: . numCellsZ - The number of local cells in the z-direction
 76: - numCells - The number of local cells

 78:   Level: developer

 80: .seealso: `DM`, `DMDA`, `DMDAGetCellPoint()`
 81: @*/
 82: PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells)
 83: {
 84:   DM_DA         *da  = (DM_DA *)dm->data;
 85:   const PetscInt dim = dm->dim;
 86:   const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
 87:   const PetscInt nC = (mx) * (dim > 1 ? (my) * (dim > 2 ? (mz) : 1) : 1);

 90:   if (numCellsX) {
 92:     *numCellsX = mx;
 93:   }
 94:   if (numCellsY) {
 96:     *numCellsY = my;
 97:   }
 98:   if (numCellsZ) {
100:     *numCellsZ = mz;
101:   }
102:   if (numCells) {
104:     *numCells = nC;
105:   }
106:   return 0;
107: }

109: /*@
110:   DMDAGetCellPoint - Get the DM point corresponding to the tuple (i, j, k) in the `DMDA`

112:   Input Parameters:
113: + dm - The `DMDA` object
114: - i,j,k - The global indices for the cell

116:   Output Parameters:
117: . point - The local `DM` point

119:   Level: developer

121: .seealso: `DM`, `DMDA`, `DMDAGetNumCells()`
122: @*/
123: PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
124: {
125:   const PetscInt dim = dm->dim;
126:   DMDALocalInfo  info;

130:   DMDAGetLocalInfo(dm, &info);
134:   *point = i + (dim > 1 ? (j + (dim > 2 ? k * info.gym : 0)) * info.gxm : 0);
135:   return 0;
136: }

138: PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
139: {
140:   DM_DA         *da  = (DM_DA *)dm->data;
141:   const PetscInt dim = dm->dim;
142:   const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
143:   const PetscInt nVx = mx + 1;
144:   const PetscInt nVy = dim > 1 ? (my + 1) : 1;
145:   const PetscInt nVz = dim > 2 ? (mz + 1) : 1;
146:   const PetscInt nV  = nVx * nVy * nVz;

148:   if (numVerticesX) {
150:     *numVerticesX = nVx;
151:   }
152:   if (numVerticesY) {
154:     *numVerticesY = nVy;
155:   }
156:   if (numVerticesZ) {
158:     *numVerticesZ = nVz;
159:   }
160:   if (numVertices) {
162:     *numVertices = nV;
163:   }
164:   return 0;
165: }

167: PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
168: {
169:   DM_DA         *da  = (DM_DA *)dm->data;
170:   const PetscInt dim = dm->dim;
171:   const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
172:   const PetscInt nxF = (dim > 1 ? (my) * (dim > 2 ? (mz) : 1) : 1);
173:   const PetscInt nXF = (mx + 1) * nxF;
174:   const PetscInt nyF = mx * (dim > 2 ? mz : 1);
175:   const PetscInt nYF = dim > 1 ? (my + 1) * nyF : 0;
176:   const PetscInt nzF = mx * (dim > 1 ? my : 0);
177:   const PetscInt nZF = dim > 2 ? (mz + 1) * nzF : 0;

179:   if (numXFacesX) {
181:     *numXFacesX = nxF;
182:   }
183:   if (numXFaces) {
185:     *numXFaces = nXF;
186:   }
187:   if (numYFacesY) {
189:     *numYFacesY = nyF;
190:   }
191:   if (numYFaces) {
193:     *numYFaces = nYF;
194:   }
195:   if (numZFacesZ) {
197:     *numZFacesZ = nzF;
198:   }
199:   if (numZFaces) {
201:     *numZFaces = nZF;
202:   }
203:   return 0;
204: }

206: PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
207: {
208:   const PetscInt dim = dm->dim;
209:   PetscInt       nC, nV, nXF, nYF, nZF;

213:   DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
214:   DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
215:   DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
216:   if (height == 0) {
217:     /* Cells */
218:     if (pStart) *pStart = 0;
219:     if (pEnd) *pEnd = nC;
220:   } else if (height == 1) {
221:     /* Faces */
222:     if (pStart) *pStart = nC + nV;
223:     if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
224:   } else if (height == dim) {
225:     /* Vertices */
226:     if (pStart) *pStart = nC;
227:     if (pEnd) *pEnd = nC + nV;
228:   } else if (height < 0) {
229:     /* All points */
230:     if (pStart) *pStart = 0;
231:     if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
232:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %" PetscInt_FMT " in the DA", height);
233:   return 0;
234: }

236: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
237: {
238:   const PetscInt dim = dm->dim;
239:   PetscInt       nC, nV, nXF, nYF, nZF;

243:   DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
244:   DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
245:   DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
246:   if (depth == dim) {
247:     /* Cells */
248:     if (pStart) *pStart = 0;
249:     if (pEnd) *pEnd = nC;
250:   } else if (depth == dim - 1) {
251:     /* Faces */
252:     if (pStart) *pStart = nC + nV;
253:     if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
254:   } else if (depth == 0) {
255:     /* Vertices */
256:     if (pStart) *pStart = nC;
257:     if (pEnd) *pEnd = nC + nV;
258:   } else if (depth < 0) {
259:     /* All points */
260:     if (pStart) *pStart = 0;
261:     if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
262:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %" PetscInt_FMT " in the DA", depth);
263:   return 0;
264: }

266: PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
267: {
268:   const PetscInt dim = dm->dim;
269:   PetscInt       nC, nV, nXF, nYF, nZF;

271:   *coneSize = 0;
272:   DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
273:   DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
274:   DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
275:   switch (dim) {
276:   case 2:
277:     if (p >= 0) {
278:       if (p < nC) {
279:         *coneSize = 4;
280:       } else if (p < nC + nV) {
281:         *coneSize = 0;
282:       } else if (p < nC + nV + nXF + nYF + nZF) {
283:         *coneSize = 2;
284:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", p, nC + nV + nXF + nYF + nZF);
285:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %" PetscInt_FMT " is invalid", p);
286:     break;
287:   case 3:
288:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
289:   }
290:   return 0;
291: }

293: PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
294: {
295:   const PetscInt dim = dm->dim;
296:   PetscInt       nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;

298:   if (!cone) DMGetWorkArray(dm, 6, MPIU_INT, cone);
299:   DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);
300:   DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
301:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
302:   switch (dim) {
303:   case 2:
304:     if (p >= 0) {
305:       if (p < nC) {
306:         const PetscInt cy = p / nCx;
307:         const PetscInt cx = p % nCx;

309:         (*cone)[0] = cy * nxF + cx + nC + nV;
310:         (*cone)[1] = cx * nyF + cy + nyF + nC + nV + nXF;
311:         (*cone)[2] = cy * nxF + cx + nxF + nC + nV;
312:         (*cone)[3] = cx * nyF + cy + nC + nV + nXF;
313:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
314:       } else if (p < nC + nV) {
315:       } else if (p < nC + nV + nXF) {
316:         const PetscInt fy = (p - nC + nV) / nxF;
317:         const PetscInt fx = (p - nC + nV) % nxF;

319:         (*cone)[0] = fy * nVx + fx + nC;
320:         (*cone)[1] = fy * nVx + fx + 1 + nC;
321:       } else if (p < nC + nV + nXF + nYF) {
322:         const PetscInt fx = (p - nC + nV + nXF) / nyF;
323:         const PetscInt fy = (p - nC + nV + nXF) % nyF;

325:         (*cone)[0] = fy * nVx + fx + nC;
326:         (*cone)[1] = fy * nVx + fx + nVx + nC;
327:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", p, nC + nV + nXF + nYF + nZF);
328:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %" PetscInt_FMT " is invalid", p);
329:     break;
330:   case 3:
331:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
332:   }
333:   return 0;
334: }

336: PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
337: {
338:   DMGetWorkArray(dm, 6, MPIU_INT, cone);
339:   return 0;
340: }

342: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
343: {
344:   DM_DA       *da = (DM_DA *)dm->data;
345:   Vec          coordinates;
346:   PetscSection section;
347:   PetscScalar *coords;
348:   PetscReal    h[3];
349:   PetscInt     dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;

352:   DMDAGetInfo(dm, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
354:   h[0] = (xu - xl) / M;
355:   h[1] = (yu - yl) / N;
356:   h[2] = (zu - zl) / P;
357:   DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);
358:   DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
359:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
360:   PetscSectionSetNumFields(section, 1);
361:   PetscSectionSetFieldComponents(section, 0, dim);
362:   PetscSectionSetChart(section, vStart, vEnd);
363:   for (v = vStart; v < vEnd; ++v) PetscSectionSetDof(section, v, dim);
364:   PetscSectionSetUp(section);
365:   PetscSectionGetStorageSize(section, &size);
366:   VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);
367:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
368:   VecGetArray(coordinates, &coords);
369:   for (k = 0; k < nVz; ++k) {
370:     PetscInt ind[3], d, off;

372:     ind[0] = 0;
373:     ind[1] = 0;
374:     ind[2] = k + da->zs;
375:     for (j = 0; j < nVy; ++j) {
376:       ind[1] = j + da->ys;
377:       for (i = 0; i < nVx; ++i) {
378:         const PetscInt vertex = (k * nVy + j) * nVx + i + vStart;

380:         PetscSectionGetOffset(section, vertex, &off);
381:         ind[0] = i + da->xs;
382:         for (d = 0; d < dim; ++d) coords[off + d] = h[d] * ind[d];
383:       }
384:     }
385:   }
386:   VecRestoreArray(coordinates, &coords);
387:   DMSetCoordinateSection(dm, PETSC_DETERMINE, section);
388:   DMSetCoordinatesLocal(dm, coordinates);
389:   PetscSectionDestroy(&section);
390:   VecDestroy(&coordinates);
391:   return 0;
392: }

394: /* ------------------------------------------------------------------- */

396: /*@C
397:      DMDAGetArray - Gets a work array for a `DMDA`

399:     Input Parameters:
400: +    da - information about my local patch
401: -    ghosted - do you want arrays for the ghosted or nonghosted patch

403:     Output Parameters:
404: .    vptr - array data structured

406:   Level: advanced

408:   Note:
409:    The vector values are NOT initialized and may have garbage in them, so you may need
410:    to zero them.

412: .seealso: `DM`, `DMDA`, `DMDARestoreArray()`
413: @*/
414: PetscErrorCode DMDAGetArray(DM da, PetscBool ghosted, void *vptr)
415: {
416:   PetscInt j, i, xs, ys, xm, ym, zs, zm;
417:   char    *iarray_start;
418:   void   **iptr = (void **)vptr;
419:   DM_DA   *dd   = (DM_DA *)da->data;

422:   if (ghosted) {
423:     for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
424:       if (dd->arrayghostedin[i]) {
425:         *iptr                 = dd->arrayghostedin[i];
426:         iarray_start          = (char *)dd->startghostedin[i];
427:         dd->arrayghostedin[i] = NULL;
428:         dd->startghostedin[i] = NULL;

430:         goto done;
431:       }
432:     }
433:     xs = dd->Xs;
434:     ys = dd->Ys;
435:     zs = dd->Zs;
436:     xm = dd->Xe - dd->Xs;
437:     ym = dd->Ye - dd->Ys;
438:     zm = dd->Ze - dd->Zs;
439:   } else {
440:     for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
441:       if (dd->arrayin[i]) {
442:         *iptr          = dd->arrayin[i];
443:         iarray_start   = (char *)dd->startin[i];
444:         dd->arrayin[i] = NULL;
445:         dd->startin[i] = NULL;

447:         goto done;
448:       }
449:     }
450:     xs = dd->xs;
451:     ys = dd->ys;
452:     zs = dd->zs;
453:     xm = dd->xe - dd->xs;
454:     ym = dd->ye - dd->ys;
455:     zm = dd->ze - dd->zs;
456:   }

458:   switch (da->dim) {
459:   case 1: {
460:     void *ptr;

462:     PetscMalloc(xm * sizeof(PetscScalar), &iarray_start);

464:     ptr   = (void *)(iarray_start - xs * sizeof(PetscScalar));
465:     *iptr = (void *)ptr;
466:     break;
467:   }
468:   case 2: {
469:     void **ptr;

471:     PetscMalloc((ym + 1) * sizeof(void *) + xm * ym * sizeof(PetscScalar), &iarray_start);

473:     ptr = (void **)(iarray_start + xm * ym * sizeof(PetscScalar) - ys * sizeof(void *));
474:     for (j = ys; j < ys + ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar) * (xm * (j - ys) - xs);
475:     *iptr = (void *)ptr;
476:     break;
477:   }
478:   case 3: {
479:     void ***ptr, **bptr;

481:     PetscMalloc((zm + 1) * sizeof(void **) + (ym * zm + 1) * sizeof(void *) + xm * ym * zm * sizeof(PetscScalar), &iarray_start);

483:     ptr  = (void ***)(iarray_start + xm * ym * zm * sizeof(PetscScalar) - zs * sizeof(void *));
484:     bptr = (void **)(iarray_start + xm * ym * zm * sizeof(PetscScalar) + zm * sizeof(void **));
485:     for (i = zs; i < zs + zm; i++) ptr[i] = bptr + ((i - zs) * ym - ys);
486:     for (i = zs; i < zs + zm; i++) {
487:       for (j = ys; j < ys + ym; j++) ptr[i][j] = iarray_start + sizeof(PetscScalar) * (xm * ym * (i - zs) + xm * (j - ys) - xs);
488:     }
489:     *iptr = (void *)ptr;
490:     break;
491:   }
492:   default:
493:     SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", da->dim);
494:   }

496: done:
497:   /* add arrays to the checked out list */
498:   if (ghosted) {
499:     for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
500:       if (!dd->arrayghostedout[i]) {
501:         dd->arrayghostedout[i] = *iptr;
502:         dd->startghostedout[i] = iarray_start;
503:         break;
504:       }
505:     }
506:   } else {
507:     for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
508:       if (!dd->arrayout[i]) {
509:         dd->arrayout[i] = *iptr;
510:         dd->startout[i] = iarray_start;
511:         break;
512:       }
513:     }
514:   }
515:   return 0;
516: }

518: /*@C
519:      DMDARestoreArray - Restores an array of derivative types for a `DMDA`

521:     Input Parameters:
522: +    da - information about my local patch
523: .    ghosted - do you want arrays for the ghosted or nonghosted patch
524: -    vptr - array data structured

526:      Level: advanced

528: .seealso: `DM`, `DMDA`, `DMDAGetArray()`
529: @*/
530: PetscErrorCode DMDARestoreArray(DM da, PetscBool ghosted, void *vptr)
531: {
532:   PetscInt i;
533:   void   **iptr = (void **)vptr, *iarray_start = NULL;
534:   DM_DA   *dd = (DM_DA *)da->data;

537:   if (ghosted) {
538:     for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
539:       if (dd->arrayghostedout[i] == *iptr) {
540:         iarray_start           = dd->startghostedout[i];
541:         dd->arrayghostedout[i] = NULL;
542:         dd->startghostedout[i] = NULL;
543:         break;
544:       }
545:     }
546:     for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
547:       if (!dd->arrayghostedin[i]) {
548:         dd->arrayghostedin[i] = *iptr;
549:         dd->startghostedin[i] = iarray_start;
550:         break;
551:       }
552:     }
553:   } else {
554:     for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
555:       if (dd->arrayout[i] == *iptr) {
556:         iarray_start    = dd->startout[i];
557:         dd->arrayout[i] = NULL;
558:         dd->startout[i] = NULL;
559:         break;
560:       }
561:     }
562:     for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
563:       if (!dd->arrayin[i]) {
564:         dd->arrayin[i] = *iptr;
565:         dd->startin[i] = iarray_start;
566:         break;
567:       }
568:     }
569:   }
570:   return 0;
571: }