Actual source code: dacorn.c

petsc-3.11.4 2019-09-28
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: .keywords: distributed array, get, component name

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

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

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

 79:    Not collective; names will contain a common value

 81:    Input Parameter:
 82: .  dm - the DMDA object

 84:    Output Parameter:
 85: .  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

 87:    Level: intermediate

 89:    Not supported from Fortran, use DMDAGetFieldName()

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

 93: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMDASetFieldNames()
 94: @*/
 95: PetscErrorCode  DMDAGetFieldNames(DM da,const char * const **names)
 96: {
 97:   DM_DA             *dd = (DM_DA*)da->data;

100:   *names = (const char * const *) dd->fieldname;
101:   return(0);
102: }

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

107:    Logically collective; names must contain a common value

109:    Input Parameters:
110: +  dm - the DMDA object
111: -  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

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

116:    Level: intermediate

118:    Not supported from Fortran, use DMDASetFieldName()

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

122: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMSetUp()
123: @*/
124: PetscErrorCode  DMDASetFieldNames(DM da,const char * const *names)
125: {
127:   DM_DA          *dd = (DM_DA*)da->data;
128:   char           **fieldname;
129:   PetscInt       nf = 0;

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

141: /*@C
142:    DMDAGetFieldName - Gets the names of individual field components in multicomponent
143:    vectors associated with a DMDA.

145:    Not collective; name will contain a common value

147:    Input Parameter:
148: +  da - the distributed array
149: -  nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
150:         number of degrees of freedom per node within the DMDA

152:    Output Parameter:
153: .  names - the name of the field (component)

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

158:   Level: intermediate

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

162: .seealso: DMDASetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMSetUp()
163: @*/
164: PetscErrorCode  DMDAGetFieldName(DM da,PetscInt nf,const char **name)
165: {
166:   DM_DA *dd = (DM_DA*)da->data;

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

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

180:    Logically collective; name must contain a common value

182:    Input Parameters:
183: +  dm - the DM
184: .  nf - coordinate number for the DMDA (0, 1, ... dim-1),
185: -  name - the name of the coordinate

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


191:   Level: intermediate

193:   Not supported from Fortran

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

197: .seealso: DMDAGetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
198: @*/
199: PetscErrorCode DMDASetCoordinateName(DM dm,PetscInt nf,const char name[])
200: {
202:   DM_DA          *dd = (DM_DA*)dm->data;

206:   if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
207:   if (!dd->coordinatename) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"You should call DMSetUp() first");
208:   PetscFree(dd->coordinatename[nf]);
209:   PetscStrallocpy(name,&dd->coordinatename[nf]);
210:   return(0);
211: }

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

216:    Not collective; name will contain a common value

218:    Input Parameter:
219: +  dm - the DM
220: -  nf -  number for the DMDA (0, 1, ... dim-1)

222:    Output Parameter:
223: .  names - the name of the coordinate direction

225:   Notes:
226:     It must be called after having called DMSetUp().
227:     

229:   Level: intermediate

231:   Not supported from Fortran

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

235: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
236: @*/
237: PetscErrorCode DMDAGetCoordinateName(DM dm,PetscInt nf,const char **name)
238: {
239:   DM_DA *dd = (DM_DA*)dm->data;

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

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

254:    Not collective

256:    Input Parameter:
257: .  da - the distributed array

259:    Output Parameters:
260: +  x,y,z - the corner indices (where y and z are optional; these are used
261:            for 2D and 3D problems)
262: -  m,n,p - widths in the corresponding directions (where n and p are optional;
263:            these are used for 2D and 3D problems)

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

272:   Level: beginner

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

276: .seealso: DMDAGetGhostCorners(), DMDAGetOwnershipRanges()
277: @*/
278: PetscErrorCode  DMDAGetCorners(DM da,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
279: {
280:   PetscInt w;
281:   DM_DA    *dd = (DM_DA*)da->data;

285:   /* since the xs, xe ... have all been multiplied by the number of degrees
286:      of freedom per cell, w = dd->w, we divide that out before returning.*/
287:   w = dd->w;
288:   if (x) *x = dd->xs/w + dd->xo;
289:   /* the y and z have NOT been multiplied by w */
290:   if (y) *y = dd->ys + dd->yo;
291:   if (z) *z = dd->zs + dd->zo;
292:   if (m) *m = (dd->xe - dd->xs)/w;
293:   if (n) *n = (dd->ye - dd->ys);
294:   if (p) *p = (dd->ze - dd->zs);
295:   return(0);
296: }

298: /*@
299:    DMDAGetLocalBoundingBox - Returns the local bounding box for the DMDA.

301:    Not collective

303:    Input Parameter:
304: .  dm - the DM

306:    Output Parameters:
307: +  lmin - local minimum coordinates (length dim, optional)
308: -  lmax - local maximim coordinates (length dim, optional)

310:   Level: beginner

312:   Not supported from Fortran

314: .keywords: distributed array, get, coordinates

316: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetBoundingBox()
317: @*/
318: PetscErrorCode DMDAGetLocalBoundingBox(DM dm,PetscReal lmin[],PetscReal lmax[])
319: {
320:   PetscErrorCode    ierr;
321:   Vec               coords = NULL;
322:   PetscInt          dim,i,j;
323:   const PetscScalar *local_coords;
324:   PetscReal         min[3]={PETSC_MAX_REAL,PETSC_MAX_REAL,PETSC_MAX_REAL},max[3]={PETSC_MIN_REAL,PETSC_MIN_REAL,PETSC_MIN_REAL};
325:   PetscInt          N,Ni;

329:   dim  = dm->dim;
330:   DMGetCoordinates(dm,&coords);
331:   if (coords) {
332:     VecGetArrayRead(coords,&local_coords);
333:     VecGetLocalSize(coords,&N);
334:     Ni   = N/dim;
335:     for (i=0; i<Ni; i++) {
336:       for (j=0; j<3; j++) {
337:         min[j] = j < dim ? PetscMin(min[j],PetscRealPart(local_coords[i*dim+j])) : 0;
338:         max[j] = j < dim ? PetscMax(max[j],PetscRealPart(local_coords[i*dim+j])) : 0;
339:       }
340:     }
341:     VecRestoreArrayRead(coords,&local_coords);
342:   } else {                      /* Just use grid indices */
343:     DMDALocalInfo info;
344:     DMDAGetLocalInfo(dm,&info);
345:     min[0] = info.xs;
346:     min[1] = info.ys;
347:     min[2] = info.zs;
348:     max[0] = info.xs + info.xm-1;
349:     max[1] = info.ys + info.ym-1;
350:     max[2] = info.zs + info.zm-1;
351:   }
352:   if (lmin) {PetscMemcpy(lmin,min,dim*sizeof(PetscReal));}
353:   if (lmax) {PetscMemcpy(lmax,max,dim*sizeof(PetscReal));}
354:   return(0);
355: }

357: /*@
358:    DMDAGetBoundingBox - Returns the global bounding box for the DMDA.

360:    Collective

362:    Input Parameter:
363: .  dm - the DM

365:    Output Parameters:
366: +  gmin - global minimum coordinates (length dim, optional)
367: -  gmax - global maximim coordinates (length dim, optional)

369:   Level: beginner

371: .keywords: distributed array, get, coordinates

373: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetLocalBoundingBox()
374: @*/
375: PetscErrorCode DMDAGetBoundingBox(DM dm,PetscReal gmin[],PetscReal gmax[])
376: {
378:   PetscMPIInt    count;
379:   PetscReal      lmin[3],lmax[3];

383:   PetscMPIIntCast(dm->dim,&count);
384:   DMDAGetLocalBoundingBox(dm,lmin,lmax);
385:   if (gmin) {MPIU_Allreduce(lmin,gmin,count,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)dm));}
386:   if (gmax) {MPIU_Allreduce(lmax,gmax,count,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)dm));}
387:   return(0);
388: }

