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:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 37:   PetscFunctionBeginUser;
 38:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 39:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 40:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

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

 49:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 50:      Create DMDA object.
 51:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
 52:   if (dim == 3) {
 53:     PetscCall(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:     PetscCall(DMSetFromOptions(da));
 55:     PetscCall(DMSetUp(da));
 56:   } else {
 57:     PetscCall(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:     PetscCall(DMSetFromOptions(da));
 59:     PetscCall(DMSetUp(da));
 60:   }

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

 68:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 69:      Populate the 3D vector
 70:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 71:   PetscCall(DMDAGetCorners(da, &ixs, &iys, &izs, &ixm, &iym, &izm));
 72:   if (dim == 3) {
 73:     PetscCall(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:     PetscCall(DMDAVecRestoreArray(da, vec_full, &vecdata3d));
 80:   } else {
 81:     PetscCall(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:     PetscCall(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:   PetscCall(DMDACreatePatchIS(da, &lower, &upper, &selectis, patchis_offproc));
120:   PetscCall(ISView(selectis, PETSC_VIEWER_STDOUT_WORLD));

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

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

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

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

143:   PetscCall(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*/