Actual source code: dacorn.c

petsc-master 2019-08-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:    Logically collective; name must contain a common value

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

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

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

 77:    Not collective; names will contain a common value

 79:    Input Parameter:
 80: .  dm - the DMDA object

 82:    Output Parameter:
 83: .  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

 85:    Level: intermediate

 87:    Not supported from Fortran, use DMDAGetFieldName()

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

 96:   *names = (const char * const *) dd->fieldname;
 97:   return(0);
 98: }

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

103:    Logically collective; names must contain a common value

105:    Input Parameters:
106: +  dm - the DMDA object
107: -  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

109:    Notes:
110:     It must be called after having called DMSetUp().

112:    Level: intermediate

114:    Not supported from Fortran, use DMDASetFieldName()

116: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMSetUp()
117: @*/
118: PetscErrorCode  DMDASetFieldNames(DM da,const char * const *names)
119: {
121:   DM_DA          *dd = (DM_DA*)da->data;
122:   char           **fieldname;
123:   PetscInt       nf = 0;

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

135: /*@C
136:    DMDAGetFieldName - Gets the names of individual field components in multicomponent
137:    vectors associated with a DMDA.

139:    Not collective; name will contain a common value

141:    Input Parameter:
142: +  da - the distributed array
143: -  nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
144:         number of degrees of freedom per node within the DMDA

146:    Output Parameter:
147: .  names - the name of the field (component)

149:   Notes:
150:     It must be called after having called DMSetUp().

152:   Level: intermediate

154: .seealso: DMDASetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMSetUp()
155: @*/
156: PetscErrorCode  DMDAGetFieldName(DM da,PetscInt nf,const char **name)
157: {
158:   DM_DA *dd = (DM_DA*)da->data;

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

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

172:    Logically collective; name must contain a common value

174:    Input Parameters:
175: +  dm - the DM
176: .  nf - coordinate number for the DMDA (0, 1, ... dim-1),
177: -  name - the name of the coordinate

179:   Notes:
180:     It must be called after having called DMSetUp().


183:   Level: intermediate

185:   Not supported from Fortran

187: .seealso: DMDAGetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
188: @*/
189: PetscErrorCode DMDASetCoordinateName(DM dm,PetscInt nf,const char name[])
190: {
192:   DM_DA          *dd = (DM_DA*)dm->data;

196:   if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
197:   if (!dd->coordinatename) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"You should call DMSetUp() first");
198:   PetscFree(dd->coordinatename[nf]);
199:   PetscStrallocpy(name,&dd->coordinatename[nf]);
200:   return(0);
201: }

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

206:    Not collective; name will contain a common value

208:    Input Parameter:
209: +  dm - the DM
210: -  nf -  number for the DMDA (0, 1, ... dim-1)

212:    Output Parameter:
213: .  names - the name of the coordinate direction

215:   Notes:
216:     It must be called after having called DMSetUp().
217:     

219:   Level: intermediate

221:   Not supported from Fortran

223: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
224: @*/
225: PetscErrorCode DMDAGetCoordinateName(DM dm,PetscInt nf,const char **name)
226: {
227:   DM_DA *dd = (DM_DA*)dm->data;

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

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

242:    Not collective

244:    Input Parameter:
245: .  da - the distributed array

247:    Output Parameters:
248: +  x,y,z - the corner indices (where y and z are optional; these are used
249:            for 2D and 3D problems)
250: -  m,n,p - widths in the corresponding directions (where n and p are optional;
251:            these are used for 2D and 3D problems)

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

260:   Level: beginner

262: .seealso: DMDAGetGhostCorners(), DMDAGetOwnershipRanges()
263: @*/
264: PetscErrorCode  DMDAGetCorners(DM da,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
265: {
266:   PetscInt w;
267:   DM_DA    *dd = (DM_DA*)da->data;

271:   /* since the xs, xe ... have all been multiplied by the number of degrees
272:      of freedom per cell, w = dd->w, we divide that out before returning.*/
273:   w = dd->w;
274:   if (x) *x = dd->xs/w + dd->xo;
275:   /* the y and z have NOT been multiplied by w */
276:   if (y) *y = dd->ys + dd->yo;
277:   if (z) *z = dd->zs + dd->zo;
278:   if (m) *m = (dd->xe - dd->xs)/w;
279:   if (n) *n = (dd->ye - dd->ys);
280:   if (p) *p = (dd->ze - dd->zs);
281:   return(0);
282: }

284: /*@
285:    DMDAGetLocalBoundingBox - Returns the local bounding box for the DMDA.

287:    Not collective

289:    Input Parameter:
290: .  dm - the DM

292:    Output Parameters:
293: +  lmin - local minimum coordinates (length dim, optional)
294: -  lmax - local maximim coordinates (length dim, optional)

296:   Level: beginner

298:   Not supported from Fortran

300: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetBoundingBox()
301: @*/
302: PetscErrorCode DMDAGetLocalBoundingBox(DM dm,PetscReal lmin[],PetscReal lmax[])
303: {
304:   PetscErrorCode    ierr;
305:   Vec               coords = NULL;
306:   PetscInt          dim,i,j;
307:   const PetscScalar *local_coords;
308:   PetscReal         min[3]={PETSC_MAX_REAL,PETSC_MAX_REAL,PETSC_MAX_REAL},max[3]={PETSC_MIN_REAL,PETSC_MIN_REAL,PETSC_MIN_REAL};
309:   PetscInt          N,Ni;

313:   dim  = dm->dim;
314:   DMGetCoordinates(dm,&coords);
315:   if (coords) {
316:     VecGetArrayRead(coords,&local_coords);
317:     VecGetLocalSize(coords,&N);
318:     Ni   = N/dim;
319:     for (i=0; i<Ni; i++) {
320:       for (j=0; j<3; j++) {
321:         min[j] = j < dim ? PetscMin(min[j],PetscRealPart(local_coords[i*dim+j])) : 0;
322:         max[j] = j < dim ? PetscMax(max[j],PetscRealPart(local_coords[i*dim+j])) : 0;
323:       }
324:     }
325:     VecRestoreArrayRead(coords,&local_coords);
326:   } else {                      /* Just use grid indices */
327:     DMDALocalInfo info;
328:     DMDAGetLocalInfo(dm,&info);
329:     min[0] = info.xs;
330:     min[1] = info.ys;
331:     min[2] = info.zs;
332:     max[0] = info.xs + info.xm-1;
333:     max[1] = info.ys + info.ym-1;
334:     max[2] = info.zs + info.zm-1;
335:   }
336:   if (lmin) {PetscArraycpy(lmin,min,dim);}
337:   if (lmax) {PetscArraycpy(lmax,max,dim);}
338:   return(0);
339: }

341: /*@
342:    DMDAGetBoundingBox - Returns the global bounding box for the DMDA.

344:    Collective

346:    Input Parameter:
347: .  dm - the DM

349:    Output Parameters:
350: +  gmin - global minimum coordinates (length dim, optional)
351: -  gmax - global maximim coordinates (length dim, optional)

353:   Level: beginner

355: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetLocalBoundingBox()
356: @*/
357: PetscErrorCode DMDAGetBoundingBox(DM dm,PetscReal gmin[],PetscReal gmax[])
358: {
360:   PetscMPIInt    count;
361:   PetscReal      lmin[3],lmax[3];

365:   PetscMPIIntCast(dm->dim,&count);
366:   DMDAGetLocalBoundingBox(dm,lmin,lmax);
367:   if (gmin) {MPIU_Allreduce(lmin,gmin,count,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)dm));}
368:   if (gmax) {MPIU_Allreduce(lmax,gmax,count,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)dm));}
369:   return(0);
370: }

