Actual source code: dacorn.c

petsc-master 2016-09-25
Report Typos and Errors
  2: /*
  3:   Code for manipulating distributed regular arrays in parallel.
  4: */

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

 10: PetscErrorCode DMCreateCoordinateDM_DA(DM dm, DM *cdm)
 11: {
 14:   DMDAGetReducedDMDA(dm,dm->dim,cdm);
 15:   return(0);
 16: }

 20: /*@C
 21:    DMDASetFieldName - Sets the names of individual field components in multicomponent
 22:    vectors associated with a DMDA.

 24:    Not Collective

 26:    Input Parameters:
 27: +  da - the distributed array
 28: .  nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
 29:         number of degrees of freedom per node within the DMDA
 30: -  names - the name of the field (component)

 32:   Level: intermediate

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

 36: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldNames()
 37: @*/
 38: PetscErrorCode  DMDASetFieldName(DM da,PetscInt nf,const char name[])
 39: {
 41:   DM_DA          *dd = (DM_DA*)da->data;

 45:   if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf);
 46:   PetscFree(dd->fieldname[nf]);
 47:   PetscStrallocpy(name,&dd->fieldname[nf]);
 48:   return(0);
 49: }

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

 56:    Collective on TS

 58:    Input Parameter:
 59: .  dm - the DMDA object

 61:    Output Parameter:
 62: .  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

 64:    Level: intermediate

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

 68: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMDASetFieldNames()
 69: @*/
 70: PetscErrorCode  DMDAGetFieldNames(DM da,const char * const **names)
 71: {
 72:   DM_DA             *dd = (DM_DA*)da->data;

 75:   *names = (const char * const *) dd->fieldname;
 76:   return(0);
 77: }

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

 84:    Collective on TS

 86:    Input Parameters:
 87: +  dm - the DMDA object
 88: -  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

 90:    Level: intermediate

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

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

102:   PetscStrArrayDestroy(&dd->fieldname);
103:   PetscStrArrayallocpy(names,&dd->fieldname);
104:   return(0);
105: }

109: /*@C
110:    DMDAGetFieldName - Gets the names of individual field components in multicomponent
111:    vectors associated with a DMDA.

113:    Not Collective

115:    Input Parameter:
116: +  da - the distributed array
117: -  nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
118:         number of degrees of freedom per node within the DMDA

120:    Output Parameter:
121: .  names - the name of the field (component)

123:   Level: intermediate

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

127: .seealso: DMDASetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName()
128: @*/
129: PetscErrorCode  DMDAGetFieldName(DM da,PetscInt nf,const char **name)
130: {
131:   DM_DA *dd = (DM_DA*)da->data;

136:   if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf);
137:   *name = dd->fieldname[nf];
138:   return(0);
139: }

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

146:    Not Collective

148:    Input Parameters:
149: +  dm - the DM
150: .  nf - coordinate number for the DMDA (0, 1, ... dim-1),
151: -  name - the name of the coordinate

153:   Level: intermediate

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

157: .seealso: DMDAGetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName()
158: @*/
159: PetscErrorCode DMDASetCoordinateName(DM dm,PetscInt nf,const char name[])
160: {
162:   DM_DA          *dd = (DM_DA*)dm->data;

166:   if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
167:   PetscFree(dd->coordinatename[nf]);
168:   PetscStrallocpy(name,&dd->coordinatename[nf]);
169:   return(0);
170: }

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

177:    Not Collective

179:    Input Parameter:
180: +  dm - the DM
181: -  nf -  number for the DMDA (0, 1, ... dim-1)

183:    Output Parameter:
184: .  names - the name of the coordinate direction

186:   Level: intermediate

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

190: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName()
191: @*/
192: PetscErrorCode DMDAGetCoordinateName(DM dm,PetscInt nf,const char **name)
193: {
194:   DM_DA *dd = (DM_DA*)dm->data;

199:   if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
200:   *name = dd->coordinatename[nf];
201:   return(0);
202: }

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

210:    Not Collective

212:    Input Parameter:
213: .  da - the distributed array

215:    Output Parameters:
216: +  x,y,z - the corner indices (where y and z are optional; these are used
217:            for 2D and 3D problems)
218: -  m,n,p - widths in the corresponding directions (where n and p are optional;
219:            these are used for 2D and 3D problems)

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

