2: static char help[] = "Demonstrates generating a slice from a DMDA Vector.\n\n";

  4: #include <petscdmda.h>
  5: #include <petscao.h>

  9: /*
 10:     Given a DMDA generates a VecScatter context that will deliver a slice
 11:   of the global vector to each processor. In this example, each processor
 12:   receives the values i=*, j=*, k=rank, i.e. one z plane.

 14:   Note: This code is written assuming only one degree of freedom per node.
 15:   For multiple degrees of freedom per node use ISCreateBlock()
 16:   instead of ISCreateGeneral().
 17: */
 18: PetscErrorCode GenerateSliceScatter(DM da,VecScatter *scatter,Vec *vslice)
 19: {
 20:   AO             ao;
 21:   PetscInt       M,N,P,nslice,*sliceindices,count,i,j;
 22:   PetscMPIInt    rank;
 24:   MPI_Comm       comm;
 25:   Vec            vglobal;
 26:   IS             isfrom,isto;

 28:   PetscObjectGetComm((PetscObject)da,&comm);
 29:   MPI_Comm_rank(comm,&rank);

 31:   DMDAGetAO(da,&ao);
 32:   DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);

 34:   /*
 35:      nslice is number of degrees of freedom in this processors slice
 36:    if there are more processors then z plans the extra processors get 0
 37:    elements in their slice.
 38:   */
 39:   if (rank < P) nslice = M*N;
 40:   else nslice = 0;

 42:   /*
 43:      Generate the local vector to hold this processors slice
 44:   */
 45:   VecCreateSeq(PETSC_COMM_SELF,nslice,vslice);
 46:   DMCreateGlobalVector(da,&vglobal);

 48:   /*
 49:        Generate the indices for the slice in the "natural" global ordering
 50:      Note: this is just an example, one could select any subset of nodes
 51:     on each processor. Just list them in the global natural ordering.

 53:   */
 54:   PetscMalloc1(nslice+1,&sliceindices);
 55:   count = 0;
 56:   if (rank < P) {
 57:     for (j=0; j<N; j++) {
 58:       for (i=0; i<M; i++) {
 59:          sliceindices[count++] = rank*M*N + j*M + i;
 60:       }
 61:     }
 62:   }
 63:   /*
 64:       Convert the indices to the "PETSc" global ordering
 65:   */
 66:   AOApplicationToPetsc(ao,nslice,sliceindices);

 68:   /* Create the "from" and "to" index set */
 69:   /* This is to scatter from the global vector */
 70:   ISCreateGeneral(PETSC_COMM_SELF,nslice,sliceindices,PETSC_OWN_POINTER,&isfrom);
 71:   /* This is to gather into the local vector */
 72:   ISCreateStride(PETSC_COMM_SELF,nslice,0,1,&isto);
 73:   VecScatterCreate(vglobal,isfrom,*vslice,isto,scatter);
 74:   ISDestroy(&isfrom);
 75:   ISDestroy(&isto);
 76:   return 0;
 77: }


 82: int main(int argc,char **argv)
 83: {
 84:   PetscMPIInt      rank;
 85:   PetscInt         m   = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,M = 3,N = 5,P=3,s=1;
 86:   PetscInt         *lx = NULL,*ly = NULL,*lz = NULL;
 87:   PetscErrorCode   ierr;
 88:   PetscBool        flg = PETSC_FALSE;
 89:   DM               da;
 90:   Vec              local,global,vslice;
 91:   PetscScalar      value;
 92:   DMBoundaryType   wrap         = DM_XYPERIODIC;
 93:   DMDAStencilType  stencil_type = DMDA_STENCIL_BOX;
 94:   VecScatter       scatter;

 96:   PetscInitialize(&argc,&argv,(char*)0,help);
 97:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 99:   /* Read options */
100:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
101:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
102:   PetscOptionsGetInt(NULL,NULL,"-P",&P,NULL);
103:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
104:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
105:   PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
106:   PetscOptionsGetInt(NULL,NULL,"-s",&s,NULL);
107:   PetscOptionsGetBool(NULL,NULL,"-star",&flg,NULL);
108:   if (flg) stencil_type =  DMDA_STENCIL_STAR;

110:   /* Create distributed array and get vectors */
111:   DMDACreate3d(PETSC_COMM_WORLD,wrap,stencil_type,M,N,P,m,n,p,1,s,
112:                       lx,ly,lz,&da);
113:   DMView(da,PETSC_VIEWER_DRAW_WORLD);
114:   DMCreateGlobalVector(da,&global);
115:   DMCreateLocalVector(da,&local);

117:   GenerateSliceScatter(da,&scatter,&vslice);

119:   /* Put the value rank+1 into all locations of vslice and transfer back to global vector */
120:   value = 1.0 + rank;
121:   VecSet(vslice,value);
122:   VecScatterBegin(scatter,vslice,global,INSERT_VALUES,SCATTER_REVERSE);
123:   VecScatterEnd(scatter,vslice,global,INSERT_VALUES,SCATTER_REVERSE);

125:   VecView(global,PETSC_VIEWER_DRAW_WORLD);

127:   VecDestroy(&local);
128:   VecDestroy(&global);
129:   DMDestroy(&da);
130:   PetscFinalize();
131:   return 0;
132: }