Actual source code: dacorn.c

petsc-3.9.4 2018-09-11
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:   Notes: It must be called after having called DMSetUp().

 30:   Level: intermediate

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

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

 43:   if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf);
 44:   if (!dd->fieldname) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ORDER,"You should call DMSetUp() first");
 45:   PetscFree(dd->fieldname[nf]);
 46:   PetscStrallocpy(name,&dd->fieldname[nf]);
 47:   return(0);
 48: }

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

 53:    Collective on TS

 55:    Input Parameter:
 56: .  dm - the DMDA object

 58:    Output Parameter:
 59: .  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

 61:    Level: intermediate

 63:    Not supported from Fortran, use DMDAGetFieldName()

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

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

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

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

 81:    Collective on TS

 83:    Input Parameters:
 84: +  dm - the DMDA object
 85: -  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

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

 89:    Level: intermediate

 91:    Not supported from Fortran, use DMDASetFieldName()

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

 95: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMSetUp()
 96: @*/
 97: PetscErrorCode  DMDASetFieldNames(DM da,const char * const *names)
 98: {
100:   DM_DA          *dd = (DM_DA*)da->data;
101:   char           **fieldname;
102:   PetscInt       nf = 0;

105:   if (!dd->fieldname) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ORDER,"You should call DMSetUp() first");
106:   while (names[nf++]) {};
107:   if (nf != dd->w+1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of fields %D",nf-1);
108:   PetscStrArrayallocpy(names,&fieldname);
109:   PetscStrArrayDestroy(&dd->fieldname);
110:   dd->fieldname = fieldname;
111:   return(0);
112: }

114: /*@C
115:    DMDAGetFieldName - Gets the names of individual field components in multicomponent
116:    vectors associated with a DMDA.

118:    Not Collective

120:    Input Parameter:
121: +  da - the distributed array
122: -  nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
123:         number of degrees of freedom per node within the DMDA

125:    Output Parameter:
126: .  names - the name of the field (component)

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

130:   Level: intermediate

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

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

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

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

152:    Not Collective

154:    Input Parameters:
155: +  dm - the DM
156: .  nf - coordinate number for the DMDA (0, 1, ... dim-1),
157: -  name - the name of the coordinate

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

161:   Level: intermediate

163:   Not supported from Fortran

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

167: .seealso: DMDAGetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
168: @*/
169: PetscErrorCode DMDASetCoordinateName(DM dm,PetscInt nf,const char name[])
170: {
172:   DM_DA          *dd = (DM_DA*)dm->data;

176:   if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
177:   if (!dd->coordinatename) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"You should call DMSetUp() first");
178:   PetscFree(dd->coordinatename[nf]);
179:   PetscStrallocpy(name,&dd->coordinatename[nf]);
180:   return(0);
181: }

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

186:    Not Collective

188:    Input Parameter:
189: +  dm - the DM
190: -  nf -  number for the DMDA (0, 1, ... dim-1)

192:    Output Parameter:
193: .  names - the name of the coordinate direction

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

197:   Level: intermediate

199:   Not supported from Fortran

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