228:   Level: beginner

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

232: .seealso: DMDAGetGhostCorners(), DMDAGetOwnershipRanges()
233: @*/
234: PetscErrorCode  DMDAGetCorners(DM da,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
235: {
236:   PetscInt w;
237:   DM_DA    *dd = (DM_DA*)da->data;

241:   /* since the xs, xe ... have all been multiplied by the number of degrees
242:      of freedom per cell, w = dd->w, we divide that out before returning.*/
243:   w = dd->w;
244:   if (x) *x = dd->xs/w + dd->xo;
245:   /* the y and z have NOT been multiplied by w */
246:   if (y) *y = dd->ys + dd->yo;
247:   if (z) *z = dd->zs + dd->zo;
248:   if (m) *m = (dd->xe - dd->xs)/w;
249:   if (n) *n = (dd->ye - dd->ys);
250:   if (p) *p = (dd->ze - dd->zs);
251:   return(0);
252: }

256: /*@
257:    DMDAGetLocalBoundingBox - Returns the local bounding box for the DMDA.

259:    Not Collective

261:    Input Parameter:
262: .  dm - the DM

264:    Output Parameters:
265: +  lmin - local minimum coordinates (length dim, optional)
266: -  lmax - local maximim coordinates (length dim, optional)

268:   Level: beginner

270: .keywords: distributed array, get, coordinates

272: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetBoundingBox()
273: @*/
274: PetscErrorCode DMDAGetLocalBoundingBox(DM dm,PetscReal lmin[],PetscReal lmax[])
275: {
276:   PetscErrorCode    ierr;
277:   Vec               coords = NULL;
278:   PetscInt          dim,i,j;
279:   const PetscScalar *local_coords;
280:   PetscReal         min[3]={PETSC_MAX_REAL,PETSC_MAX_REAL,PETSC_MAX_REAL},max[3]={PETSC_MIN_REAL,PETSC_MIN_REAL,PETSC_MIN_REAL};
281:   PetscInt          N,Ni;

285:   dim  = dm->dim;
286:   DMGetCoordinates(dm,&coords);
287:   if (coords) {
288:     VecGetArrayRead(coords,&local_coords);
289:     VecGetLocalSize(coords,&N);
290:     Ni   = N/dim;
291:     for (i=0; i<Ni; i++) {
292:       for (j=0; j<3; j++) {
293:         min[j] = j < dim ? PetscMin(min[j],PetscRealPart(local_coords[i*dim+j])) : 0;
294:         max[j] = j < dim ? PetscMax(max[j],PetscRealPart(local_coords[i*dim+j])) : 0;
295:       }
296:     }
297:     VecRestoreArrayRead(coords,&local_coords);
298:   } else {                      /* Just use grid indices */
299:     DMDALocalInfo info;
300:     DMDAGetLocalInfo(dm,&info);
301:     min[0] = info.xs;
302:     min[1] = info.ys;
303:     min[2] = info.zs;
304:     max[0] = info.xs + info.xm-1;
305:     max[1] = info.ys + info.ym-1;
306:     max[2] = info.zs + info.zm-1;
307:   }
308:   if (lmin) {PetscMemcpy(lmin,min,dim*sizeof(PetscReal));}
309:   if (lmax) {PetscMemcpy(lmax,max,dim*sizeof(PetscReal));}
310:   return(0);
311: }

315: /*@
316:    DMDAGetBoundingBox - Returns the global bounding box for the DMDA.

318:    Collective on DMDA

320:    Input Parameter:
321: .  dm - the DM

323:    Output Parameters:
324: +  gmin - global minimum coordinates (length dim, optional)
325: -  gmax - global maximim coordinates (length dim, optional)

327:   Level: beginner

329: .keywords: distributed array, get, coordinates

331: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetLocalBoundingBox()
332: @*/
333: PetscErrorCode DMDAGetBoundingBox(DM dm,PetscReal gmin[],PetscReal gmax[])
334: {
336:   PetscMPIInt    count;
337:   PetscReal      lmin[3],lmax[3];

341:   PetscMPIIntCast(dm->dim,&count);
342:   DMDAGetLocalBoundingBox(dm,lmin,lmax);
343:   if (gmin) {MPIU_Allreduce(lmin,gmin,count,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)dm));}
344:   if (gmax) {MPIU_Allreduce(lmax,gmax,count,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)dm));}
345:   return(0);
346: }

