Actual source code: dacorn.c

petsc-master 2018-05-22
Report Typos and Errors

  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: {
 13:   DMDACreateCompatibleDMDA(dm,dm->dim,cdm);
 14:   return(0);
 15: }

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

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

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

 45:    Not Collective

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

 53:   Notes:
 54:     It must be called after having called DMSetUp().

 56:   Level: intermediate

 58: .keywords: distributed array, get, component name

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

 69:   if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf);
 70:   if (!dd->fieldname) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ORDER,"You should call DMSetUp() first");
 71:   PetscFree(dd->fieldname[nf]);
 72:   PetscStrallocpy(name,&dd->fieldname[nf]);
 73:   return(0);
 74: }

 76: /*@C
 77:    DMDAGetFieldNames - Gets the name of each component in the vector associated with the DMDA

 79:    Collective on TS

 81:    Input Parameter:
 82: .  dm - the DMDA object

 84:    Output Parameter:
 85: .  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

 87:    Level: intermediate

 89:    Not supported from Fortran, use DMDAGetFieldName()

 91: .keywords: distributed array, get, component name

 93: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMDASetFieldNames()
 94: @*/
 95: PetscErrorCode  DMDAGetFieldNames(DM da,const char * const **names)
 96: {
 97:   DM_DA             *dd = (DM_DA*)da->data;

100:   *names = (const char * const *) dd->fieldname;
101:   return(0);
102: }

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

107:    Collective on TS

109:    Input Parameters:
110: +  dm - the DMDA object
111: -  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

113:    Notes:
114:     It must be called after having called DMSetUp().

116:    Level: intermediate

118:    Not supported from Fortran, use DMDASetFieldName()

120: .keywords: distributed array, get, component name

122: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMSetUp()
123: @*/
124: PetscErrorCode  DMDASetFieldNames(DM da,const char * const *names)
125: {
127:   DM_DA          *dd = (DM_DA*)da->data;
128:   char           **fieldname;
129:   PetscInt       nf = 0;

132:   if (!dd->fieldname) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ORDER,"You should call DMSetUp() first");
133:   while (names[nf++]) {};
134:   if (nf != dd->w+1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of fields %D",nf-1);
135:   PetscStrArrayallocpy(names,&fieldname);
136:   PetscStrArrayDestroy(&dd->fieldname);
137:   dd->fieldname = fieldname;
138:   return(0);
139: }

141: /*@C
142:    DMDAGetFieldName - Gets the names of individual field components in multicomponent
143:    vectors associated with a DMDA.

145:    Not Collective

147:    Input Parameter:
148: +  da - the distributed array
149: -  nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
150:         number of degrees of freedom per node within the DMDA

152:    Output Parameter:
153: .  names - the name of the field (component)

155:   Notes:
156:     It must be called after having called DMSetUp().

158:   Level: intermediate

160: .keywords: distributed array, get, component name

162: .seealso: DMDASetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMSetUp()
163: @*/
164: PetscErrorCode  DMDAGetFieldName(DM da,PetscInt nf,const char **name)
165: {
166:   DM_DA *dd = (DM_DA*)da->data;

171:   if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf);
172:   if (!dd->fieldname) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ORDER,"You should call DMSetUp() first");
173:   *name = dd->fieldname[nf];
174:   return(0);
175: }

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

180:    Not Collective

182:    Input Parameters:
183: +  dm - the DM
184: .  nf - coordinate number for the DMDA (0, 1, ... dim-1),
185: -  name - the name of the coordinate

187:   Notes:
188:     It must be called after having called DMSetUp().

190:   Level: intermediate

192:   Not supported from Fortran

194: .keywords: distributed array, get, component name

196: .seealso: DMDAGetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
197: @*/
198: PetscErrorCode DMDASetCoordinateName(DM dm,PetscInt nf,const char name[])
199: {
201:   DM_DA          *dd = (DM_DA*)dm->data;

205:   if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
206:   if (!dd->coordinatename) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"You should call DMSetUp() first");
207:   PetscFree(dd->coordinatename[nf]);
208:   PetscStrallocpy(name,&dd->coordinatename[nf]);
209:   return(0);
210: }

212: /*@C
213:    DMDAGetCoordinateName - Gets the name of a coodinate direction associated with a DMDA.

215:    Not Collective

217:    Input Parameter:
218: +  dm - the DM
219: -  nf -  number for the DMDA (0, 1, ... dim-1)

221:    Output Parameter:
222: .  names - the name of the coordinate direction

224:   Notes:
225:     It must be called after having called DMSetUp().

227:   Level: intermediate

229:   Not supported from Fortran

231: .keywords: distributed array, get, component name

233: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
234: @*/
235: PetscErrorCode DMDAGetCoordinateName(DM dm,PetscInt nf,const char **name)
236: {
237:   DM_DA *dd = (DM_DA*)dm->data;

242:   if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
243:   if (!dd->coordinatename) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"You should call DMSetUp() first");
244:   *name = dd->coordinatename[nf];
245:   return(0);
246: }

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