203: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
204: @*/
205: PetscErrorCode DMDAGetCoordinateName(DM dm,PetscInt nf,const char **name)
206: {
207:   DM_DA *dd = (DM_DA*)dm->data;

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

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

222:    Not Collective

224:    Input Parameter:
225: .  da - the distributed array

227:    Output Parameters:
228: +  x,y,z - the corner indices (where y and z are optional; these are used
229:            for 2D and 3D problems)
230: -  m,n,p - widths in the corresponding directions (where n and p are optional;
231:            these are used for 2D and 3D problems)

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

240:   Level: beginner

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

244: .seealso: DMDAGetGhostCorners(), DMDAGetOwnershipRanges()
245: @*/
246: PetscErrorCode  DMDAGetCorners(DM da,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
247: {
248:   PetscInt w;
249:   DM_DA    *dd = (DM_DA*)da->data;

253:   /* since the xs, xe ... have all been multiplied by the number of degrees
254:      of freedom per cell, w = dd->w, we divide that out before returning.*/
255:   w = dd->w;
256:   if (x) *x = dd->xs/w + dd->xo;
257:   /* the y and z have NOT been multiplied by w */
258:   if (y) *y = dd->ys + dd->yo;
259:   if (z) *z = dd->zs + dd->zo;
260:   if (m) *m = (dd->xe - dd->xs)/w;
261:   if (n) *n = (dd->ye - dd->ys);
262:   if (p) *p = (dd->ze - dd->zs);
263:   return(0);
264: }

266: /*@
267:    DMDAGetLocalBoundingBox - Returns the local bounding box for the DMDA.

269:    Not Collective

271:    Input Parameter:
272: .  dm - the DM

274:    Output Parameters:
275: +  lmin - local minimum coordinates (length dim, optional)
276: -  lmax - local maximim coordinates (length dim, optional)

278:   Level: beginner

280:   Not supported from Fortran

282: .keywords: distributed array, get, coordinates

284: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetBoundingBox()
285: @*/
286: PetscErrorCode DMDAGetLocalBoundingBox(DM dm,PetscReal lmin[],PetscReal lmax[])
287: {
288:   PetscErrorCode    ierr;
289:   Vec               coords = NULL;
290:   PetscInt          dim,i,j;
291:   const PetscScalar *local_coords;
292:   PetscReal         min[3]={PETSC_MAX_REAL,PETSC_MAX_REAL,PETSC_MAX_REAL},max[3]={PETSC_MIN_REAL,PETSC_MIN_REAL,PETSC_MIN_REAL};
293:   PetscInt          N,Ni;

297:   dim  = dm->dim;
298:   DMGetCoordinates(dm,&coords);
299:   if (coords) {
300:     VecGetArrayRead(coords,&local_coords);
301:     VecGetLocalSize(coords,&N);
302:     Ni   = N/dim;
303:     for (i=0; i<Ni; i++) {
304:       for (j=0; j<3; j++) {
305:         min[j] = j < dim ? PetscMin(min[j],PetscRealPart(local_coords[i*dim+j])) : 0;
306:         max[j] = j < dim ? PetscMax(max[j],PetscRealPart(local_coords[i*dim+j])) : 0;
307:       }
308:     }
309:     VecRestoreArrayRead(coords,&local_coords);
310:   } else {                      /* Just use grid indices */
311:     DMDALocalInfo info;
312:     DMDAGetLocalInfo(dm,&info);
313:     min[0] = info.xs;
314:     min[1] = info.ys;
315:     min[2] = info.zs;
316:     max[0] = info.xs + info.xm-1;
317:     max[1] = info.ys + info.ym-1;
318:     max[2] = info.zs + info.zm-1;
319:   }
320:   if (lmin) {PetscMemcpy(lmin,min,dim*sizeof(PetscReal));}
321:   if (lmax) {PetscMemcpy(lmax,max,dim*sizeof(PetscReal));}
322:   return(0);
323: }

325: /*@
326:    DMDAGetBoundingBox - Returns the global bounding box for the DMDA.

328:    Collective on DMDA

330:    Input Parameter:
331: .  dm - the DM

333:    Output Parameters:
334: +  gmin - global minimum coordinates (length dim, optional)
335: -  gmax - global maximim coordinates (length dim, optional)

337:   Level: beginner

339: .keywords: distributed array, get, coordinates

341: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetLocalBoundingBox()
342: @*/
343: PetscErrorCode DMDAGetBoundingBox(DM dm,PetscReal gmin[],PetscReal gmax[])
344: {
346:   PetscMPIInt    count;
347:   PetscReal      lmin[3],lmax[3];

351:   PetscMPIIntCast(dm->dim,&count);
352:   DMDAGetLocalBoundingBox(dm,lmin,lmax);
353:   if (gmin) {MPIU_Allreduce(lmin,gmin,count,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)dm));}
354:   if (gmax) {MPIU_Allreduce(lmax,gmax,count,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)dm));}
355:   return(0);
356: }

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