390: /*@
391:    DMDAGetReducedDMDA - Deprecated; use DMDACreateCompatibleDMDA()

393:    Level: deprecated
394: @*/
395: PetscErrorCode DMDAGetReducedDMDA(DM da,PetscInt nfields,DM *nda)
396: {

400:   DMDACreateCompatibleDMDA(da,nfields,nda);
401:   return(0);
402: }

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

407:    Collective

409:    Input Parameters:
410: +  da - the distributed array
411: -  nfields - number of fields in new DMDA

413:    Output Parameter:
414: .  nda - the new DMDA

416:   Level: intermediate

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

420: .seealso: DMDAGetGhostCorners(), DMSetCoordinates(), DMDASetUniformCoordinates(), DMGetCoordinates(), DMDAGetGhostedCoordinates()
421: @*/
422: PetscErrorCode  DMDACreateCompatibleDMDA(DM da,PetscInt nfields,DM *nda)
423: {
424:   PetscErrorCode   ierr;
425:   DM_DA            *dd = (DM_DA*)da->data;
426:   PetscInt         s,m,n,p,M,N,P,dim,Mo,No,Po;
427:   const PetscInt   *lx,*ly,*lz;
428:   DMBoundaryType   bx,by,bz;
429:   DMDAStencilType  stencil_type;
430:   PetscInt         ox,oy,oz;
431:   PetscInt         cl,rl;

434:   dim = da->dim;
435:   M   = dd->M;
436:   N   = dd->N;
437:   P   = dd->P;
438:   m   = dd->m;
439:   n   = dd->n;
440:   p   = dd->p;
441:   s   = dd->s;
442:   bx  = dd->bx;
443:   by  = dd->by;
444:   bz  = dd->bz;

446:   stencil_type = dd->stencil_type;

448:   DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
449:   if (dim == 1) {
450:     DMDACreate1d(PetscObjectComm((PetscObject)da),bx,M,nfields,s,dd->lx,nda);
451:   } else if (dim == 2) {
452:     DMDACreate2d(PetscObjectComm((PetscObject)da),bx,by,stencil_type,M,N,m,n,nfields,s,lx,ly,nda);
453:   } else if (dim == 3) {
454:     DMDACreate3d(PetscObjectComm((PetscObject)da),bx,by,bz,stencil_type,M,N,P,m,n,p,nfields,s,lx,ly,lz,nda);
455:   }
456:   DMSetUp(*nda);
457:   if (da->coordinates) {
458:     PetscObjectReference((PetscObject)da->coordinates);
459:     (*nda)->coordinates = da->coordinates;
460:   }

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

466:   /* allow for getting a reduced DA corresponding to a coarsened DA */
467:   DMGetCoarsenLevel(da,&cl);
468:   DMGetRefineLevel(da,&rl);

470:   (*nda)->levelup   = rl;
471:   (*nda)->leveldown = cl;
472:   return(0);
473: }

