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"; 6 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 */ 12 13 #include <petscdmda.h> 14 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 */ 33 PetscErrorCode ierr; /* error checking */ 34 35 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 36 Initialize program and set problem parameters 37 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 38 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 39 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 40 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 41 42 ierr = PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL);CHKERRQ(ierr); 43 ierr = PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL);CHKERRQ(ierr); 44 ierr = PetscOptionsGetInt(NULL, NULL, "-mz", &mz, NULL);CHKERRQ(ierr); 45 ierr = PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL);CHKERRQ(ierr); 46 ierr = PetscOptionsGetInt(NULL, NULL, "-sliceid", &sliceid, NULL);CHKERRQ(ierr); 47 ierr = PetscOptionsGetInt(NULL, NULL, "-sliceaxis", &sliceaxis, NULL);CHKERRQ(ierr); 48 49 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 50 Create DMDA object. 51 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 52 if (dim==3) { 53 ierr = 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);CHKERRQ(ierr); 61 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 62 ierr = DMSetUp(da);CHKERRQ(ierr); 63 } else { 64 ierr = 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);CHKERRQ(ierr); 72 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 73 ierr = DMSetUp(da);CHKERRQ(ierr); 74 } 75 76 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 77 Create the parent vector 78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 79 ierr = DMCreateGlobalVector(da, &vec_full);CHKERRQ(ierr); 80 ierr = PetscObjectSetName((PetscObject) vec_full, "full_vector");CHKERRQ(ierr); 81 82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 83 Populate the 3D vector 84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 85 ierr = DMDAGetCorners(da, &ixs, &iys, &izs, &ixm, &iym, &izm);CHKERRQ(ierr); 86 if (dim==3) { 87 ierr = DMDAVecGetArray(da, vec_full, &vecdata3d);CHKERRQ(ierr); 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 ierr = DMDAVecRestoreArray(da, vec_full, &vecdata3d);CHKERRQ(ierr); 96 } else { 97 ierr = DMDAVecGetArray(da, vec_full, &vecdata2d);CHKERRQ(ierr); 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 ierr = DMDAVecRestoreArray(da, vec_full, &vecdata2d);CHKERRQ(ierr); 104 } 105 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 ierr = DMDACreatePatchIS(da, &lower, &upper, &selectis, patchis_offproc);CHKERRQ(ierr); 120 ierr = ISView(selectis, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 121 122 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 123 Use the obtained IS to extract the slice as a subvector 124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 125 ierr = VecGetSubVector(vec_full, selectis, &vec_slice);CHKERRQ(ierr); 126 127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 128 View the extracted subvector 129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 130 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_DENSE);CHKERRQ(ierr); 131 ierr = VecView(vec_slice, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 132 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 133 134 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 135 Restore subvector, destroy data structures and exit. 136 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 137 ierr = VecRestoreSubVector(vec_full, selectis, &vec_slice);CHKERRQ(ierr); 138 139 ierr = ISDestroy(&selectis);CHKERRQ(ierr); 140 ierr = DMDestroy(&da);CHKERRQ(ierr); 141 ierr = VecDestroy(&vec_full);CHKERRQ(ierr); 142 143 ierr = PetscFinalize(); 144 return ierr; 145 } 146 147 /*TEST 148 149 test: 150 nsize: 1 151 args: -sliceaxis 0 152 153 test: 154 suffix: 2 155 nsize: 2 156 args: -sliceaxis 1 157 158 test: 159 suffix: 3 160 nsize: 4 161 args: -sliceaxis 2 162 163 test: 164 suffix: 4 165 nsize: 2 166 args: -sliceaxis 1 -dim 2 167 168 test: 169 suffix: 5 170 nsize: 3 171 args: -sliceaxis 0 -dim 2 172 173 TEST*/ 174