350: /*@
351:    DMDAGetReducedDMDA - Gets the DMDA with the same layout but with fewer or more fields

353:    Collective on DMDA

355:    Input Parameters:
356: +  da - the distributed array
357: -  nfields - number of fields in new DMDA

359:    Output Parameter:
360: .  nda - the new DMDA

362:   Level: intermediate

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

366: .seealso: DMDAGetGhostCorners(), DMSetCoordinates(), DMDASetUniformCoordinates(), DMGetCoordinates(), DMDAGetGhostedCoordinates()
367: @*/
368: PetscErrorCode  DMDAGetReducedDMDA(DM da,PetscInt nfields,DM *nda)
369: {
370:   PetscErrorCode   ierr;
371:   DM_DA            *dd = (DM_DA*)da->data;
372:   PetscInt         s,m,n,p,M,N,P,dim,Mo,No,Po;
373:   const PetscInt   *lx,*ly,*lz;
374:   DMBoundaryType   bx,by,bz;
375:   DMDAStencilType  stencil_type;
376:   PetscInt         ox,oy,oz;
377:   PetscInt         cl,rl;

380:   dim = da->dim;
381:   M   = dd->M;
382:   N   = dd->N;
383:   P   = dd->P;
384:   m   = dd->m;
385:   n   = dd->n;
386:   p   = dd->p;
387:   s   = dd->s;
388:   bx  = dd->bx;
389:   by  = dd->by;
390:   bz  = dd->bz;

392:   stencil_type = dd->stencil_type;

394:   DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
395:   if (dim == 1) {
396:     DMDACreate1d(PetscObjectComm((PetscObject)da),bx,M,nfields,s,dd->lx,nda);
397:   } else if (dim == 2) {
398:     DMDACreate2d(PetscObjectComm((PetscObject)da),bx,by,stencil_type,M,N,m,n,nfields,s,lx,ly,nda);
399:   } else if (dim == 3) {
400:     DMDACreate3d(PetscObjectComm((PetscObject)da),bx,by,bz,stencil_type,M,N,P,m,n,p,nfields,s,lx,ly,lz,nda);
401:   }
402:   DMSetUp(*nda);
403:   if (da->coordinates) {
404:     PetscObjectReference((PetscObject)da->coordinates);
405:     (*nda)->coordinates = da->coordinates;
406:   }

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

412:   /* allow for getting a reduced DA corresponding to a coarsened DA */
413:   DMGetCoarsenLevel(da,&cl);
414:   DMGetRefineLevel(da,&rl);

416:   (*nda)->levelup   = rl;
417:   (*nda)->leveldown = cl;
418:   return(0);
419: }

423: /*@C
424:    DMDAGetCoordinateArray - Gets an array containing the coordinates of the DMDA

426:    Not Collective

428:    Input Parameter:
429: .  dm - the DM

431:    Output Parameter:
432: .  xc - the coordinates

434:   Level: intermediate

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

438: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDARestoreCoordinateArray()
439: @*/
440: PetscErrorCode DMDAGetCoordinateArray(DM dm,void *xc)
441: {
443:   DM             cdm;
444:   Vec            x;

448:   DMGetCoordinates(dm,&x);
449:   DMGetCoordinateDM(dm,&cdm);
450:   DMDAVecGetArray(cdm,x,xc);
451:   return(0);
452: }

456: /*@C
457:    DMDARestoreCoordinateArray - Sets an array containing the coordinates of the DMDA

459:    Not Collective

461:    Input Parameter:
462: +  dm - the DM
463: -  xc - the coordinates

465:   Level: intermediate

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

469: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDAGetCoordinateArray()
470: @*/
471: PetscErrorCode DMDARestoreCoordinateArray(DM dm,void *xc)
472: {
474:   DM             cdm;
475:   Vec            x;

479:   DMGetCoordinates(dm,&x);
480:   DMGetCoordinateDM(dm,&cdm);
481:   DMDAVecRestoreArray(cdm,x,xc);
482:   return(0);
483: }