Actual source code: dacorn.c

petsc-dev 2014-08-22
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: {
 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()
 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:    DMDAGetFieldName - Gets the names of individual field components in multicomponent
 55:    vectors associated with a DMDA.

 57:    Not Collective

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

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

 67:   Level: intermediate

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

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

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

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

 90:    Not Collective

 92:    Input Parameters:
 93: +  dm - the DM
 94: .  nf - coordinate number for the DMDA (0, 1, ... dim-1),
 95: -  name - the name of the coordinate

 97:   Level: intermediate

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

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

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

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

121:    Not Collective

123:    Input Parameter:
124: +  dm - the DM
125: -  nf -  number for the DMDA (0, 1, ... dim-1)

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

130:   Level: intermediate

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

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

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

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

154:    Not Collective

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

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

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

172:   Level: beginner

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

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

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

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

200:    Not Collective

202:    Input Parameter:
203: .  dm - the DM

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

209:   Level: beginner

211: .keywords: distributed array, get, coordinates

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

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

256: /*@
257:    DMDAGetBoundingBox - Returns the global bounding box for the DMDA.

259:    Collective on DMDA

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

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

268:   Level: beginner

270: .keywords: distributed array, get, coordinates

272: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetLocalBoundingBox()
273: @*/
274: PetscErrorCode DMDAGetBoundingBox(DM dm,PetscReal gmin[],PetscReal gmax[])
275: {
277:   PetscMPIInt    count;
278:   PetscReal      lmin[3],lmax[3];

282:   PetscMPIIntCast(dm->dim,&count);
283:   DMDAGetLocalBoundingBox(dm,lmin,lmax);
284:   if (gmin) {MPI_Allreduce(lmin,gmin,count,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)dm));}
285:   if (gmax) {MPI_Allreduce(lmax,gmax,count,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)dm));}
286:   return(0);
287: }

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

294:    Collective on DMDA

296:    Input Parameter:
297: +  da - the distributed array
298: .  nfields - number of fields in new DMDA

300:    Output Parameter:
301: .  nda - the new DMDA

303:   Level: intermediate

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

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

321:   dim = da->dim;
322:   M   = dd->M;
323:   N   = dd->N;
324:   P   = dd->P;
325:   m   = dd->m;
326:   n   = dd->n;
327:   p   = dd->p;
328:   s   = dd->s;
329:   bx  = dd->bx;
330:   by  = dd->by;
331:   bz  = dd->bz;

333:   stencil_type = dd->stencil_type;

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

346:     (*nda)->coordinates = da->coordinates;
347:   }

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

353:   /* allow for getting a reduced DA corresponding to a coarsened DA */
354:   DMGetCoarsenLevel(da,&cl);
355:   DMGetRefineLevel(da,&rl);

357:   (*nda)->levelup   = rl;
358:   (*nda)->leveldown = cl;
359:   return(0);
360: }