372: /*@
373:    DMDAGetReducedDMDA - Deprecated; use DMDACreateCompatibleDMDA()

375:    Level: deprecated
376: @*/
377: PetscErrorCode DMDAGetReducedDMDA(DM da,PetscInt nfields,DM *nda)
378: {

382:   DMDACreateCompatibleDMDA(da,nfields,nda);
383:   return(0);
384: }

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

389:    Collective

391:    Input Parameters:
392: +  da - the distributed array
393: -  nfields - number of fields in new DMDA

395:    Output Parameter:
396: .  nda - the new DMDA

398:   Level: intermediate

400: .seealso: DMDAGetGhostCorners(), DMSetCoordinates(), DMDASetUniformCoordinates(), DMGetCoordinates(), DMDAGetGhostedCoordinates()
401: @*/
402: PetscErrorCode  DMDACreateCompatibleDMDA(DM da,PetscInt nfields,DM *nda)
403: {
404:   PetscErrorCode   ierr;
405:   DM_DA            *dd = (DM_DA*)da->data;
406:   PetscInt         s,m,n,p,M,N,P,dim,Mo,No,Po;
407:   const PetscInt   *lx,*ly,*lz;
408:   DMBoundaryType   bx,by,bz;
409:   DMDAStencilType  stencil_type;
410:   PetscInt         ox,oy,oz;
411:   PetscInt         cl,rl;

414:   dim = da->dim;
415:   M   = dd->M;
416:   N   = dd->N;
417:   P   = dd->P;
418:   m   = dd->m;
419:   n   = dd->n;
420:   p   = dd->p;
421:   s   = dd->s;
422:   bx  = dd->bx;
423:   by  = dd->by;
424:   bz  = dd->bz;

426:   stencil_type = dd->stencil_type;

428:   DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
429:   if (dim == 1) {
430:     DMDACreate1d(PetscObjectComm((PetscObject)da),bx,M,nfields,s,dd->lx,nda);
431:   } else if (dim == 2) {
432:     DMDACreate2d(PetscObjectComm((PetscObject)da),bx,by,stencil_type,M,N,m,n,nfields,s,lx,ly,nda);
433:   } else if (dim == 3) {
434:     DMDACreate3d(PetscObjectComm((PetscObject)da),bx,by,bz,stencil_type,M,N,P,m,n,p,nfields,s,lx,ly,lz,nda);
435:   }
436:   DMSetUp(*nda);
437:   if (da->coordinates) {
438:     PetscObjectReference((PetscObject)da->coordinates);
439:     (*nda)->coordinates = da->coordinates;
440:   }

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

446:   /* allow for getting a reduced DA corresponding to a coarsened DA */
447:   DMGetCoarsenLevel(da,&cl);
448:   DMGetRefineLevel(da,&rl);

450:   (*nda)->levelup   = rl;
451:   (*nda)->leveldown = cl;
452:   return(0);
453: }

455: /*@C
456:    DMDAGetCoordinateArray - Gets an array containing the coordinates of the DMDA

458:    Not collective

460:    Input Parameter:
461: .  dm - the DM

463:    Output Parameter:
464: .  xc - the coordinates

466:   Level: intermediate

468:   Not supported from Fortran

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

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

486: /*@C
487:    DMDARestoreCoordinateArray - Sets an array containing the coordinates of the DMDA

489:    Not collective

491:    Input Parameter:
492: +  dm - the DM
493: -  xc - the coordinates

495:   Level: intermediate

497:   Not supported from Fortran

499: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDAGetCoordinateArray()
500: @*/
501: PetscErrorCode DMDARestoreCoordinateArray(DM dm,void *xc)
502: {
504:   DM             cdm;
505:   Vec            x;

509:   DMGetCoordinates(dm,&x);
510:   DMGetCoordinateDM(dm,&cdm);
511:   DMDAVecRestoreArray(cdm,x,xc);
512:   return(0);
513: }