Actual source code: dacorn.c

petsc-master 2017-02-25
Report Typos and Errors

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

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

  8: PetscErrorCode DMCreateCoordinateDM_DA(DM dm, DM *cdm)
  9: {
 12:   DMDAGetReducedDMDA(dm,dm->dim,cdm);
 13:   return(0);
 14: }

 16: /*@C
 17:    DMDASetFieldName - Sets the names of individual field components in multicomponent
 18:    vectors associated with a DMDA.

 20:    Not Collective

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

 28:   Level: intermediate

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

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

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

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

 50:    Collective on TS

 52:    Input Parameter:
 53: .  dm - the DMDA object

 55:    Output Parameter:
 56: .  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

 58:    Level: intermediate

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

 62: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMDASetFieldNames()
 63: @*/
 64: PetscErrorCode  DMDAGetFieldNames(DM da,const char * const **names)
 65: {
 66:   DM_DA             *dd = (DM_DA*)da->data;

 69:   *names = (const char * const *) dd->fieldname;
 70:   return(0);
 71: }

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

 76:    Collective on TS

 78:    Input Parameters:
 79: +  dm - the DMDA object
 80: -  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

 82:    Level: intermediate

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

 86: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName()
 87: @*/
 88: PetscErrorCode  DMDASetFieldNames(DM da,const char * const *names)
 89: {
 90:   PetscErrorCode    ierr;
 91:   DM_DA             *dd = (DM_DA*)da->data;

 94:   PetscStrArrayDestroy(&dd->fieldname);
 95:   PetscStrArrayallocpy(names,&dd->fieldname);
 96:   return(0);
 97: }

 99: /*@C
100:    DMDAGetFieldName - Gets the names of individual field components in multicomponent
101:    vectors associated with a DMDA.

103:    Not Collective

105:    Input Parameter:
106: +  da - the distributed array
107: -  nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
108:         number of degrees of freedom per node within the DMDA

110:    Output Parameter:
111: .  names - the name of the field (component)

113:   Level: intermediate

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

117: .seealso: DMDASetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName()
118: @*/
119: PetscErrorCode  DMDAGetFieldName(DM da,PetscInt nf,const char **name)
120: {
121:   DM_DA *dd = (DM_DA*)da->data;

126:   if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf);
127:   *name = dd->fieldname[nf];
128:   return(0);
129: }

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

134:    Not Collective

136:    Input Parameters:
137: +  dm - the DM
138: .  nf - coordinate number for the DMDA (0, 1, ... dim-1),
139: -  name - the name of the coordinate

141:   Level: intermediate

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

145: .seealso: DMDAGetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName()
146: @*/
147: PetscErrorCode DMDASetCoordinateName(DM dm,PetscInt nf,const char name[])
148: {
150:   DM_DA          *dd = (DM_DA*)dm->data;

154:   if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
155:   PetscFree(dd->coordinatename[nf]);
156:   PetscStrallocpy(name,&dd->coordinatename[nf]);
157:   return(0);
158: }

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

163:    Not Collective

165:    Input Parameter:
166: +  dm - the DM
167: -  nf -  number for the DMDA (0, 1, ... dim-1)

169:    Output Parameter:
170: .  names - the name of the coordinate direction

172:   Level: intermediate

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

