Actual source code: dacorn.c

petsc-3.5.4 2015-05-23
Report Typos and Errors
  2: /*
  3:   Code for manipulating distributed regular arrays in parallel.
  4: */

  6: #include <petsc-private/dmdaimpl.h>    /*I   "petscdmda.h"   I*/

 10: PetscErrorCode DMCreateCoordinateDM_DA(DM dm, DM *cdm)
 11: {
 13:   DM_DA          *da = (DM_DA*) dm->data;
 15:   DMDAGetReducedDMDA(dm,da->dim,cdm);
 16:   return(0);
 17: }

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

 25:    Not Collective

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

 33:   Level: intermediate

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

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

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

 54: /*@C
 55:    DMDAGetFieldName - Gets the names of individual field components in multicomponent
 56:    vectors associated with a DMDA.

 58:    Not Collective

 60:    Input Parameter:
 61: +  da - the distributed array
 62: -  nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
 63:         number of degrees of freedom per node within the DMDA

 65:    Output Parameter:
 66: .  names - the name of the field (component)

 68:   Level: intermediate

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

 72: .seealso: DMDASetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName()
 73: @*/
 74: PetscErrorCode  DMDAGetFieldName(DM da,PetscInt nf,const char **name)
 75: {
 76:   DM_DA *dd = (DM_DA*)da->data;

 81:   if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf);
 82:   *name = dd->fieldname[nf];
 83:   return(0);
 84: }

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

 91:    Not Collective

 93:    Input Parameters:
 94: +  da - the distributed array
 95: .  nf - coordinate number for the DMDA (0, 1, ... dim-1),
 96: -  name - the name of the coordinate

 98:   Level: intermediate

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

102: .seealso: DMDAGetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName()
103: @*/
104: PetscErrorCode  DMDASetCoordinateName(DM da,PetscInt nf,const char name[])
105: {
107:   DM_DA          *dd = (DM_DA*)da->data;

111:   if (nf < 0 || nf >= dd->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
112:   PetscFree(dd->coordinatename[nf]);
113:   PetscStrallocpy(name,&dd->coordinatename[nf]);
114:   return(0);
115: }

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

122:    Not Collective

124:    Input Parameter:
125: +  da - the distributed array
126: -  nf -  number for the DMDA (0, 1, ... dim-1)

128:    Output Parameter:
129: .  names - the name of the coordinate direction

131:   Level: intermediate

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

135: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName()
136: @*/
137: PetscErrorCode  DMDAGetCoordinateName(DM da,PetscInt nf,const char **name)
138: {
139:   DM_DA *dd = (DM_DA*)da->data;

144:   if (nf < 0 || nf >= dd->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
145:   *name = dd->coordinatename[nf];
146:   return(0);
147: }

151: /*@
152:    DMDAGetCorners - Returns the global (x,y,z) indices of the lower left
153:    corner of the local region, excluding ghost points.

155:    Not Collective

157:    Input Parameter:
158: .  da - the distributed array

160:    Output Parameters:
161: +  x,y,z - the corner indices (where y and z are optional; these are used
162:            for 2D and 3D problems)
163: -  m,n,p - widths in the corresponding directions (where n and p are optional;
164:            these are used for 2D and 3D problems)

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

173:   Level: beginner

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

177: .seealso: DMDAGetGhostCorners(), DMDAGetOwnershipRanges()
178: @*/
179: PetscErrorCode  DMDAGetCorners(DM da,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
180: {
181:   PetscInt w;
182:   DM_DA    *dd = (DM_DA*)da->data;

186:   /* since the xs, xe ... have all been multiplied by the number of degrees
187:      of freedom per cell, w = dd->w, we divide that out before returning.*/
188:   w = dd->w;
189:   if (x) *x = dd->xs/w + dd->xo; if (m) *m = (dd->xe - dd->xs)/w;
190:   /* the y and z have NOT been multiplied by w */
191:   if (y) *y = dd->ys + dd->yo; if (n) *n = (dd->ye - dd->ys);
192:   if (z) *z = dd->zs + dd->zo; if (p) *p = (dd->ze - dd->zs);
193:   return(0);
194: }

198: /*@
199:    DMDAGetLocalBoundingBox - Returns the local bounding box for the DMDA.

201:    Not Collective

203:    Input Parameter:
204: .  da - the distributed array

206:    Output Parameters:
207: +  lmin - local minimum coordinates (length dim, optional)
208: -  lmax - local maximim coordinates (length dim, optional)

210:   Level: beginner

212: .keywords: distributed array, get, coordinates

214: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetBoundingBox()
215: @*/
216: PetscErrorCode  DMDAGetLocalBoundingBox(DM da,PetscReal lmin[],PetscReal lmax[])
217: {
218:   PetscErrorCode    ierr;
219:   Vec               coords = NULL;
220:   PetscInt          dim,i,j;
221:   const PetscScalar *local_coords;
222:   PetscReal         min[3]={PETSC_MAX_REAL,PETSC_MAX_REAL,PETSC_MAX_REAL},max[3]={PETSC_MIN_REAL,PETSC_MIN_REAL,PETSC_MIN_REAL};
223:   PetscInt          N,Ni;
224:   DM_DA             *dd = (DM_DA*)da->data;

228:   dim  = dd->dim;
229:   DMGetCoordinates(da,&coords);
230:   if (coords) {
231:     VecGetArrayRead(coords,&local_coords);
232:     VecGetLocalSize(coords,&N);
233:     Ni   = N/dim;
234:     for (i=0; i<Ni; i++) {
235:       for (j=0; j<3; j++) {
236:         min[j] = j < dim ? PetscMin(min[j],PetscRealPart(local_coords[i*dim+j])) : 0;
237:         max[j] = j < dim ? PetscMax(max[j],PetscRealPart(local_coords[i*dim+j])) : 0;
238:       }
239:     }
240:     VecRestoreArrayRead(coords,&local_coords);
241:   } else {                      /* Just use grid indices */
242:     DMDALocalInfo info;
243:     DMDAGetLocalInfo(da,&info);
244:     min[0] = info.xs;
245:     min[1] = info.ys;
246:     min[2] = info.zs;
247:     max[0] = info.xs + info.xm-1;
248:     max[1] = info.ys + info.ym-1;
249:     max[2] = info.zs + info.zm-1;
250:   }
251:   if (lmin) {PetscMemcpy(lmin,min,dim*sizeof(PetscReal));}
252:   if (lmax) {PetscMemcpy(lmax,max,dim*sizeof(PetscReal));}
253:   return(0);
254: }

258: /*@
259:    DMDAGetBoundingBox - Returns the global bounding box for the DMDA.

261:    Collective on DMDA

263:    Input Parameter:
264: .  da - the distributed array

266:    Output Parameters:
267: +  gmin - global minimum coordinates (length dim, optional)
268: -  gmax - global maximim coordinates (length dim, optional)

270:   Level: beginner

272: .keywords: distributed array, get, coordinates

274: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetLocalBoundingBox()
275: @*/
276: PetscErrorCode  DMDAGetBoundingBox(DM da,PetscReal gmin[],PetscReal gmax[])
277: {
279:   PetscMPIInt    count;
280:   PetscReal      lmin[3],lmax[3];
281:   DM_DA          *dd = (DM_DA*)da->data;

285:   PetscMPIIntCast(dd->dim,&count);
286:   DMDAGetLocalBoundingBox(da,lmin,lmax);
287:   if (gmin) {MPI_Allreduce(lmin,gmin,count,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)da));}
288:   if (gmax) {MPI_Allreduce(lmax,gmax,count,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));}
289:   return(0);
290: }

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

