2: static char help[] = "\n\n";

  4: /*
  5:      Demonstrates using DM_BOUNDARY_GHOSTED how to handle a rotated boundary conditions where one edge
  6:     is connected to its immediate neighbor

  8:     Consider the domain (with natural numbering)

 10:      6   7   8
 11:      3   4   5
 12:      0   1   2

 14:     The ghost points along the bottom (directly below the three columns above) should be 0 3 and 6
 15:     while the ghost points along the left side should be 0 1 2

 17:     Note that the ghosted local vectors extend in both the x and y directions so, for example if we have a
 18:     single MPI process the ghosted vector has (in the original natural numbering)

 20:      x  x  x  x  x
 21:      2  6  7  8  x
 22:      1  3  4  5  x
 23:      0  0  1  2  x
 24:      x  0  3  6  x

 26:     where x indicates a location that is not updated by the communication and should be used.

 28:     For this to make sense the number of grid points in the x and y directions must be the same

 30:     This ghost point mapping was suggested by: Wenbo Zhao <zhaowenbo.npic@gmail.com>
 31: */

 33: #include <petscdm.h>
 34: #include <petscdmda.h>

 36: int main(int argc,char **argv)
 37: {
 38:   PetscInt         M = 6;
 39:   DM               da;
 40:   Vec              local,global,natural;
 41:   PetscInt         i,start,end,*ifrom,x,y,xm,ym;
 42:   PetscScalar      *xnatural;
 43:   IS               from,to;
 44:   AO               ao;
 45:   VecScatter       scatter1,scatter2;
 46:   PetscViewer      subviewer;

 48:   PetscInitialize(&argc,&argv,(char*)0,help);

 50:   /* Create distributed array and get vectors */
 51:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_GHOSTED,DM_BOUNDARY_GHOSTED,DMDA_STENCIL_STAR,M,M,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
 52:   DMSetUp(da);
 53:   DMCreateGlobalVector(da,&global);
 54:   DMCreateLocalVector(da,&local);

 56:   /* construct global to local scatter for the left side of the domain to the ghost on the bottom */
 57:   DMDAGetCorners(da,&x,&y,NULL,&xm,&ym,NULL);
 58:   if (!y) { /* only processes on the bottom of the domain fill up the ghost locations */
 59:     ISCreateStride(PETSC_COMM_SELF,xm,1,1,&to);
 60:   } else {
 61:     ISCreateStride(PETSC_COMM_SELF,0,0,0,&to);
 62:   }
 63:   PetscMalloc1(xm,&ifrom);
 64:   for (i=x;i<x+xm;i++) {
 65:     ifrom[i-x] = M*i;
 66:   }
 67:   DMDAGetAO(da,&ao);
 68:   AOApplicationToPetsc(ao,xm,ifrom);
 69:   if (!y) {
 70:     ISCreateGeneral(PETSC_COMM_WORLD,xm,ifrom,PETSC_OWN_POINTER,&from);
 71:   } else {
 72:     PetscFree(ifrom);
 73:     ISCreateGeneral(PETSC_COMM_WORLD,0,NULL,PETSC_COPY_VALUES,&from);
 74:   }
 75:   VecScatterCreate(global,from,local,to,&scatter1);
 76:   ISDestroy(&to);
 77:   ISDestroy(&from);

 79:   /* construct global to local scatter for the bottom side of the domain to the ghost on the right */
 80:   if (!x) { /* only processes on the left side of the domain fill up the ghost locations */
 81:     ISCreateStride(PETSC_COMM_SELF,ym,xm+2,xm+2,&to);
 82:   } else {
 83:     ISCreateStride(PETSC_COMM_SELF,0,0,0,&to);
 84:   }
 85:   PetscMalloc1(ym,&ifrom);
 86:   for (i=y;i<y+ym;i++) {
 87:     ifrom[i-y] = i;
 88:   }
 89:   DMDAGetAO(da,&ao);
 90:   AOApplicationToPetsc(ao,ym,ifrom);
 91:   if (!x) {
 92:     ISCreateGeneral(PETSC_COMM_WORLD,ym,ifrom,PETSC_OWN_POINTER,&from);
 93:   } else {
 94:     PetscFree(ifrom);
 95:     ISCreateGeneral(PETSC_COMM_WORLD,0,NULL,PETSC_COPY_VALUES,&from);
 96:   }
 97:   VecScatterCreate(global,from,local,to,&scatter2);
 98:   ISDestroy(&to);
 99:   ISDestroy(&from);

101:   /*
102:      fill the global vector with the natural global numbering for each local entry
103:      this is only done for testing purposes since it is easy to see if the scatter worked correctly
104:   */
105:   DMDACreateNaturalVector(da,&natural);
106:   VecGetOwnershipRange(natural,&start,&end);
107:   VecGetArray(natural,&xnatural);
108:   for (i=start; i<end; i++) {
109:     xnatural[i-start] = i;
110:   }
111:   VecRestoreArray(natural,&xnatural);
112:   DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,global);
113:   DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,global);
114:   VecDestroy(&natural);

116:   /* scatter from global to local */
117:   VecScatterBegin(scatter1,global,local,INSERT_VALUES,SCATTER_FORWARD);
118:   VecScatterEnd(scatter1,global,local,INSERT_VALUES,SCATTER_FORWARD);
119:   VecScatterBegin(scatter2,global,local,INSERT_VALUES,SCATTER_FORWARD);
120:   VecScatterEnd(scatter2,global,local,INSERT_VALUES,SCATTER_FORWARD);
121:   /*
122:      normally here you would also call
123:   DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);
124:   DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);
125:     to update all the interior ghost cells between neighboring processes.
126:     We don't do it here since this is only a test of "special" ghost points.
127:   */

129:   /* view each local ghosted vector */
130:   PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);
131:   VecView(local,subviewer);
132:   PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);

134:   VecScatterDestroy(&scatter1);
135:   VecScatterDestroy(&scatter2);
136:   VecDestroy(&local);
137:   VecDestroy(&global);
138:   DMDestroy(&da);
139:   PetscFinalize();
140:   return 0;
141: }

143: /*TEST

145:    test:

147:    test:
148:       suffix: 2
149:       nsize: 2

151:    test:
152:       suffix: 4
153:       nsize: 4

155:    test:
156:       suffix: 9
157:       nsize: 9

159: TEST*/