252:    Not Collective

254:    Input Parameter:
255: .  da - the distributed array

257:    Output Parameters:
258: +  x,y,z - the corner indices (where y and z are optional; these are used
259:            for 2D and 3D problems)
260: -  m,n,p - widths in the corresponding directions (where n and p are optional;
261:            these are used for 2D and 3D problems)

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

270:   Level: beginner

272: .keywords: distributed array, get, corners, nodes, local indices

274: .seealso: DMDAGetGhostCorners(), DMDAGetOwnershipRanges()
275: @*/
276: PetscErrorCode  DMDAGetCorners(DM da,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
277: {
278:   PetscInt w;
279:   DM_DA    *dd = (DM_DA*)da->data;

283:   /* since the xs, xe ... have all been multiplied by the number of degrees
284:      of freedom per cell, w = dd->w, we divide that out before returning.*/
285:   w = dd->w;
286:   if (x) *x = dd->xs/w + dd->xo;
287:   /* the y and z have NOT been multiplied by w */
288:   if (y) *y = dd->ys + dd->yo;
289:   if (z) *z = dd->zs + dd->zo;
290:   if (m) *m = (dd->xe - dd->xs)/w;
291:   if (n) *n = (dd->ye - dd->ys);
292:   if (p) *p = (dd->ze - dd->zs);
293:   return(0);
294: }

296: /*@
297:    DMDAGetLocalBoundingBox - Returns the local bounding box for the DMDA.

299:    Not Collective

301:    Input Parameter:
302: .  dm - the DM

304:    Output Parameters:
305: +  lmin - local minimum coordinates (length dim, optional)
306: -  lmax - local maximim coordinates (length dim, optional)

308:   Level: beginner

310:   Not supported from Fortran

312: .keywords: distributed array, get, coordinates

314: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetBoundingBox()
315: @*/
316: PetscErrorCode DMDAGetLocalBoundingBox(DM dm,PetscReal lmin[],PetscReal lmax[])
317: {
318:   PetscErrorCode    ierr;
319:   Vec               coords = NULL;
320:   PetscInt          dim,i,j;
321:   const PetscScalar *local_coords;
322:   PetscReal         min[3]={PETSC_MAX_REAL,PETSC_MAX_REAL,PETSC_MAX_REAL},max[3]={PETSC_MIN_REAL,PETSC_MIN_REAL,PETSC_MIN_REAL};
323:   PetscInt          N,Ni;

327:   dim  = dm->dim;
328:   DMGetCoordinates(dm,&coords);
329:   if (coords) {
330:     VecGetArrayRead(coords,&local_coords);
331:     VecGetLocalSize(coords,&N);
332:     Ni   = N/dim;
333:     for (i=0; i<Ni; i++) {
334:       for (j=0; j<3; j++) {
335:         min[j] = j < dim ? PetscMin(min[j],PetscRealPart(local_coords[i*dim+j])) : 0;
336:         max[j] = j < dim ? PetscMax(max[j],PetscRealPart(local_coords[i*dim+j])) : 0;
337:       }
338:     }
339:     VecRestoreArrayRead(coords,&local_coords);
340:   } else {                      /* Just use grid indices */
341:     DMDALocalInfo info;
342:     DMDAGetLocalInfo(dm,&info);
343:     min[0] = info.xs;
344:     min[1] = info.ys;
345:     min[2] = info.zs;
346:     max[0] = info.xs + info.xm-1;
347:     max[1] = info.ys + info.ym-1;
348:     max[2] = info.zs + info.zm-1;
349:   }
350:   if (lmin) {PetscMemcpy(lmin,min,dim*sizeof(PetscReal));}
351:   if (lmax) {PetscMemcpy(lmax,max,dim*sizeof(PetscReal));}
352:   return(0);
353: }

355: /*@
356:    DMDAGetBoundingBox - Returns the global bounding box for the DMDA.

358:    Collective on DMDA

360:    Input Parameter:
361: .  dm - the DM

363:    Output Parameters:
364: +  gmin - global minimum coordinates (length dim, optional)
365: -  gmax - global maximim coordinates (length dim, optional)

367:   Level: beginner

369: .keywords: distributed array, get, coordinates

371: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetLocalBoundingBox()
372: @*/
373: PetscErrorCode DMDAGetBoundingBox(DM dm,PetscReal gmin[],PetscReal gmax[])
374: {
376:   PetscMPIInt    count;
377:   PetscReal      lmin[3],lmax[3];

381:   PetscMPIIntCast(dm->dim,&count);
382:   DMDAGetLocalBoundingBox(dm,lmin,lmax);
383:   if (gmin) {MPIU_Allreduce(lmin,gmin,count,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)dm));}
384:   if (gmax) {MPIU_Allreduce(lmax,gmax,count,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)dm));}
385:   return(0);
386: }

