Actual source code: ex53.c

  1: static char help[] = "Use DMDACreatePatchIS  to extract a slice from a vector, Command line options :\n\
  2: mx/my/mz - set the dimensions of the parent vector\n\
  3: dim - set the dimensionality of the parent vector (2,3)\n\
  4: sliceaxis - Integer describing the axis along which the sice will be selected (0-X, 1-Y, 2-Z)\n\
  5: sliceid - set the location where the slice will be extracted from the parent vector\n";

  7: /*
  8:    This test checks the functionality of DMDACreatePatchIS when
  9:    extracting a 2D vector from a 3D vector and 1D vector from a
 10:    2D vector.
 11:    */

 13: #include <petscdmda.h>

 15: int main(int argc, char **argv)
 16: {
 17:   PetscMPIInt    rank, size;                    /* MPI rank and size */
 18:   PetscInt       mx = 4, my = 4, mz = 4;        /* Dimensions of parent vector */
 19:   PetscInt       sliceid   = 2;                 /* k (z) index to pick the slice */
 20:   PetscInt       sliceaxis = 2;                 /* Select axis along which the slice will be extracted */
 21:   PetscInt       dim       = 3;                 /* Dimension of the parent vector */
 22:   PetscInt       i, j, k;                       /* Iteration indices */
 23:   PetscInt       ixs, iys, izs;                 /* Corner indices for 3D vector */
 24:   PetscInt       ixm, iym, izm;                 /* Widths of parent vector */
 25:   PetscScalar ***vecdata3d;                     /* Pointer to access 3d parent vector */
 26:   PetscScalar  **vecdata2d;                     /* Pointer to access 2d parent vector */
 27:   DM             da;                            /* 2D/3D DMDA object */
 28:   Vec            vec_full;                      /* Parent vector */
 29:   Vec            vec_slice;                     /* Slice vector */
 30:   MatStencil     lower, upper;                  /* Stencils to select slice */
 31:   IS             selectis;                      /* IS to select slice and extract subvector */
 32:   PetscBool      patchis_offproc = PETSC_FALSE; /* flag to DMDACreatePatchIS indicating that off-proc values are to be ignored */

 34:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 35:      Initialize program and set problem parameters
 36:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 38:   PetscInitialize(&argc, &argv, (char *)0, help);
 39:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 40:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

 42:   PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL);
 43:   PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL);
 44:   PetscOptionsGetInt(NULL, NULL, "-mz", &mz, NULL);
 45:   PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL);
 46:   PetscOptionsGetInt(NULL, NULL, "-sliceid", &sliceid, NULL);
 47:   PetscOptionsGetInt(NULL, NULL, "-sliceaxis", &sliceaxis, NULL);

 49:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 50:      Create DMDA object.
 51:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
 52:   if (dim == 3) {
 53:     DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, mx, my, mz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, NULL, &da);
 54:     DMSetFromOptions(da);
 55:     DMSetUp(da);
 56:   } else {
 57:     DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, mx, my, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da);
 58:     DMSetFromOptions(da);
 59:     DMSetUp(da);
 60:   }

 62:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 63:      Create the parent vector
 64:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
 65:   DMCreateGlobalVector(da, &vec_full);
 66:   PetscObjectSetName((PetscObject)vec_full, "full_vector");

 68:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 69:      Populate the 3D vector
 70:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 71:   DMDAGetCorners(da, &ixs, &iys, &izs, &ixm, &iym, &izm);
 72:   if (dim == 3) {
 73:     DMDAVecGetArray(da, vec_full, &vecdata3d);
 74:     for (k = izs; k < izs + izm; k++) {
 75:       for (j = iys; j < iys + iym; j++) {
 76:         for (i = ixs; i < ixs + ixm; i++) vecdata3d[k][j][i] = ((i - mx / 2) * (j + mx / 2)) + k * 100;
 77:       }
 78:     }
 79:     DMDAVecRestoreArray(da, vec_full, &vecdata3d);
 80:   } else {
 81:     DMDAVecGetArray(da, vec_full, &vecdata2d);
 82:     for (j = iys; j < iys + iym; j++) {
 83:       for (i = ixs; i < ixs + ixm; i++) vecdata2d[j][i] = ((i - mx / 2) * (j + mx / 2));
 84:     }
 85:     DMDAVecRestoreArray(da, vec_full, &vecdata2d);
 86:   }

 88:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 89:      Get an IS corresponding to a 2D slice
 90:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 91:   if (sliceaxis == 0) {
 92:     lower.i = sliceid;
 93:     lower.j = 0;
 94:     lower.k = 0;
 95:     lower.c = 1;
 96:     upper.i = sliceid;
 97:     upper.j = my;
 98:     upper.k = mz;
 99:     upper.c = 1;
100:   } else if (sliceaxis == 1) {
101:     lower.i = 0;
102:     lower.j = sliceid;
103:     lower.k = 0;
104:     lower.c = 1;
105:     upper.i = mx;
106:     upper.j = sliceid;
107:     upper.k = mz;
108:     upper.c = 1;
109:   } else if (sliceaxis == 2 && dim == 3) {
110:     lower.i = 0;
111:     lower.j = 0;
112:     lower.k = sliceid;
113:     lower.c = 1;
114:     upper.i = mx;
115:     upper.j = my;
116:     upper.k = sliceid;
117:     upper.c = 1;
118:   }
119:   DMDACreatePatchIS(da, &lower, &upper, &selectis, patchis_offproc);
120:   ISView(selectis, PETSC_VIEWER_STDOUT_WORLD);

122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123:      Use the obtained IS to extract the slice as a subvector
124:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125:   VecGetSubVector(vec_full, selectis, &vec_slice);

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:      View the extracted subvector
129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_DENSE);
131:   VecView(vec_slice, PETSC_VIEWER_STDOUT_WORLD);
132:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);

134:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135:      Restore subvector, destroy data structures and exit.
136:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137:   VecRestoreSubVector(vec_full, selectis, &vec_slice);

139:   ISDestroy(&selectis);
140:   DMDestroy(&da);
141:   VecDestroy(&vec_full);

143:   PetscFinalize();
144:   return 0;
145: }

147: /*TEST

149:     test:
150:       nsize: 1
151:       args: -sliceaxis 0

153:     test:
154:       suffix: 2
155:       nsize:  2
156:       args: -sliceaxis 1

158:     test:
159:       suffix: 3
160:       nsize:  4
161:       args:  -sliceaxis 2

163:     test:
164:       suffix: 4
165:       nsize:  2
166:       args: -sliceaxis 1 -dim 2

168:     test:
169:       suffix: 5
170:       nsize:  3
171:       args: -sliceaxis 0 -dim 2

173: TEST*/