Actual source code: dacorn.c


  2: /*
  3:   Code for manipulating distributed regular arrays in parallel.
  4: */

  6: #include <petsc/private/dmdaimpl.h>
  7: #include <petscdmfield.h>

  9: PetscErrorCode DMCreateCoordinateDM_DA(DM dm, DM *cdm)
 10: {
 11:   const char *prefix;

 13:   PetscFunctionBegin;
 14:   PetscCall(DMDACreateCompatibleDMDA(dm, dm->dim, cdm));
 15:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
 16:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*cdm, prefix));
 17:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*cdm, "cdm_"));
 18:   PetscFunctionReturn(PETSC_SUCCESS);
 19: }

 21: PetscErrorCode DMCreateCoordinateField_DA(DM dm, DMField *field)
 22: {
 23:   PetscReal   gmin[3], gmax[3];
 24:   PetscScalar corners[24];
 25:   PetscInt    dim;
 26:   PetscInt    i, j;
 27:   DM          cdm;

 29:   PetscFunctionBegin;
 30:   PetscCall(DMGetDimension(dm, &dim));
 31:   /* TODO: this is wrong if coordinates are not rectilinear */
 32:   PetscCall(DMGetBoundingBox(dm, gmin, gmax));
 33:   for (i = 0; i < (1 << dim); i++) {
 34:     for (j = 0; j < dim; j++) corners[i * dim + j] = (i & (1 << j)) ? gmax[j] : gmin[j];
 35:   }
 36:   PetscCall(DMClone(dm, &cdm));
 37:   PetscCall(DMFieldCreateDA(cdm, dim, corners, field));
 38:   PetscCall(DMDestroy(&cdm));
 39:   PetscFunctionReturn(PETSC_SUCCESS);
 40: }

 42: /*@C
 43:    DMDASetFieldName - Sets the names of individual field components in multicomponent
 44:    vectors associated with a `DMDA`.

 46:    Logically Collective; name must contain a common value

 48:    Input Parameters:
 49: +  da - the distributed array
 50: .  nf - field number for the `DMDA` (0, 1, ... dof-1), where dof indicates the
 51:         number of degrees of freedom per node within the `DMDA`
 52: -  names - the name of the field (component)

 54:   Level: intermediate

 56:   Note:
 57:     It must be called after having called `DMSetUp()`.

 59: .seealso: `DM`, `DMDA`, `DMDAGetFieldName()`, `DMDASetCoordinateName()`, `DMDAGetCoordinateName()`, `DMDASetFieldNames()`, `DMSetUp()`
 60: @*/
 61: PetscErrorCode DMDASetFieldName(DM da, PetscInt nf, const char name[])
 62: {
 63:   DM_DA *dd = (DM_DA *)da->data;

 65:   PetscFunctionBegin;
 67:   PetscCheck(nf >= 0 && nf < dd->w, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid field number: %" PetscInt_FMT, nf);
 68:   PetscCheck(dd->fieldname, PetscObjectComm((PetscObject)da), PETSC_ERR_ORDER, "You should call DMSetUp() first");
 69:   PetscCall(PetscFree(dd->fieldname[nf]));
 70:   PetscCall(PetscStrallocpy(name, &dd->fieldname[nf]));
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 72: }

 74: /*@C
 75:    DMDAGetFieldNames - Gets the name of each component in the vector associated with the `DMDA`

 77:    Not Collective; names will contain a common value; No Fortran Support

 79:    Input Parameter:
 80: .  dm - the `DMDA` object

 82:    Output Parameter:
 83: .  names - the names of the components, final string is `NULL`, will have the same number of entries as the dof used in creating the `DMDA`

 85:    Level: intermediate

 87:    Fortran Note:
 88:    Use `DMDAGetFieldName()`

 90: .seealso: `DM`, `DMDA`, `DMDAGetFieldName()`, `DMDASetCoordinateName()`, `DMDAGetCoordinateName()`, `DMDASetFieldName()`, `DMDASetFieldNames()`
 91: @*/
 92: PetscErrorCode DMDAGetFieldNames(DM da, const char *const **names)
 93: {
 94:   DM_DA *dd = (DM_DA *)da->data;

 96:   PetscFunctionBegin;
 97:   *names = (const char *const *)dd->fieldname;
 98:   PetscFunctionReturn(PETSC_SUCCESS);
 99: }

101: /*@C
102:    DMDASetFieldNames - Sets the name of each component in the vector associated with the DMDA

104:    Logically Collective; names must contain a common value; No Fortran Support

106:    Input Parameters:
107: +  dm - the `DMDA` object
108: -  names - the names of the components, final string must be NULL, must have the same number of entries as the dof used in creating the `DMDA`

110:    Level: intermediate

112:    Note:
113:     It must be called after having called `DMSetUp()`.

115:    Fortran Note:
116:    Use `DMDASetFieldName()`

118: .seealso: `DM`, `DMDA`, `DMDAGetFieldName()`, `DMDASetCoordinateName()`, `DMDAGetCoordinateName()`, `DMDASetFieldName()`, `DMSetUp()`
119: @*/
120: PetscErrorCode DMDASetFieldNames(DM da, const char *const *names)
121: {
122:   DM_DA   *dd = (DM_DA *)da->data;
123:   char   **fieldname;
124:   PetscInt nf = 0;

126:   PetscFunctionBegin;
127:   PetscCheck(dd->fieldname, PetscObjectComm((PetscObject)da), PETSC_ERR_ORDER, "You should call DMSetUp() first");
128:   while (names[nf++]) { };
129:   PetscCheck(nf == dd->w + 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of fields %" PetscInt_FMT, nf - 1);
130:   PetscCall(PetscStrArrayallocpy(names, &fieldname));
131:   PetscCall(PetscStrArrayDestroy(&dd->fieldname));
132:   dd->fieldname = fieldname;
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

136: /*@C
137:    DMDAGetFieldName - Gets the names of individual field components in multicomponent
138:    vectors associated with a `DMDA`.

140:    Not Collective; name will contain a common value

142:    Input Parameters:
143: +  da - the distributed array
144: -  nf - field number for the `DMDA` (0, 1, ... dof-1), where dof indicates the
145:         number of degrees of freedom per node within the `DMDA`

147:    Output Parameter:
148: .  names - the name of the field (component)

150:   Level: intermediate

152:   Note:
153:     It must be called after having called `DMSetUp()`.

155: .seealso: `DM`, `DMDA`, `DMDASetFieldName()`, `DMDASetCoordinateName()`, `DMDAGetCoordinateName()`, `DMSetUp()`
156: @*/
157: PetscErrorCode DMDAGetFieldName(DM da, PetscInt nf, const char **name)
158: {
159:   DM_DA *dd = (DM_DA *)da->data;

161:   PetscFunctionBegin;
164:   PetscCheck(nf >= 0 && nf < dd->w, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid field number: %" PetscInt_FMT, nf);
165:   PetscCheck(dd->fieldname, PetscObjectComm((PetscObject)da), PETSC_ERR_ORDER, "You should call DMSetUp() first");
166:   *name = dd->fieldname[nf];
167:   PetscFunctionReturn(PETSC_SUCCESS);
168: }

170: /*@C
171:    DMDASetCoordinateName - Sets the name of the coordinate directions associated with a `DMDA`, for example "x" or "y"

173:    Logically Collective; name must contain a common value; No Fortran Support

175:    Input Parameters:
176: +  dm - the `DMDA`
177: .  nf - coordinate number for the DMDA (0, 1, ... dim-1),
178: -  name - the name of the coordinate

180:   Level: intermediate

182:   Note:
183:     Must be called after having called `DMSetUp()`.

185: .seealso: `DM`, `DMDA`, `DMDAGetCoordinateName()`, `DMDASetFieldName()`, `DMDAGetFieldName()`, `DMSetUp()`
186: @*/
187: PetscErrorCode DMDASetCoordinateName(DM dm, PetscInt nf, const char name[])
188: {
189:   DM_DA *dd = (DM_DA *)dm->data;

191:   PetscFunctionBegin;
193:   PetscCheck(nf >= 0 && nf < dm->dim, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid coordinate number: %" PetscInt_FMT, nf);
194:   PetscCheck(dd->coordinatename, PetscObjectComm((PetscObject)dm), PETSC_ERR_ORDER, "You should call DMSetUp() first");
195:   PetscCall(PetscFree(dd->coordinatename[nf]));
196:   PetscCall(PetscStrallocpy(name, &dd->coordinatename[nf]));
197:   PetscFunctionReturn(PETSC_SUCCESS);
198: }

200: /*@C
201:    DMDAGetCoordinateName - Gets the name of a coordinate direction associated with a `DMDA`.

203:    Not Collective; name will contain a common value; No Fortran Support

205:    Input Parameters:
206: +  dm - the `DMDA`
207: -  nf -  number for the `DMDA` (0, 1, ... dim-1)

209:    Output Parameter:
210: .  names - the name of the coordinate direction

212:   Level: intermediate

214:   Note:
215:     It must be called after having called `DMSetUp()`.

217: .seealso: `DM`, `DMDA`, `DMDASetCoordinateName()`, `DMDASetFieldName()`, `DMDAGetFieldName()`, `DMSetUp()`
218: @*/
219: PetscErrorCode DMDAGetCoordinateName(DM dm, PetscInt nf, const char **name)
220: {
221:   DM_DA *dd = (DM_DA *)dm->data;

223:   PetscFunctionBegin;
226:   PetscCheck(nf >= 0 && nf < dm->dim, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid coordinate number: %" PetscInt_FMT, nf);
227:   PetscCheck(dd->coordinatename, PetscObjectComm((PetscObject)dm), PETSC_ERR_ORDER, "You should call DMSetUp() first");
228:   *name = dd->coordinatename[nf];
229:   PetscFunctionReturn(PETSC_SUCCESS);
230: }

232: /*@C
233:    DMDAGetCorners - Returns the global (x,y,z) indices of the lower left
234:    corner and size of the local region, excluding ghost points.

236:    Not Collective

238:    Input Parameter:
239: .  da - the distributed array

241:    Output Parameters:
242: +  x - the corner index for the first dimension
243: .  y - the corner index for the second dimension (only used in 2D and 3D problems)
244: .  z - the corner index for the third dimension (only used in 3D problems)
245: .  m - the width in the first dimension
246: .  n - the width in the second dimension (only used in 2D and 3D problems)
247: -  p - the width in the third dimension (only used in 3D problems)

249:   Level: beginner

251:    Note:
252:    The corner information is independent of the number of degrees of
253:    freedom per node set with the `DMDACreateXX()` routine. Thus the x, y, z, and
254:    m, n, p can be thought of as coordinates on a logical grid, where each
255:    grid point has (potentially) several degrees of freedom.
256:    Any of y, z, n, and p can be passed in as NULL if not needed.

258: .seealso: `DM`, `DMDA`, `DMDAGetGhostCorners()`, `DMDAGetOwnershipRanges()`, `DMStagGetCorners()`
259: @*/
260: PetscErrorCode DMDAGetCorners(DM da, PetscInt *x, PetscInt *y, PetscInt *z, PetscInt *m, PetscInt *n, PetscInt *p)
261: {
262:   PetscInt w;
263:   DM_DA   *dd = (DM_DA *)da->data;

265:   PetscFunctionBegin;
267:   /* since the xs, xe ... have all been multiplied by the number of degrees
268:      of freedom per cell, w = dd->w, we divide that out before returning.*/
269:   w = dd->w;
270:   if (x) *x = dd->xs / w + dd->xo;
271:   /* the y and z have NOT been multiplied by w */
272:   if (y) *y = dd->ys + dd->yo;
273:   if (z) *z = dd->zs + dd->zo;
274:   if (m) *m = (dd->xe - dd->xs) / w;
275:   if (n) *n = (dd->ye - dd->ys);
276:   if (p) *p = (dd->ze - dd->zs);
277:   PetscFunctionReturn(PETSC_SUCCESS);
278: }

280: PetscErrorCode DMGetLocalBoundingIndices_DMDA(DM dm, PetscReal lmin[], PetscReal lmax[])
281: {
282:   DMDALocalInfo info;

284:   PetscFunctionBegin;
285:   PetscCall(DMDAGetLocalInfo(dm, &info));
286:   lmin[0] = info.xs;
287:   lmin[1] = info.ys;
288:   lmin[2] = info.zs;
289:   lmax[0] = info.xs + info.xm - 1;
290:   lmax[1] = info.ys + info.ym - 1;
291:   lmax[2] = info.zs + info.zm - 1;
292:   PetscFunctionReturn(PETSC_SUCCESS);
293: }

295: /*@
296:    DMDAGetReducedDMDA - Deprecated; use DMDACreateCompatibleDMDA()

298:    Level: deprecated
299: @*/
300: PetscErrorCode DMDAGetReducedDMDA(DM da, PetscInt nfields, DM *nda)
301: {
302:   PetscFunctionBegin;
303:   PetscCall(DMDACreateCompatibleDMDA(da, nfields, nda));
304:   PetscFunctionReturn(PETSC_SUCCESS);
305: }

307: /*@
308:    DMDACreateCompatibleDMDA - Creates a `DMDA` with the same layout but with fewer or more fields

310:    Collective

312:    Input Parameters:
313: +  da - the distributed array
314: -  nfields - number of fields in new `DMDA`

316:    Output Parameter:
317: .  nda - the new `DMDA`

319:   Level: intermediate

321: .seealso: `DM`, `DMDA`, `DMDAGetGhostCorners()`, `DMSetCoordinates()`, `DMDASetUniformCoordinates()`, `DMGetCoordinates()`, `DMDAGetGhostedCoordinates()`, `DMStagCreateCompatibleDMStag()`
322: @*/
323: PetscErrorCode DMDACreateCompatibleDMDA(DM da, PetscInt nfields, DM *nda)
324: {
325:   DM_DA          *dd = (DM_DA *)da->data;
326:   PetscInt        s, m, n, p, M, N, P, dim, Mo, No, Po;
327:   const PetscInt *lx, *ly, *lz;
328:   DMBoundaryType  bx, by, bz;
329:   DMDAStencilType stencil_type;
330:   Vec             coords;
331:   PetscInt        ox, oy, oz;
332:   PetscInt        cl, rl;

334:   PetscFunctionBegin;
335:   dim = da->dim;
336:   M   = dd->M;
337:   N   = dd->N;
338:   P   = dd->P;
339:   m   = dd->m;
340:   n   = dd->n;
341:   p   = dd->p;
342:   s   = dd->s;
343:   bx  = dd->bx;
344:   by  = dd->by;
345:   bz  = dd->bz;

347:   stencil_type = dd->stencil_type;

349:   PetscCall(DMDAGetOwnershipRanges(da, &lx, &ly, &lz));
350:   if (dim == 1) {
351:     PetscCall(DMDACreate1d(PetscObjectComm((PetscObject)da), bx, M, nfields, s, dd->lx, nda));
352:   } else if (dim == 2) {
353:     PetscCall(DMDACreate2d(PetscObjectComm((PetscObject)da), bx, by, stencil_type, M, N, m, n, nfields, s, lx, ly, nda));
354:   } else if (dim == 3) {
355:     PetscCall(DMDACreate3d(PetscObjectComm((PetscObject)da), bx, by, bz, stencil_type, M, N, P, m, n, p, nfields, s, lx, ly, lz, nda));
356:   }
357:   PetscCall(DMSetUp(*nda));
358:   PetscCall(DMGetCoordinates(da, &coords));
359:   PetscCall(DMSetCoordinates(*nda, coords));

361:   /* allow for getting a reduced DA corresponding to a domain decomposition */
362:   PetscCall(DMDAGetOffset(da, &ox, &oy, &oz, &Mo, &No, &Po));
363:   PetscCall(DMDASetOffset(*nda, ox, oy, oz, Mo, No, Po));

365:   /* allow for getting a reduced DA corresponding to a coarsened DA */
366:   PetscCall(DMGetCoarsenLevel(da, &cl));
367:   PetscCall(DMGetRefineLevel(da, &rl));

369:   (*nda)->levelup   = rl;
370:   (*nda)->leveldown = cl;
371:   PetscFunctionReturn(PETSC_SUCCESS);
372: }

374: /*@C
375:    DMDAGetCoordinateArray - Gets an array containing the coordinates of the `DMDA`

377:    Not Collective; No Fortran Support

379:    Input Parameter:
380: .  dm - the `DMDA`

382:    Output Parameter:
383: .  xc - the coordinates

385:   Level: intermediate

387: .seealso: `DM`, `DMDA`, `DMDASetCoordinateName()`, `DMDASetFieldName()`, `DMDAGetFieldName()`, `DMDARestoreCoordinateArray()`
388: @*/
389: PetscErrorCode DMDAGetCoordinateArray(DM dm, void *xc)
390: {
391:   DM  cdm;
392:   Vec x;

394:   PetscFunctionBegin;
396:   PetscCall(DMGetCoordinates(dm, &x));
397:   PetscCall(DMGetCoordinateDM(dm, &cdm));
398:   PetscCall(DMDAVecGetArray(cdm, x, xc));
399:   PetscFunctionReturn(PETSC_SUCCESS);
400: }

402: /*@C
403:    DMDARestoreCoordinateArray - Sets an array containing the coordinates of the `DMDA`

405:    Not Collective; No Fortran Support

407:    Input Parameters:
408: +  dm - the `DMDA`
409: -  xc - the coordinates

411:   Level: intermediate

413: .seealso: `DM`, `DMDA`, `DMDASetCoordinateName()`, `DMDASetFieldName()`, `DMDAGetFieldName()`, `DMDAGetCoordinateArray()`
414: @*/
415: PetscErrorCode DMDARestoreCoordinateArray(DM dm, void *xc)
416: {
417:   DM  cdm;
418:   Vec x;

420:   PetscFunctionBegin;
422:   PetscCall(DMGetCoordinates(dm, &x));
423:   PetscCall(DMGetCoordinateDM(dm, &cdm));
424:   PetscCall(DMDAVecRestoreArray(cdm, x, xc));
425:   PetscFunctionReturn(PETSC_SUCCESS);
426: }