176: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName()
177: @*/
178: PetscErrorCode DMDAGetCoordinateName(DM dm,PetscInt nf,const char **name)
179: {
180:   DM_DA *dd = (DM_DA*)dm->data;

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

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

194:    Not Collective

196:    Input Parameter:
197: .  da - the distributed array

199:    Output Parameters:
200: +  x,y,z - the corner indices (where y and z are optional; these are used
201:            for 2D and 3D problems)
202: -  m,n,p - widths in the corresponding directions (where n and p are optional;
203:            these are used for 2D and 3D problems)

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

212:   Level: beginner

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

216: .seealso: DMDAGetGhostCorners(), DMDAGetOwnershipRanges()
217: @*/
218: PetscErrorCode  DMDAGetCorners(DM da,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
219: {
220:   PetscInt w;
221:   DM_DA    *dd = (DM_DA*)da->data;

225:   /* since the xs, xe ... have all been multiplied by the number of degrees
226:      of freedom per cell, w = dd->w, we divide that out before returning.*/
227:   w = dd->w;
228:   if (x) *x = dd->xs/w + dd->xo;
229:   /* the y and z have NOT been multiplied by w */
230:   if (y) *y = dd->ys + dd->yo;
231:   if (z) *z = dd->zs + dd->zo;
232:   if (m) *m = (dd->xe - dd->xs)/w;
233:   if (n) *n = (dd->ye - dd->ys);
234:   if (p) *p = (dd->ze - dd->zs);
235:   return(0);
236: }

238: /*@
239:    DMDAGetLocalBoundingBox - Returns the local bounding box for the DMDA.

241:    Not Collective

243:    Input Parameter:
244: .  dm - the DM

246:    Output Parameters:
247: +  lmin - local minimum coordinates (length dim, optional)
248: -  lmax - local maximim coordinates (length dim, optional)

250:   Level: beginner

252: .keywords: distributed array, get, coordinates

254: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetBoundingBox()
255: @*/
256: PetscErrorCode DMDAGetLocalBoundingBox(DM dm,PetscReal lmin[],PetscReal lmax[])
257: {
258:   PetscErrorCode    ierr;
259:   Vec               coords = NULL;
260:   PetscInt          dim,i,j;
261:   const PetscScalar *local_coords;
262:   PetscReal         min[3]={PETSC_MAX_REAL,PETSC_MAX_REAL,PETSC_MAX_REAL},max[3]={PETSC_MIN_REAL,PETSC_MIN_REAL,PETSC_MIN_REAL};
263:   PetscInt          N,Ni;

267:   dim  = dm->dim;
268:   DMGetCoordinates(dm,&coords);
269:   if (coords) {
270:     VecGetArrayRead(coords,&local_coords);
271:     VecGetLocalSize(coords,&N);
272:     Ni   = N/dim;
273:     for (i=0; i<Ni; i++) {
274:       for (j=0; j<3; j++) {
275:         min[j] = j < dim ? PetscMin(min[j],PetscRealPart(local_coords[i*dim+j])) : 0;
276:         max[j] = j < dim ? PetscMax(max[j],PetscRealPart(local_coords[i*dim+j])) : 0;
277:       }
278:     }
279:     VecRestoreArrayRead(coords,&local_coords);
280:   } else {                      /* Just use grid indices */
281:     DMDALocalInfo info;
282:     DMDAGetLocalInfo(dm,&info);
283:     min[0] = info.xs;
284:     min[1] = info.ys;
285:     min[2] = info.zs;
286:     max[0] = info.xs + info.xm-1;
287:     max[1] = info.ys + info.ym-1;
288:     max[2] = info.zs + info.zm-1;
289:   }
290:   if (lmin) {PetscMemcpy(lmin,min,dim*sizeof(PetscReal));}
291:   if (lmax) {PetscMemcpy(lmax,max,dim*sizeof(PetscReal));}
292:   return(0);
293: }

295: /*@
296:    DMDAGetBoundingBox - Returns the global bounding box for the DMDA.

298:    Collective on DMDA

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

303:    Output Parameters:
304: +  gmin - global minimum coordinates (length dim, optional)
305: -  gmax - global maximim coordinates (length dim, optional)

307:   Level: beginner

309: .keywords: distributed array, get, coordinates

311: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetLocalBoundingBox()
312: @*/
313: PetscErrorCode DMDAGetBoundingBox(DM dm,PetscReal gmin[],PetscReal gmax[])
314: {
316:   PetscMPIInt    count;
317:   PetscReal      lmin[3],lmax[3];

321:   PetscMPIIntCast(dm->dim,&count);
322:   DMDAGetLocalBoundingBox(dm,lmin,lmax);
323:   if (gmin) {MPIU_Allreduce(lmin,gmin,count,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)dm));}
324:   if (gmax) {MPIU_Allreduce(lmax,gmax,count,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)dm));}
325:   return(0);
326: }

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

