Actual source code: ex53.c

petsc-3.15.0 2021-03-30
Report Typos and Errors
  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 extraced 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 */

 35:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 36:      Initialize program and set problem parameters
 37:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 38:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 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,
 54:                         DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,
 55:                         DMDA_STENCIL_STAR,
 56:                         mx, my, mz,
 57:                         PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
 58:                         1, 1,
 59:                         NULL, NULL, NULL,
 60:                         &da);
 61:     DMSetFromOptions(da);
 62:     DMSetUp(da);
 63:   } else {
 64:     DMDACreate2d(PETSC_COMM_WORLD,
 65:                         DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,
 66:                         DMDA_STENCIL_STAR,
 67:                         mx, my,
 68:                         PETSC_DECIDE, PETSC_DECIDE,
 69:                         1, 1,
 70:                         NULL, NULL,
 71:                         &da);
 72:     DMSetFromOptions(da);
 73:     DMSetUp(da);
 74:   }

 76:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 77:      Create the parent vector
 78:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
 79:   DMCreateGlobalVector(da, &vec_full);
 80:   PetscObjectSetName((PetscObject) vec_full, "full_vector");

 82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 83:      Populate the 3D vector
 84:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 85:   DMDAGetCorners(da, &ixs, &iys, &izs, &ixm, &iym, &izm);
 86:   if (dim==3) {
 87:     DMDAVecGetArray(da, vec_full, &vecdata3d);
 88:     for (k=izs; k<izs+izm; k++){
 89:       for (j=iys; j<iys+iym; j++) {
 90:         for (i=ixs; i<ixs+ixm; i++) {
 91:           vecdata3d[k][j][i] = ((i-mx/2)*(j+mx/2))+k*100;
 92:         }
 93:       }
 94:     }
 95:     DMDAVecRestoreArray(da, vec_full, &vecdata3d);
 96:   } else {
 97:     DMDAVecGetArray(da, vec_full, &vecdata2d);
 98:     for (j=iys; j<iys+iym; j++) {
 99:       for (i=ixs; i<ixs+ixm; i++) {
100:         vecdata2d[j][i] = ((i-mx/2)*(j+mx/2));
101:       }
102:     }
103:     DMDAVecRestoreArray(da, vec_full, &vecdata2d);
104:   }

106:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107:      Get an IS corresponding to a 2D slice
108:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109:   if (sliceaxis==0) {
110:     lower.i = sliceid; lower.j = 0;  lower.k = 0;  lower.c = 1;
111:     upper.i = sliceid; upper.j = my; upper.k = mz; upper.c = 1;
112:   } else if (sliceaxis==1) {
113:     lower.i = 0;  lower.j = sliceid; lower.k = 0;  lower.c = 1;
114:     upper.i = mx; upper.j = sliceid; upper.k = mz; upper.c = 1;
115:   } else if (sliceaxis==2 && dim==3) {
116:     lower.i = 0;  lower.j = 0;  lower.k = sliceid; lower.c = 1;
117:     upper.i = mx; upper.j = my; upper.k = sliceid; 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 ierr;
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*/