297:    Collective on DMDA

299:    Input Parameter:
300: +  da - the distributed array
301: .  nfields - number of fields in new DMDA

303:    Output Parameter:
304: .  nda - the new DMDA

306:   Level: intermediate

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

310: .seealso: DMDAGetGhostCorners(), DMSetCoordinates(), DMDASetUniformCoordinates(), DMGetCoordinates(), DMDAGetGhostedCoordinates()
311: @*/
312: PetscErrorCode  DMDAGetReducedDMDA(DM da,PetscInt nfields,DM *nda)
313: {
314:   PetscErrorCode   ierr;
315:   DM_DA            *dd = (DM_DA*)da->data;
316:   PetscInt         s,m,n,p,M,N,P,dim,Mo,No,Po;
317:   const PetscInt   *lx,*ly,*lz;
318:   DMBoundaryType   bx,by,bz;
319:   DMDAStencilType  stencil_type;
320:   PetscInt         ox,oy,oz;
321:   PetscInt         cl,rl;

324:   dim = dd->dim;
325:   M   = dd->M;
326:   N   = dd->N;
327:   P   = dd->P;
328:   m   = dd->m;
329:   n   = dd->n;
330:   p   = dd->p;
331:   s   = dd->s;
332:   bx  = dd->bx;
333:   by  = dd->by;
334:   bz  = dd->bz;

336:   stencil_type = dd->stencil_type;

338:   DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
339:   if (dim == 1) {
340:     DMDACreate1d(PetscObjectComm((PetscObject)da),bx,M,nfields,s,dd->lx,nda);
341:   } else if (dim == 2) {
342:     DMDACreate2d(PetscObjectComm((PetscObject)da),bx,by,stencil_type,M,N,m,n,nfields,s,lx,ly,nda);
343:   } else if (dim == 3) {
344:     DMDACreate3d(PetscObjectComm((PetscObject)da),bx,by,bz,stencil_type,M,N,P,m,n,p,nfields,s,lx,ly,lz,nda);
345:   }
346:   if (da->coordinates) {
347:     PetscObjectReference((PetscObject)da->coordinates);

349:     (*nda)->coordinates = da->coordinates;
350:   }

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

356:   /* allow for getting a reduced DA corresponding to a coarsened DA */
357:   DMGetCoarsenLevel(da,&cl);
358:   DMGetRefineLevel(da,&rl);

360:   (*nda)->levelup   = rl;
361:   (*nda)->leveldown = cl;
362:   return(0);
363: }