388: /*@
389:    DMDAGetReducedDMDA - Deprecated; use DMDACreateCompatibleDMDA()

391:    Level: deprecated
392: @*/
393: PetscErrorCode DMDAGetReducedDMDA(DM da,PetscInt nfields,DM *nda)
394: {

398:   DMDACreateCompatibleDMDA(da,nfields,nda);
399:   return(0);
400: }

402: /*@
403:    DMDACreateCompatibleDMDA - Creates a DMDA with the same layout but with fewer or more fields

405:    Collective on DMDA

407:    Input Parameters:
408: +  da - the distributed array
409: -  nfields - number of fields in new DMDA

411:    Output Parameter:
412: .  nda - the new DMDA

414:   Level: intermediate

416: .keywords: distributed array, get, corners, nodes, local indices, coordinates

418: .seealso: DMDAGetGhostCorners(), DMSetCoordinates(), DMDASetUniformCoordinates(), DMGetCoordinates(), DMDAGetGhostedCoordinates()
419: @*/
420: PetscErrorCode  DMDACreateCompatibleDMDA(DM da,PetscInt nfields,DM *nda)
421: {
422:   PetscErrorCode   ierr;
423:   DM_DA            *dd = (DM_DA*)da->data;
424:   PetscInt         s,m,n,p,M,N,P,dim,Mo,No,Po;
425:   const PetscInt   *lx,*ly,*lz;
426:   DMBoundaryType   bx,by,bz;
427:   DMDAStencilType  stencil_type;
428:   PetscInt         ox,oy,oz;
429:   PetscInt         cl,rl;

432:   dim = da->dim;
433:   M   = dd->M;
434:   N   = dd->N;
435:   P   = dd->P;
436:   m   = dd->m;
437:   n   = dd->n;
438:   p   = dd->p;
439:   s   = dd->s;
440:   bx  = dd->bx;
441:   by  = dd->by;
442:   bz  = dd->bz;

444:   stencil_type = dd->stencil_type;

446:   DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
447:   if (dim == 1) {
448:     DMDACreate1d(PetscObjectComm((PetscObject)da),bx,M,nfields,s,dd->lx,nda);
449:   } else if (dim == 2) {
450:     DMDACreate2d(PetscObjectComm((PetscObject)da),bx,by,stencil_type,M,N,m,n,nfields,s,lx,ly,nda);
451:   } else if (dim == 3) {
452:     DMDACreate3d(PetscObjectComm((PetscObject)da),bx,by,bz,stencil_type,M,N,P,m,n,p,nfields,s,lx,ly,lz,nda);
453:   }
454:   DMSetUp(*nda);
455:   if (da->coordinates) {
456:     PetscObjectReference((PetscObject)da->coordinates);
457:     (*nda)->coordinates = da->coordinates;
458:   }

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

464:   /* allow for getting a reduced DA corresponding to a coarsened DA */
465:   DMGetCoarsenLevel(da,&cl);
466:   DMGetRefineLevel(da,&rl);

468:   (*nda)->levelup   = rl;
469:   (*nda)->leveldown = cl;
470:   return(0);
471: }

473: /*@C
474:    DMDAGetCoordinateArray - Gets an array containing the coordinates of the DMDA

476:    Not Collective

478:    Input Parameter:
479: .  dm - the DM

481:    Output Parameter:
482: .  xc - the coordinates

484:   Level: intermediate

486:   Not supported from Fortran

488: .keywords: distributed array, get, component name

490: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDARestoreCoordinateArray()
491: @*/
492: PetscErrorCode DMDAGetCoordinateArray(DM dm,void *xc)
493: {
495:   DM             cdm;
496:   Vec            x;

500:   DMGetCoordinates(dm,&x);
501:   DMGetCoordinateDM(dm,&cdm);
502:   DMDAVecGetArray(cdm,x,xc);
503:   return(0);
504: }

506: /*@C
507:    DMDARestoreCoordinateArray - Sets an array containing the coordinates of the DMDA

509:    Not Collective

511:    Input Parameter:
512: +  dm - the DM
513: -  xc - the coordinates

515:   Level: intermediate

517:   Not supported from Fortran

519: .keywords: distributed array, get, component name

521: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDAGetCoordinateArray()
522: @*/
523: PetscErrorCode DMDARestoreCoordinateArray(DM dm,void *xc)
524: {
526:   DM             cdm;
527:   Vec            x;

531:   DMGetCoordinates(dm,&x);
532:   DMGetCoordinateDM(dm,&cdm);
533:   DMDAVecRestoreArray(cdm,x,xc);
534:   return(0);
535: }