331:    Collective on DMDA

333:    Input Parameters:
334: +  da - the distributed array
335: -  nfields - number of fields in new DMDA

337:    Output Parameter:
338: .  nda - the new DMDA

340:   Level: intermediate

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

344: .seealso: DMDAGetGhostCorners(), DMSetCoordinates(), DMDASetUniformCoordinates(), DMGetCoordinates(), DMDAGetGhostedCoordinates()
345: @*/
346: PetscErrorCode  DMDAGetReducedDMDA(DM da,PetscInt nfields,DM *nda)
347: {
348:   PetscErrorCode   ierr;
349:   DM_DA            *dd = (DM_DA*)da->data;
350:   PetscInt         s,m,n,p,M,N,P,dim,Mo,No,Po;
351:   const PetscInt   *lx,*ly,*lz;
352:   DMBoundaryType   bx,by,bz;
353:   DMDAStencilType  stencil_type;
354:   PetscInt         ox,oy,oz;
355:   PetscInt         cl,rl;

358:   dim = da->dim;
359:   M   = dd->M;
360:   N   = dd->N;
361:   P   = dd->P;
362:   m   = dd->m;
363:   n   = dd->n;
364:   p   = dd->p;
365:   s   = dd->s;
366:   bx  = dd->bx;
367:   by  = dd->by;
368:   bz  = dd->bz;

370:   stencil_type = dd->stencil_type;

372:   DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
373:   if (dim == 1) {
374:     DMDACreate1d(PetscObjectComm((PetscObject)da),bx,M,nfields,s,dd->lx,nda);
375:   } else if (dim == 2) {
376:     DMDACreate2d(PetscObjectComm((PetscObject)da),bx,by,stencil_type,M,N,m,n,nfields,s,lx,ly,nda);
377:   } else if (dim == 3) {
378:     DMDACreate3d(PetscObjectComm((PetscObject)da),bx,by,bz,stencil_type,M,N,P,m,n,p,nfields,s,lx,ly,lz,nda);
379:   }
380:   DMSetUp(*nda);
381:   if (da->coordinates) {
382:     PetscObjectReference((PetscObject)da->coordinates);
383:     (*nda)->coordinates = da->coordinates;
384:   }

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

390:   /* allow for getting a reduced DA corresponding to a coarsened DA */
391:   DMGetCoarsenLevel(da,&cl);
392:   DMGetRefineLevel(da,&rl);

394:   (*nda)->levelup   = rl;
395:   (*nda)->leveldown = cl;
396:   return(0);
397: }

399: /*@C
400:    DMDAGetCoordinateArray - Gets an array containing the coordinates of the DMDA

402:    Not Collective

404:    Input Parameter:
405: .  dm - the DM

407:    Output Parameter:
408: .  xc - the coordinates

410:   Level: intermediate

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

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

424:   DMGetCoordinates(dm,&x);
425:   DMGetCoordinateDM(dm,&cdm);
426:   DMDAVecGetArray(cdm,x,xc);
427:   return(0);
428: }

430: /*@C
431:    DMDARestoreCoordinateArray - Sets an array containing the coordinates of the DMDA

433:    Not Collective

435:    Input Parameter:
436: +  dm - the DM
437: -  xc - the coordinates

439:   Level: intermediate

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

443: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDAGetCoordinateArray()
444: @*/
445: PetscErrorCode DMDARestoreCoordinateArray(DM dm,void *xc)
446: {
448:   DM             cdm;
449:   Vec            x;

453:   DMGetCoordinates(dm,&x);
454:   DMGetCoordinateDM(dm,&cdm);
455:   DMDAVecRestoreArray(cdm,x,xc);
456:   return(0);
457: }