Actual source code: dagetelem.c

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

      4: static PetscErrorCode DMDAGetElements_1D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
      5: {
      7:   DM_DA          *da = (DM_DA*)dm->data;
      8:   PetscInt       i,xs,xe,Xs,Xe;
      9:   PetscInt       cnt=0;

     12:   if (!da->e) {
     13:     PetscInt corners[2];

     15:     if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
     16:     DMDAGetCorners(dm,&xs,NULL,NULL,&xe,NULL,NULL);
     17:     DMDAGetGhostCorners(dm,&Xs,NULL,NULL,&Xe,NULL,NULL);
     18:     xe    += xs; Xe += Xs; if (xs != Xs) xs -= 1;
     19:     da->ne = 1*(xe - xs - 1);
     20:     PetscMalloc1(1 + 2*da->ne,&da->e);
     21:     for (i=xs; i<xe-1; i++) {
     22:       da->e[cnt++] = (i-Xs);
     23:       da->e[cnt++] = (i-Xs+1);
     24:     }
     25:     da->nen = 2;

     27:     corners[0] = (xs  -Xs);
     28:     corners[1] = (xe-1-Xs);
     29:     ISCreateGeneral(PETSC_COMM_SELF,2,corners,PETSC_COPY_VALUES,&da->ecorners);
     30:   }
     31:   *nel = da->ne;
     32:   *nen = da->nen;
     33:   *e   = da->e;
     34:   return(0);
     35: }

     37: static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
     38: {
     40:   DM_DA          *da = (DM_DA*)dm->data;
     41:   PetscInt       i,xs,xe,Xs,Xe;
     42:   PetscInt       j,ys,ye,Ys,Ye;
     43:   PetscInt       cnt=0, cell[4], ns=2;
     44:   PetscInt       c, split[] = {0,1,3,
     45:                                2,3,1};

     48:   if (!da->e) {
     49:     PetscInt corners[4],nn = 0;

     51:     if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");

     53:     switch (da->elementtype) {
     54:     case DMDA_ELEMENT_Q1:
     55:       da->nen = 4;
     56:       break;
     57:     case DMDA_ELEMENT_P1:
     58:       da->nen = 3;
     59:       break;
     60:     default:
     61:       SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unknown element type %d",da->elementtype);
     62:     }
     63:     nn = da->nen;

     65:     if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;}
     66:     if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
     67:     DMDAGetCorners(dm,&xs,&ys,NULL,&xe,&ye,NULL);
     68:     DMDAGetGhostCorners(dm,&Xs,&Ys,NULL,&Xe,&Ye,NULL);
     69:     xe    += xs; Xe += Xs; if (xs != Xs) xs -= 1;
     70:     ye    += ys; Ye += Ys; if (ys != Ys) ys -= 1;
     71:     da->ne = ns*(xe - xs - 1)*(ye - ys - 1);
     72:     PetscMalloc1(1 + nn*da->ne,&da->e);
     73:     for (j=ys; j<ye-1; j++) {
     74:       for (i=xs; i<xe-1; i++) {
     75:         cell[0] = (i-Xs)   + (j-Ys)*(Xe-Xs);
     76:         cell[1] = (i-Xs+1) + (j-Ys)*(Xe-Xs);
     77:         cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs);
     78:         cell[3] = (i-Xs)   + (j-Ys+1)*(Xe-Xs);
     79:         if (da->elementtype == DMDA_ELEMENT_P1) {
     80:           for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
     81:         }
     82:         if (da->elementtype == DMDA_ELEMENT_Q1) {
     83:           for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
     84:         }
     85:       }
     86:     }

     88:     corners[0] = (xs  -Xs) + (ys  -Ys)*(Xe-Xs);
     89:     corners[1] = (xe-1-Xs) + (ys  -Ys)*(Xe-Xs);
     90:     corners[2] = (xs  -Xs) + (ye-1-Ys)*(Xe-Xs);
     91:     corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs);
     92:     ISCreateGeneral(PETSC_COMM_SELF,4,corners,PETSC_COPY_VALUES,&da->ecorners);
     93:   }
     94:   *nel = da->ne;
     95:   *nen = da->nen;
     96:   *e   = da->e;
     97:   return(0);
     98: }

    100: static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
    101: {
    103:   DM_DA          *da = (DM_DA*)dm->data;
    104:   PetscInt       i,xs,xe,Xs,Xe;
    105:   PetscInt       j,ys,ye,Ys,Ye;
    106:   PetscInt       k,zs,ze,Zs,Ze;
    107:   PetscInt       cnt=0, cell[8], ns=6;
    108:   PetscInt       c, split[] = {0,1,3,7,
    109:                                0,1,7,4,
    110:                                1,2,3,7,
    111:                                1,2,7,6,
    112:                                1,4,5,7,
    113:                                1,5,6,7};

    116:   if (!da->e) {
    117:     PetscInt corners[8],nn = 0;

    119:     if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");

    121:     switch (da->elementtype) {
    122:     case DMDA_ELEMENT_Q1:
    123:       da->nen = 8;
    124:       break;
    125:     case DMDA_ELEMENT_P1:
    126:       da->nen = 4;
    127:       break;
    128:     default:
    129:       SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unknown element type %d",da->elementtype);
    130:     }
    131:     nn = da->nen;

    133:     if (da->elementtype == DMDA_ELEMENT_P1) {ns=6;}
    134:     if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
    135:     DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);
    136:     DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
    137:     xe    += xs; Xe += Xs; if (xs != Xs) xs -= 1;
    138:     ye    += ys; Ye += Ys; if (ys != Ys) ys -= 1;
    139:     ze    += zs; Ze += Zs; if (zs != Zs) zs -= 1;
    140:     da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
    141:     PetscMalloc1(1 + nn*da->ne,&da->e);
    142:     for (k=zs; k<ze-1; k++) {
    143:       for (j=ys; j<ye-1; j++) {
    144:         for (i=xs; i<xe-1; i++) {
    145:           cell[0] = (i  -Xs) + (j  -Ys)*(Xe-Xs) + (k  -Zs)*(Xe-Xs)*(Ye-Ys);
    146:           cell[1] = (i+1-Xs) + (j  -Ys)*(Xe-Xs) + (k  -Zs)*(Xe-Xs)*(Ye-Ys);
    147:           cell[2] = (i+1-Xs) + (j+1-Ys)*(Xe-Xs) + (k  -Zs)*(Xe-Xs)*(Ye-Ys);
    148:           cell[3] = (i  -Xs) + (j+1-Ys)*(Xe-Xs) + (k  -Zs)*(Xe-Xs)*(Ye-Ys);
    149:           cell[4] = (i  -Xs) + (j  -Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
    150:           cell[5] = (i+1-Xs) + (j  -Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
    151:           cell[6] = (i+1-Xs) + (j+1-Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
    152:           cell[7] = (i  -Xs) + (j+1-Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
    153:           if (da->elementtype == DMDA_ELEMENT_P1) {
    154:             for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
    155:           }
    156:           if (da->elementtype == DMDA_ELEMENT_Q1) {
    157:             for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
    158:           }
    159:         }
    160:       }
    161:     }

    163:     corners[0] = (xs  -Xs) + (ys  -Ys)*(Xe-Xs) + (zs-  Zs)*(Xe-Xs)*(Ye-Ys);
    164:     corners[1] = (xe-1-Xs) + (ys  -Ys)*(Xe-Xs) + (zs-  Zs)*(Xe-Xs)*(Ye-Ys);
    165:     corners[2] = (xs  -Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-  Zs)*(Xe-Xs)*(Ye-Ys);
    166:     corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-  Zs)*(Xe-Xs)*(Ye-Ys);
    167:     corners[4] = (xs  -Xs) + (ys  -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
    168:     corners[5] = (xe-1-Xs) + (ys  -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
    169:     corners[6] = (xs  -Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
    170:     corners[7] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
    171:     ISCreateGeneral(PETSC_COMM_SELF,8,corners,PETSC_COPY_VALUES,&da->ecorners);
    172:   }
    173:   *nel = da->ne;
    174:   *nen = da->nen;
    175:   *e   = da->e;
    176:   return(0);
    177: }

    179: /*@
    180:    DMDAGetElementsCorners - Returns the global (x,y,z) indices of the lower left
    181:    corner of the non-overlapping decomposition identified by DMDAGetElements()

    183:     Not Collective

    185:    Input Parameter:
    186: .     da - the DM object

    188:    Output Parameters:
    189: +     gx - the x index
    190: .     gy - the y index
    191: -     gz - the z index

    193:    Level: intermediate

    195:    Notes:

    197: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
    198: @*/
    199: PetscErrorCode  DMDAGetElementsCorners(DM da, PetscInt *gx, PetscInt *gy, PetscInt *gz)
    200: {
    201:   PetscInt       xs,Xs;
    202:   PetscInt       ys,Ys;
    203:   PetscInt       zs,Zs;
    204:   PetscBool      isda;

    212:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
    213:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
    214:   DMDAGetCorners(da,&xs,&ys,&zs,NULL,NULL,NULL);
    215:   DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
    216:   if (xs != Xs) xs -= 1;
    217:   if (ys != Ys) ys -= 1;
    218:   if (zs != Zs) zs -= 1;
    219:   if (gx) *gx  = xs;
    220:   if (gy) *gy  = ys;
    221:   if (gz) *gz  = zs;
    222:   return(0);
    223: }

    225: /*@
    226:       DMDAGetElementsSizes - Gets the local number of elements per direction for the non-overlapping decomposition identified by DMDAGetElements()

    228:     Not Collective

    230:    Input Parameter:
    231: .     da - the DM object

    233:    Output Parameters:
    234: +     mx - number of local elements in x-direction
    235: .     my - number of local elements in y-direction
    236: -     mz - number of local elements in z-direction

    238:    Level: intermediate

    240:    Notes:
    241:     It returns the same number of elements, irrespective of the DMDAElementType

    243: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements
    244: @*/
    245: PetscErrorCode  DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz)
    246: {
    247:   PetscInt       xs,xe,Xs;
    248:   PetscInt       ys,ye,Ys;
    249:   PetscInt       zs,ze,Zs;
    250:   PetscInt       dim;
    251:   PetscBool      isda;

    259:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
    260:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
    261:   DMDAGetCorners(da,&xs,&ys,&zs,&xe,&ye,&ze);
    262:   DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
    263:   xe  += xs; if (xs != Xs) xs -= 1;
    264:   ye  += ys; if (ys != Ys) ys -= 1;
    265:   ze  += zs; if (zs != Zs) zs -= 1;
    266:   if (mx) *mx  = 0;
    267:   if (my) *my  = 0;
    268:   if (mz) *mz  = 0;
    269:   DMGetDimension(da,&dim);
    270:   switch (dim) {
    271:   case 3:
    272:     if (mz) *mz = ze - zs - 1;
    273:   case 2:
    274:     if (my) *my = ye - ys - 1;
    275:   case 1:
    276:     if (mx) *mx = xe - xs - 1;
    277:     break;
    278:   }
    279:   return(0);
    280: }

    282: /*@
    283:       DMDASetElementType - Sets the element type to be returned by DMDAGetElements()

    285:     Not Collective

    287:    Input Parameter:
    288: .     da - the DMDA object

    290:    Output Parameters:
    291: .     etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1

    293:    Level: intermediate

    295: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
    296: @*/
    297: PetscErrorCode  DMDASetElementType(DM da, DMDAElementType etype)
    298: {
    299:   DM_DA          *dd = (DM_DA*)da->data;
    301:   PetscBool      isda;

    306:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
    307:   if (!isda) return(0);
    308:   if (dd->elementtype != etype) {
    309:     PetscFree(dd->e);
    310:     ISDestroy(&dd->ecorners);

    312:     dd->elementtype = etype;
    313:     dd->ne          = 0;
    314:     dd->nen         = 0;
    315:     dd->e           = NULL;
    316:   }
    317:   return(0);
    318: }

    320: /*@
    321:       DMDAGetElementType - Gets the element type to be returned by DMDAGetElements()

    323:     Not Collective

    325:    Input Parameter:
    326: .     da - the DMDA object

    328:    Output Parameters:
    329: .     etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1

    331:    Level: intermediate

    333: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
    334: @*/
    335: PetscErrorCode  DMDAGetElementType(DM da, DMDAElementType *etype)
    336: {
    337:   DM_DA          *dd = (DM_DA*)da->data;
    339:   PetscBool      isda;

    344:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
    345:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
    346:   *etype = dd->elementtype;
    347:   return(0);
    348: }

    350: /*@C
    351:       DMDAGetElements - Gets an array containing the indices (in local coordinates)
    352:                  of all the local elements

    354:     Not Collective

    356:    Input Parameter:
    357: .     dm - the DM object

    359:    Output Parameters:
    360: +     nel - number of local elements
    361: .     nen - number of element nodes
    362: -     e - the local indices of the elements' vertices

    364:    Level: intermediate

    366:    Notes:
    367:      Call DMDARestoreElements() once you have finished accessing the elements.

    369:      Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes.

    371:      If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result.

    373:      Not supported in Fortran

    375: .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
    376: @*/
    377: PetscErrorCode  DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
    378: {
    379:   PetscInt       dim;
    381:   DM_DA          *dd = (DM_DA*)dm->data;
    382:   PetscBool      isda;

    389:   PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
    390:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
    391:   if (dd->stencil_type == DMDA_STENCIL_STAR) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMDAGetElements() requires you use a stencil type of DMDA_STENCIL_BOX");
    392:   DMGetDimension(dm, &dim);
    393:   if (dd->e) {
    394:     *nel = dd->ne;
    395:     *nen = dd->nen;
    396:     *e   = dd->e;
    397:     return(0);
    398:   }
    399:   if (dim==-1) {
    400:     *nel = 0; *nen = 0; *e = NULL;
    401:   } else if (dim==1) {
    402:     DMDAGetElements_1D(dm,nel,nen,e);
    403:   } else if (dim==2) {
    404:     DMDAGetElements_2D(dm,nel,nen,e);
    405:   } else if (dim==3) {
    406:     DMDAGetElements_3D(dm,nel,nen,e);
    407:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
    408:   return(0);
    409: }

    411: /*@
    412:       DMDAGetSubdomainCornersIS - Gets an index set containing the corner indices (in local coordinates)
    413:                                  of the non-overlapping decomposition identified by DMDAGetElements

    415:     Not Collective

    417:    Input Parameter:
    418: .     dm - the DM object

    420:    Output Parameters:
    421: .     is - the index set

    423:    Level: intermediate

    425:    Notes:
    426:     Call DMDARestoreSubdomainCornersIS() once you have finished accessing the index set.

    428: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElementsCornersIS()
    429: @*/
    430: PetscErrorCode  DMDAGetSubdomainCornersIS(DM dm,IS *is)
    431: {
    433:   DM_DA          *dd = (DM_DA*)dm->data;
    434:   PetscBool      isda;

    439:   PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
    440:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
    441:   if (dd->stencil_type == DMDA_STENCIL_STAR) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMDAGetElement() requires you use a stencil type of DMDA_STENCIL_BOX");
    442:   if (!dd->ecorners) { /* compute elements if not yet done */
    443:     const PetscInt *e;
    444:     PetscInt       nel,nen;

    446:     DMDAGetElements(dm,&nel,&nen,&e);
    447:     DMDARestoreElements(dm,&nel,&nen,&e);
    448:   }
    449:   *is = dd->ecorners;
    450:   return(0);
    451: }

    453: /*@C
    454:       DMDARestoreElements - Restores the array obtained with DMDAGetElements()

    456:     Not Collective

    458:    Input Parameter:
    459: +     dm - the DM object
    460: .     nel - number of local elements
    461: .     nen - number of element nodes
    462: -     e - the local indices of the elements' vertices

    464:    Level: intermediate

    466:    Note: You should not access these values after you have called this routine.

    468:          This restore signals the DMDA object that you no longer need access to the array information.

    470:          Not supported in Fortran

    472: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
    473: @*/
    474: PetscErrorCode  DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
    475: {
    481:   *nel = 0;
    482:   *nen = -1;
    483:   *e = NULL;
    484:   return(0);
    485: }

    487: /*@
    488:       DMDARestoreSubdomainCornersIS - Restores the IS obtained with DMDAGetSubdomainCornersIS()

    490:     Not Collective

    492:    Input Parameter:
    493: +     dm - the DM object
    494: -     is - the index set

    496:    Level: intermediate

    498:    Note:

    500: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetSubdomainCornersIS()
    501: @*/
    502: PetscErrorCode  DMDARestoreSubdomainCornersIS(DM dm,IS *is)
    503: {
    507:   *is = NULL;
    508:   return(0);
    509: }