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:   PetscErrorCode   ierr;
 40:   DM               da;
 41:   Vec              local,global,natural;
 42:   PetscInt         i,start,end,*ifrom,x,y,xm,ym;
 43:   PetscScalar      *xnatural;
 44:   IS               from,to;
 45:   AO               ao;
 46:   VecScatter       scatter1,scatter2;
 47:   PetscViewer      subviewer;

 49:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

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

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

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

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

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

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

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



146: /*TEST

148:    test:

150:    test:
151:       suffix: 2
152:       nsize: 2

154:    test:
155:       suffix: 4
156:       nsize: 4

158:    test:
159:       suffix: 9
160:       nsize: 9

162: TEST*/