Actual source code: dacorn.c

      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:   DMGetBoundingBox(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().

    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(), DMStagGetCorners()
    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: PetscErrorCode DMGetLocalBoundingIndices_DMDA(DM dm, PetscReal lmin[], PetscReal lmax[])
    285: {
    286:   DMDALocalInfo  info;

    290:   DMDAGetLocalInfo(dm, &info);
    291:   lmin[0] = info.xs;
    292:   lmin[1] = info.ys;
    293:   lmin[2] = info.zs;
    294:   lmax[0] = info.xs + info.xm-1;
    295:   lmax[1] = info.ys + info.ym-1;
    296:   lmax[2] = info.zs + info.zm-1;
    297:   return(0);
    298: }

    300: /*@
    301:    DMDAGetReducedDMDA - Deprecated; use DMDACreateCompatibleDMDA()

    303:    Level: deprecated
    304: @*/
    305: PetscErrorCode DMDAGetReducedDMDA(DM da,PetscInt nfields,DM *nda)
    306: {

    310:   DMDACreateCompatibleDMDA(da,nfields,nda);
    311:   return(0);
    312: }

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

    317:    Collective

    319:    Input Parameters:
    320: +  da - the distributed array
    321: -  nfields - number of fields in new DMDA

    323:    Output Parameter:
    324: .  nda - the new DMDA

    326:   Level: intermediate

    328: .seealso: DMDAGetGhostCorners(), DMSetCoordinates(), DMDASetUniformCoordinates(), DMGetCoordinates(), DMDAGetGhostedCoordinates(), DMStagCreateCompatibleDMStag()
    329: @*/
    330: PetscErrorCode  DMDACreateCompatibleDMDA(DM da,PetscInt nfields,DM *nda)
    331: {
    332:   PetscErrorCode   ierr;
    333:   DM_DA            *dd = (DM_DA*)da->data;
    334:   PetscInt         s,m,n,p,M,N,P,dim,Mo,No,Po;
    335:   const PetscInt   *lx,*ly,*lz;
    336:   DMBoundaryType   bx,by,bz;
    337:   DMDAStencilType  stencil_type;
    338:   PetscInt         ox,oy,oz;
    339:   PetscInt         cl,rl;

    342:   dim = da->dim;
    343:   M   = dd->M;
    344:   N   = dd->N;
    345:   P   = dd->P;
    346:   m   = dd->m;
    347:   n   = dd->n;
    348:   p   = dd->p;
    349:   s   = dd->s;
    350:   bx  = dd->bx;
    351:   by  = dd->by;
    352:   bz  = dd->bz;

    354:   stencil_type = dd->stencil_type;

    356:   DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
    357:   if (dim == 1) {
    358:     DMDACreate1d(PetscObjectComm((PetscObject)da),bx,M,nfields,s,dd->lx,nda);
    359:   } else if (dim == 2) {
    360:     DMDACreate2d(PetscObjectComm((PetscObject)da),bx,by,stencil_type,M,N,m,n,nfields,s,lx,ly,nda);
    361:   } else if (dim == 3) {
    362:     DMDACreate3d(PetscObjectComm((PetscObject)da),bx,by,bz,stencil_type,M,N,P,m,n,p,nfields,s,lx,ly,lz,nda);
    363:   }
    364:   DMSetUp(*nda);
    365:   if (da->coordinates) {
    366:     PetscObjectReference((PetscObject)da->coordinates);
    367:     (*nda)->coordinates = da->coordinates;
    368:   }

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

    374:   /* allow for getting a reduced DA corresponding to a coarsened DA */
    375:   DMGetCoarsenLevel(da,&cl);
    376:   DMGetRefineLevel(da,&rl);

    378:   (*nda)->levelup   = rl;
    379:   (*nda)->leveldown = cl;
    380:   return(0);
    381: }

    383: /*@C
    384:    DMDAGetCoordinateArray - Gets an array containing the coordinates of the DMDA

    386:    Not collective

    388:    Input Parameter:
    389: .  dm - the DM

    391:    Output Parameter:
    392: .  xc - the coordinates

    394:   Level: intermediate

    396:   Not supported from Fortran

    398: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDARestoreCoordinateArray()
    399: @*/
    400: PetscErrorCode DMDAGetCoordinateArray(DM dm,void *xc)
    401: {
    403:   DM             cdm;
    404:   Vec            x;

    408:   DMGetCoordinates(dm,&x);
    409:   DMGetCoordinateDM(dm,&cdm);
    410:   DMDAVecGetArray(cdm,x,xc);
    411:   return(0);
    412: }

    414: /*@C
    415:    DMDARestoreCoordinateArray - Sets an array containing the coordinates of the DMDA

    417:    Not collective

    419:    Input Parameter:
    420: +  dm - the DM
    421: -  xc - the coordinates

    423:   Level: intermediate

    425:   Not supported from Fortran

    427: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDAGetCoordinateArray()
    428: @*/
    429: PetscErrorCode DMDARestoreCoordinateArray(DM dm,void *xc)
    430: {
    432:   DM             cdm;
    433:   Vec            x;

    437:   DMGetCoordinates(dm,&x);
    438:   DMGetCoordinateDM(dm,&cdm);
    439:   DMDAVecRestoreArray(cdm,x,xc);
    440:   return(0);
    441: }