475: /*@C
476:    DMDAGetCoordinateArray - Gets an array containing the coordinates of the DMDA

478:    Not collective

480:    Input Parameter:
481: .  dm - the DM

483:    Output Parameter:
484: .  xc - the coordinates

486:   Level: intermediate

488:   Not supported from Fortran

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

492: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDARestoreCoordinateArray()
493: @*/
494: PetscErrorCode DMDAGetCoordinateArray(DM dm,void *xc)
495: {
497:   DM             cdm;
498:   Vec            x;

502:   DMGetCoordinates(dm,&x);
503:   DMGetCoordinateDM(dm,&cdm);
504:   DMDAVecGetArray(cdm,x,xc);
505:   return(0);
506: }

508: /*@C
509:    DMDARestoreCoordinateArray - Sets an array containing the coordinates of the DMDA

511:    Not collective

513:    Input Parameter:
514: +  dm - the DM
515: -  xc - the coordinates

517:   Level: intermediate

519:   Not supported from Fortran

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

523: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDAGetCoordinateArray()
524: @*/
525: PetscErrorCode DMDARestoreCoordinateArray(DM dm,void *xc)
526: {
528:   DM             cdm;
529:   Vec            x;

533:   DMGetCoordinates(dm,&x);
534:   DMGetCoordinateDM(dm,&cdm);
535:   DMDAVecRestoreArray(cdm,x,xc);
536:   return(0);
537: }