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;

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

 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++) ifrom[i - x] = M * i;
 66:   DMDAGetAO(da, &ao);
 67:   AOApplicationToPetsc(ao, xm, ifrom);
 68:   if (!y) {
 69:     ISCreateGeneral(PETSC_COMM_WORLD, xm, ifrom, PETSC_OWN_POINTER, &from);
 70:   } else {
 71:     PetscFree(ifrom);
 72:     ISCreateGeneral(PETSC_COMM_WORLD, 0, NULL, PETSC_COPY_VALUES, &from);
 73:   }
 74:   VecScatterCreate(global, from, local, to, &scatter1);
 75:   ISDestroy(&to);
 76:   ISDestroy(&from);

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

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

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

124:   /* view each local ghosted vector */
125:   PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &subviewer);
126:   VecView(local, subviewer);
127:   PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &subviewer);

129:   VecScatterDestroy(&scatter1);
130:   VecScatterDestroy(&scatter2);
131:   VecDestroy(&local);
132:   VecDestroy(&global);
133:   DMDestroy(&da);
134:   PetscFinalize();
135:   return 0;
136: }

138: /*TEST

140:    test:

142:    test:
143:       suffix: 2
144:       nsize: 2

146:    test:
147:       suffix: 4
148:       nsize: 4

150:    test:
151:       suffix: 9
152:       nsize: 9

154: TEST*/