361:    Collective on DMDA

363:    Input Parameters:
364: +  da - the distributed array
365: -  nfields - number of fields in new DMDA

367:    Output Parameter:
368: .  nda - the new DMDA

370:   Level: intermediate

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

374: .seealso: DMDAGetGhostCorners(), DMSetCoordinates(), DMDASetUniformCoordinates(), DMGetCoordinates(), DMDAGetGhostedCoordinates()
375: @*/
376: PetscErrorCode  DMDAGetReducedDMDA(DM da,PetscInt nfields,DM *nda)
377: {
378:   PetscErrorCode   ierr;
379:   DM_DA            *dd = (DM_DA*)da->data;
380:   PetscInt         s,m,n,p,M,N,P,dim,Mo,No,Po;
381:   const PetscInt   *lx,*ly,*lz;
382:   DMBoundaryType   bx,by,bz;
383:   DMDAStencilType  stencil_type;
384:   PetscInt         ox,oy,oz;
385:   PetscInt         cl,rl;

388:   dim = da->dim;
389:   M   = dd->M;
390:   N   = dd->N;
391:   P   = dd->P;
392:   m   = dd->m;
393:   n   = dd->n;
394:   p   = dd->p;
395:   s   = dd->s;
396:   bx  = dd->bx;
397:   by  = dd->by;
398:   bz  = dd->bz;

400:   stencil_type = dd->stencil_type;

402:   DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
403:   if (dim == 1) {
404:     DMDACreate1d(PetscObjectComm((PetscObject)da),bx,M,nfields,s,dd->lx,nda);
405:   } else if (dim == 2) {
406:     DMDACreate2d(PetscObjectComm((PetscObject)da),bx,by,stencil_type,M,N,m,n,nfields,s,lx,ly,nda);
407:   } else if (dim == 3) {
408:     DMDACreate3d(PetscObjectComm((PetscObject)da),bx,by,bz,stencil_type,M,N,P,m,n,p,nfields,s,lx,ly,lz,nda);
409:   }
410:   DMSetUp(*nda);
411:   if (da->coordinates) {
412:     PetscObjectReference((PetscObject)da->coordinates);
413:     (*nda)->coordinates = da->coordinates;
414:   }

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

420:   /* allow for getting a reduced DA corresponding to a coarsened DA */
421:   DMGetCoarsenLevel(da,&cl);
422:   DMGetRefineLevel(da,&rl);

424:   (*nda)->levelup   = rl;
425:   (*nda)->leveldown = cl;
426:   return(0);
427: }

429: /*@C
430:    DMDAGetCoordinateArray - Gets an array containing the coordinates of the DMDA

432:    Not Collective

434:    Input Parameter:
435: .  dm - the DM

437:    Output Parameter:
438: .  xc - the coordinates

440:   Level: intermediate

442:   Not supported from Fortran

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

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

456:   DMGetCoordinates(dm,&x);
457:   DMGetCoordinateDM(dm,&cdm);
458:   DMDAVecGetArray(cdm,x,xc);
459:   return(0);
460: }

462: /*@C
463:    DMDARestoreCoordinateArray - Sets an array containing the coordinates of the DMDA

465:    Not Collective

467:    Input Parameter:
468: +  dm - the DM
469: -  xc - the coordinates

471:   Level: intermediate

473:   Not supported from Fortran

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

477: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDAGetCoordinateArray()
478: @*/
479: PetscErrorCode DMDARestoreCoordinateArray(DM dm,void *xc)
480: {
482:   DM             cdm;
483:   Vec            x;

487:   DMGetCoordinates(dm,&x);
488:   DMGetCoordinateDM(dm,&cdm);
489:   DMDAVecRestoreArray(cdm,x,xc);
490:   return(0);
491: }