Actual source code: stag.c

      1: /*
      2:    Implementation of DMStag, defining dimension-independent functions in the
      3:    DM API. stag1d.c, stag2d.c, and stag3d.c may include dimension-specific
      4:    implementations of DM API functions, and other files here contain additional
      5:    DMStag-specific API functions, as well as internal functions.
      6: */
      7: #include <petsc/private/dmstagimpl.h>
      8: #include <petscsf.h>

     10: static PetscErrorCode DMClone_Stag(DM dm,DM *newdm)
     11: {

     15:   /* Destroy the DM created by generic logic in DMClone() */
     16:   if (*newdm) {
     17:     DMDestroy(newdm);
     18:   }
     19:   DMStagDuplicateWithoutSetup(dm,PetscObjectComm((PetscObject)dm),newdm);
     20:   DMSetUp(*newdm);
     21:   return(0);
     22: }

     24: static PetscErrorCode DMDestroy_Stag(DM dm)
     25: {
     27:   DM_Stag        *stag;
     28:   PetscInt       i;

     31:   stag = (DM_Stag*)dm->data;
     32:   for (i=0; i<DMSTAG_MAX_DIM; ++i) {
     33:     PetscFree(stag->l[i]);
     34:   }
     35:   VecScatterDestroy(&stag->gtol);
     36:   VecScatterDestroy(&stag->ltog_injective);
     37:   PetscFree(stag->neighbors);
     38:   PetscFree(stag->locationOffsets);
     39:   PetscFree(stag->coordinateDMType);
     40:   PetscFree(stag);
     41:   return(0);
     42: }

     44: static PetscErrorCode DMCreateGlobalVector_Stag(DM dm,Vec *vec)
     45: {
     46:   PetscErrorCode  ierr;
     47:   DM_Stag * const stag = (DM_Stag*)dm->data;

     50:   VecCreateMPI(PetscObjectComm((PetscObject)dm),stag->entries,PETSC_DECIDE,vec);
     51:   VecSetDM(*vec,dm);
     52:   /* Could set some ops, as DMDA does */
     53:   VecSetLocalToGlobalMapping(*vec,dm->ltogmap);
     54:   return(0);
     55: }

     57: static PetscErrorCode DMCreateLocalVector_Stag(DM dm,Vec *vec)
     58: {
     59:   PetscErrorCode  ierr;
     60:   DM_Stag * const stag = (DM_Stag*)dm->data;

     63:   VecCreateSeq(PETSC_COMM_SELF,stag->entriesGhost,vec);
     64:   VecSetBlockSize(*vec,stag->entriesPerElement);
     65:   VecSetDM(*vec,dm);
     66:   return(0);
     67: }

     69: static PetscErrorCode DMCreateMatrix_Stag(DM dm,Mat *mat)
     70: {
     71:   PetscErrorCode         ierr;
     72:   MatType                matType;
     73:   PetscBool              isaij,isshell;
     74:   PetscInt               entries,width,nNeighbors,dim,dof[DMSTAG_MAX_STRATA],stencilWidth;
     75:   DMStagStencilType      stencilType;
     76:   ISLocalToGlobalMapping ltogmap;

     79:   DMGetDimension(dm,&dim);
     80:   DMGetMatType(dm,&matType);
     81:   PetscStrcmp(matType,MATAIJ,&isaij);
     82:   PetscStrcmp(matType,MATSHELL,&isshell);
     83:   DMStagGetEntries(dm,&entries);
     84:   DMStagGetDOF(dm,&dof[0],&dof[1],&dof[2],&dof[3]);
     85:   DMStagGetStencilWidth(dm,&stencilWidth);
     86:   DMStagGetStencilType(dm,&stencilType);

     88:   if (isaij) {
     89:     /* This implementation gives a very dense stencil, which is likely unsuitable for
     90:        real applications. */
     91:     switch (stencilType) {
     92:       case DMSTAG_STENCIL_NONE:
     93:         nNeighbors = 1;
     94:         break;
     95:       case DMSTAG_STENCIL_STAR:
     96:         switch (dim) {
     97:           case 1 :
     98:             nNeighbors = 2*stencilWidth + 1;
     99:             break;
    100:           case 2 :
    101:             nNeighbors = 4*stencilWidth + 3;
    102:             break;
    103:           case 3 :
    104:             nNeighbors = 6*stencilWidth + 5;
    105:             break;
    106:           default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
    107:         }
    108:         break;
    109:       case DMSTAG_STENCIL_BOX:
    110:         switch (dim) {
    111:           case 1 :
    112:             nNeighbors = (2*stencilWidth + 1);
    113:             break;
    114:           case 2 :
    115:             nNeighbors = (2*stencilWidth + 1) * (2*stencilWidth + 1);
    116:             break;
    117:           case 3 :
    118:             nNeighbors = (2*stencilWidth + 1) * (2*stencilWidth + 1) * (2*stencilWidth + 1);
    119:             break;
    120:           default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
    121:         }
    122:         break;
    123:       default : SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported stencil");
    124:     }
    125:     width = (dof[0] + dof[1] + dof[2] + dof[3]) * nNeighbors;
    126:     MatCreateAIJ(PetscObjectComm((PetscObject)dm),entries,entries,PETSC_DETERMINE,PETSC_DETERMINE,width,NULL,width,NULL,mat);
    127:   } else if (isshell) {
    128:     MatCreate(PetscObjectComm((PetscObject)dm),mat);
    129:     MatSetSizes(*mat,entries,entries,PETSC_DETERMINE,PETSC_DETERMINE);
    130:     MatSetType(*mat,MATSHELL);
    131:     MatSetUp(*mat);
    132:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for Mattype %s",matType);

    134:   DMGetLocalToGlobalMapping(dm,&ltogmap);
    135:   MatSetLocalToGlobalMapping(*mat,ltogmap,ltogmap);
    136:   MatSetDM(*mat,dm);
    137:   return(0);
    138: }

    140: static PetscErrorCode DMGetCompatibility_Stag(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
    141: {
    142:   PetscErrorCode  ierr;
    143:   const DM_Stag * const stag  = (DM_Stag*)dm->data;
    144:   const DM_Stag * const stag2 = (DM_Stag*)dm2->data;
    145:   PetscInt              dim,dim2,i;
    146:   MPI_Comm              comm;
    147:   PetscMPIInt           sameComm;
    148:   DMType                type2;
    149:   PetscBool             sameType;

    152:   DMGetType(dm2,&type2);
    153:   PetscStrcmp(DMSTAG,type2,&sameType);
    154:   if (!sameType) {
    155:     PetscInfo1((PetscObject)dm,"DMStag compatibility check not implemented with DM of type %s\n",type2);
    156:     *set = PETSC_FALSE;
    157:     return(0);
    158:   }

    160:   PetscObjectGetComm((PetscObject)dm,&comm);
    161:   MPI_Comm_compare(comm,PetscObjectComm((PetscObject)dm2),&sameComm);
    162:   if (sameComm != MPI_IDENT) {
    163:     PetscInfo2((PetscObject)dm,"DMStag objects have different communicators: %d != %d\n",comm,PetscObjectComm((PetscObject)dm2));
    164:     *set = PETSC_FALSE;
    165:     return(0);
    166:   }
    167:   DMGetDimension(dm ,&dim);
    168:   DMGetDimension(dm2,&dim2);
    169:   if (dim != dim2) {
    170:     PetscInfo((PetscObject)dm,"DMStag objects have different dimensions");
    171:     *set = PETSC_TRUE;
    172:     *compatible = PETSC_FALSE;
    173:     return(0);
    174:   }
    175:   for (i=0; i<dim; ++i) {
    176:     if (stag->N[i] != stag2->N[i]) {
    177:       PetscInfo3((PetscObject)dm,"DMStag objects have different global numbers of elements in dimension %D: %D != %D\n",i,stag->n[i],stag2->n[i]);
    178:       *set = PETSC_TRUE;
    179:       *compatible = PETSC_FALSE;
    180:       return(0);
    181:     }
    182:     if (stag->n[i] != stag2->n[i]) {
    183:       PetscInfo3((PetscObject)dm,"DMStag objects have different local numbers of elements in dimension %D: %D != %D\n",i,stag->n[i],stag2->n[i]);
    184:       *set = PETSC_TRUE;
    185:       *compatible = PETSC_FALSE;
    186:       return(0);
    187:     }
    188:     if (stag->boundaryType[i] != stag2->boundaryType[i]) {
    189:       PetscInfo3((PetscObject)dm,"DMStag objects have different boundary types in dimension %d: %s != %s\n",i,stag->boundaryType[i],stag2->boundaryType[i]);
    190:       *set = PETSC_TRUE;
    191:       *compatible = PETSC_FALSE;
    192:       return(0);
    193:     }
    194:   }
    195:   /* Note: we include stencil type and width in the notion of compatibility, as this affects
    196:      the "atlas" (local subdomains). This might be irritating in legitimate cases
    197:      of wanting to transfer between two other-wise compatible DMs with different
    198:      stencil characteristics. */
    199:   if (stag->stencilType != stag2->stencilType) {
    200:     PetscInfo2((PetscObject)dm,"DMStag objects have different ghost stencil types: %s != %s\n",DMStagStencilTypes[stag->stencilType],DMStagStencilTypes[stag2->stencilType]);
    201:     *set = PETSC_TRUE;
    202:     *compatible = PETSC_FALSE;
    203:     return(0);
    204:   }
    205:   if (stag->stencilWidth != stag2->stencilWidth) {
    206:     PetscInfo2((PetscObject)dm,"DMStag objects have different ghost stencil widths: %D != %D\n",stag->stencilWidth,stag->stencilWidth);
    207:     *set = PETSC_TRUE;
    208:     *compatible = PETSC_FALSE;
    209:     return(0);
    210:   }
    211:   *set = PETSC_TRUE;
    212:   *compatible = PETSC_TRUE;
    213:   return(0);
    214: }

    217: /*
    218: Note there are several orderings in play here.
    219: In all cases, non-element dof are associated with the element that they are below/left/behind, and the order in 2D proceeds vertex/bottom edge/left edge/element (with all dof on each together).
    220: Also in all cases, only subdomains which are the last in their dimension have partial elements.

    222: 1) "Natural" Ordering (not used). Number adding each full or partial (on the right or top) element, starting at the bottom left (i=0,j=0) and proceeding across the entire domain, row by row to get a global numbering.
    223: 2) Global ("PETSc") ordering. The same as natural, but restricted to each domain. So, traverse all elements (again starting at the bottom left and going row-by-row) on rank 0, then continue numbering with rank 1, and so on.
    224: 3) Local ordering. Including ghost elements (both interior and on the right/top/front to complete partial elements), use the same convention to create a local numbering.
    225: */

    227: static PetscErrorCode DMLocalToGlobalBegin_Stag(DM dm,Vec l,InsertMode mode,Vec g)
    228: {
    229:   PetscErrorCode  ierr;
    230:   DM_Stag * const stag = (DM_Stag*)dm->data;

    233:   if (mode == ADD_VALUES) {
    234:     VecScatterBegin(stag->gtol,l,g,mode,SCATTER_REVERSE);
    235:   } else if (mode == INSERT_VALUES) {
    236:     if (stag->ltog_injective) {
    237:       VecScatterBegin(stag->ltog_injective,l,g,mode,SCATTER_FORWARD);
    238:     } else {
    239:       VecScatterBegin(stag->gtol,l,g,mode,SCATTER_REVERSE_LOCAL);
    240:     }
    241:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported InsertMode");
    242:   return(0);
    243: }

    245: static PetscErrorCode DMLocalToGlobalEnd_Stag(DM dm,Vec l,InsertMode mode,Vec g)
    246: {
    247:   PetscErrorCode  ierr;
    248:   DM_Stag * const stag = (DM_Stag*)dm->data;

    251:   if (mode == ADD_VALUES) {
    252:     VecScatterEnd(stag->gtol,l,g,mode,SCATTER_REVERSE);
    253:   } else if (mode == INSERT_VALUES) {
    254:     if (stag->ltog_injective) {
    255:       VecScatterEnd(stag->ltog_injective,l,g,mode,SCATTER_FORWARD);
    256:     } else {
    257:       VecScatterEnd(stag->gtol,l,g,mode,SCATTER_REVERSE_LOCAL);
    258:     }
    259:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported InsertMode");
    260:   return(0);
    261: }

    263: static PetscErrorCode DMGlobalToLocalBegin_Stag(DM dm,Vec g,InsertMode mode,Vec l)
    264: {
    265:   PetscErrorCode  ierr;
    266:   DM_Stag * const stag = (DM_Stag*)dm->data;

    269:   VecScatterBegin(stag->gtol,g,l,mode,SCATTER_FORWARD);
    270:   return(0);
    271: }

    273: static PetscErrorCode DMGlobalToLocalEnd_Stag(DM dm,Vec g,InsertMode mode,Vec l)
    274: {
    275:   PetscErrorCode  ierr;
    276:   DM_Stag * const stag = (DM_Stag*)dm->data;

    279:   VecScatterEnd(stag->gtol,g,l,mode,SCATTER_FORWARD);
    280:   return(0);
    281: }

    283: /*
    284: If a stratum is active (non-zero dof), make it active in the coordinate DM.
    285: */
    286: static PetscErrorCode DMCreateCoordinateDM_Stag(DM dm,DM *dmc)
    287: {
    288:   PetscErrorCode  ierr;
    289:   DM_Stag * const stag = (DM_Stag*)dm->data;
    290:   PetscInt        dim;
    291:   PetscBool       isstag,isproduct;

    295:   if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Before creating a coordinate DM, a type must be specified with DMStagSetCoordinateDMType()");

    297:   DMGetDimension(dm,&dim);
    298:   PetscStrcmp(stag->coordinateDMType,DMSTAG,&isstag);
    299:   PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&isproduct);
    300:   if (isstag) {
    301:     DMStagCreateCompatibleDMStag(dm,
    302:         stag->dof[0] > 0 ? dim : 0,
    303:         stag->dof[1] > 0 ? dim : 0,
    304:         stag->dof[2] > 0 ? dim : 0,
    305:         stag->dof[3] > 0 ? dim : 0,
    306:         dmc);
    307:   } else if (isproduct) {
    308:     DMCreate(PETSC_COMM_WORLD,dmc);
    309:     DMSetType(*dmc,DMPRODUCT);
    310:     DMSetDimension(*dmc,dim);
    311:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported coordinate DM type %s",stag->coordinateDMType);
    312:   return(0);
    313: }

    315: static PetscErrorCode DMGetNeighbors_Stag(DM dm,PetscInt *nRanks,const PetscMPIInt *ranks[])
    316: {
    317:   PetscErrorCode  ierr;
    318:   DM_Stag * const stag = (DM_Stag*)dm->data;
    319:   PetscInt        dim;

    322:   DMGetDimension(dm,&dim);
    323:   switch (dim) {
    324:     case 1: *nRanks = 3; break;
    325:     case 2: *nRanks = 9; break;
    326:     case 3: *nRanks = 27; break;
    327:     default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Get neighbors not implemented for dim = %D",dim);
    328:   }
    329:   *ranks = stag->neighbors;
    330:   return(0);
    331: }

    333: static PetscErrorCode DMView_Stag(DM dm,PetscViewer viewer)
    334: {
    335:   PetscErrorCode  ierr;
    336:   DM_Stag * const stag = (DM_Stag*)dm->data;
    337:   PetscBool       isascii,viewAllRanks;
    338:   PetscMPIInt     rank,size;
    339:   PetscInt        dim,maxRanksToView,i;

    342:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
    343:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
    344:   DMGetDimension(dm,&dim);
    345:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
    346:   if (isascii) {
    347:     PetscViewerASCIIPrintf(viewer,"Dimension: %D\n",dim);
    348:     switch (dim) {
    349:       case 1:
    350:         PetscViewerASCIIPrintf(viewer,"Global size: %D\n",stag->N[0]);
    351:         break;
    352:       case 2:
    353:         PetscViewerASCIIPrintf(viewer,"Global sizes: %D x %D\n",stag->N[0],stag->N[1]);
    354:         PetscViewerASCIIPrintf(viewer,"Parallel decomposition: %D x %D ranks\n",stag->nRanks[0],stag->nRanks[1]);
    355:         break;
    356:       case 3:
    357:         PetscViewerASCIIPrintf(viewer,"Global sizes: %D x %D x %D\n",stag->N[0],stag->N[1],stag->N[2]);
    358:         PetscViewerASCIIPrintf(viewer,"Parallel decomposition: %D x %D x %D ranks\n",stag->nRanks[0],stag->nRanks[1],stag->nRanks[2]);
    359:         break;
    360:       default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"not implemented for dim==%D",dim);
    361:     }
    362:     PetscViewerASCIIPrintf(viewer,"Boundary ghosting:");
    363:     for (i=0; i<dim; ++i) {
    364:       PetscViewerASCIIPrintf(viewer," %s",DMBoundaryTypes[stag->boundaryType[i]]);
    365:     }
    366:     PetscViewerASCIIPrintf(viewer,"\n");
    367:     PetscViewerASCIIPrintf(viewer,"Elementwise ghost stencil: %s, width %D\n",DMStagStencilTypes[stag->stencilType],stag->stencilWidth);
    368:     PetscViewerASCIIPrintf(viewer,"Stratum dof:");
    369:     for (i=0; i<dim+1; ++i) {
    370:       PetscViewerASCIIPrintf(viewer," %D:%D",i,stag->dof[i]);
    371:     }
    372:     PetscViewerASCIIPrintf(viewer,"\n");
    373:     if (dm->coordinateDM) {
    374:       PetscViewerASCIIPrintf(viewer,"Has coordinate DM\n");
    375:     }
    376:     maxRanksToView = 16;
    377:     viewAllRanks = (PetscBool)(size <= maxRanksToView);
    378:     if (viewAllRanks) {
    379:       PetscViewerASCIIPushSynchronized(viewer);
    380:       switch (dim) {
    381:         case 1:
    382:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local elements : %D (%D with ghosts)\n",rank,stag->n[0],stag->nGhost[0]);
    383:           break;
    384:         case 2:
    385:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Rank coordinates (%d,%d)\n",rank,stag->rank[0],stag->rank[1]);
    386:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local elements : %D x %D (%D x %D with ghosts)\n",rank,stag->n[0],stag->n[1],stag->nGhost[0],stag->nGhost[1]);
    387:           break;
    388:         case 3:
    389:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Rank coordinates (%d,%d,%d)\n",rank,stag->rank[0],stag->rank[1],stag->rank[2]);
    390:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local elements : %D x %D x %D (%D x %D x %D with ghosts)\n",rank,stag->n[0],stag->n[1],stag->n[2],stag->nGhost[0],stag->nGhost[1],stag->nGhost[2]);
    391:           break;
    392:         default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"not implemented for dim==%D",dim);
    393:       }
    394:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local native entries: %d\n",rank,stag->entries);
    395:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local entries total : %d\n",rank,stag->entriesGhost);
    396:       PetscViewerFlush(viewer);
    397:       PetscViewerASCIIPopSynchronized(viewer);
    398:     } else {
    399:       PetscViewerASCIIPrintf(viewer,"(Per-rank information omitted since >%D ranks used)\n",maxRanksToView);
    400:     }
    401:   }
    402:   return(0);
    403: }

    405: static PetscErrorCode DMSetFromOptions_Stag(PetscOptionItems *PetscOptionsObject,DM dm)
    406: {
    407:   PetscErrorCode  ierr;
    408:   DM_Stag * const stag = (DM_Stag*)dm->data;
    409:   PetscInt        dim;

    412:   DMGetDimension(dm,&dim);
    413:   PetscOptionsHead(PetscOptionsObject,"DMStag Options");
    414:   PetscOptionsInt("-stag_grid_x","Number of grid points in x direction","DMStagSetGlobalSizes",stag->N[0],&stag->N[0],NULL);
    415:   if (dim > 1) { PetscOptionsInt("-stag_grid_y","Number of grid points in y direction","DMStagSetGlobalSizes",stag->N[1],&stag->N[1],NULL); }
    416:   if (dim > 2) { PetscOptionsInt("-stag_grid_z","Number of grid points in z direction","DMStagSetGlobalSizes",stag->N[2],&stag->N[2],NULL); }
    417:   PetscOptionsInt("-stag_ranks_x","Number of ranks in x direction","DMStagSetNumRanks",stag->nRanks[0],&stag->nRanks[0],NULL);
    418:   if (dim > 1) { PetscOptionsInt("-stag_ranks_y","Number of ranks in y direction","DMStagSetNumRanks",stag->nRanks[1],&stag->nRanks[1],NULL); }
    419:   if (dim > 2) { PetscOptionsInt("-stag_ranks_z","Number of ranks in z direction","DMStagSetNumRanks",stag->nRanks[2],&stag->nRanks[2],NULL); }
    420:   PetscOptionsInt("-stag_stencil_width","Elementwise stencil width","DMStagSetStencilWidth",stag->stencilWidth,&stag->stencilWidth,NULL);
    421:   PetscOptionsEnum("-stag_stencil_type","Elementwise stencil stype","DMStagSetStencilType",DMStagStencilTypes,(PetscEnum)stag->stencilType,(PetscEnum*)&stag->stencilType,NULL);
    422:   PetscOptionsEnum("-stag_boundary_type_x","Treatment of (physical) boundaries in x direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[0],(PetscEnum*)&stag->boundaryType[0],NULL);
    423:   PetscOptionsEnum("-stag_boundary_type_y","Treatment of (physical) boundaries in y direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[1],(PetscEnum*)&stag->boundaryType[1],NULL);
    424:   PetscOptionsEnum("-stag_boundary_type_z","Treatment of (physical) boundaries in z direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[2],(PetscEnum*)&stag->boundaryType[2],NULL);
    425:   PetscOptionsInt("-stag_dof_0","Number of dof per 0-cell (vertex/corner)","DMStagSetDOF",stag->dof[0],&stag->dof[0],NULL);
    426:   PetscOptionsInt("-stag_dof_1","Number of dof per 1-cell (edge)",         "DMStagSetDOF",stag->dof[1],&stag->dof[1],NULL);
    427:   PetscOptionsInt("-stag_dof_2","Number of dof per 2-cell (face)",         "DMStagSetDOF",stag->dof[2],&stag->dof[2],NULL);
    428:   PetscOptionsInt("-stag_dof_3","Number of dof per 3-cell (hexahedron)",   "DMStagSetDOF",stag->dof[3],&stag->dof[3],NULL);
    429:   PetscOptionsTail();
    430:   return(0);
    431: }

    433: /*MC
    434:   DMSTAG = "stag" - A DM object representing a "staggered grid" or a structured cell complex.

    436:   This implementation parallels the DMDA implementation in many ways, but allows degrees of freedom
    437:   to be associated with all "strata" in a logically-rectangular grid: vertices, edges, faces, and elements.

    439:   Level: beginner

    441: .seealso: DM, DMPRODUCT, DMDA, DMPLEX, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMType, DMCreate(), DMSetType()
    442: M*/

    444: PETSC_EXTERN PetscErrorCode DMCreate_Stag(DM dm)
    445: {
    447:   DM_Stag        *stag;
    448:   PetscInt       i,dim;

    452:   PetscNewLog(dm,&stag);
    453:   dm->data = stag;

    455:   stag->gtol                                          = NULL;
    456:   stag->ltog_injective                                = NULL;
    457:   for (i=0; i<DMSTAG_MAX_STRATA; ++i) stag->dof[i]    = 0;
    458:   for (i=0; i<DMSTAG_MAX_DIM;    ++i) stag->l[i]      = NULL;
    459:   stag->stencilType                                   = DMSTAG_STENCIL_NONE;
    460:   stag->stencilWidth                                  = 0;
    461:   for (i=0; i<DMSTAG_MAX_DIM;    ++i) stag->nRanks[i] = -1;
    462:   stag->coordinateDMType                              = NULL;

    464:   DMGetDimension(dm,&dim);
    465:   if (dim != 1 && dim != 2 && dim != 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DMSetDimension() must be called to set a dimension with value 1, 2, or 3");

    467:   PetscMemzero(dm->ops,sizeof(*(dm->ops)));
    468:   dm->ops->createcoordinatedm  = DMCreateCoordinateDM_Stag;
    469:   dm->ops->createglobalvector  = DMCreateGlobalVector_Stag;
    470:   dm->ops->createinterpolation = NULL;
    471:   dm->ops->createlocalvector   = DMCreateLocalVector_Stag;
    472:   dm->ops->creatematrix        = DMCreateMatrix_Stag;
    473:   dm->ops->destroy             = DMDestroy_Stag;
    474:   dm->ops->getneighbors        = DMGetNeighbors_Stag;
    475:   dm->ops->globaltolocalbegin  = DMGlobalToLocalBegin_Stag;
    476:   dm->ops->globaltolocalend    = DMGlobalToLocalEnd_Stag;
    477:   dm->ops->localtoglobalbegin  = DMLocalToGlobalBegin_Stag;
    478:   dm->ops->localtoglobalend    = DMLocalToGlobalEnd_Stag;
    479:   dm->ops->setfromoptions      = DMSetFromOptions_Stag;
    480:   switch (dim) {
    481:     case 1: dm->ops->setup     = DMSetUp_Stag_1d; break;
    482:     case 2: dm->ops->setup     = DMSetUp_Stag_2d; break;
    483:     case 3: dm->ops->setup     = DMSetUp_Stag_3d; break;
    484:     default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
    485:   }
    486:   dm->ops->clone               = DMClone_Stag;
    487:   dm->ops->view                = DMView_Stag;
    488:   dm->ops->getcompatibility    = DMGetCompatibility_Stag;
    489:   return(0);
    490: }