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 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 39 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 40 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 41 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)); 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);PetscCall(ierr); 61 PetscCall(DMSetFromOptions(da)); 62 PetscCall(DMSetUp(da)); 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);PetscCall(ierr); 72 PetscCall(DMSetFromOptions(da)); 73 PetscCall(DMSetUp(da)); 74 } 75 76 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 77 Create the parent vector 78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 79 PetscCall(DMCreateGlobalVector(da, &vec_full)); 80 PetscCall(PetscObjectSetName((PetscObject) vec_full, "full_vector")); 81 82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 83 Populate the 3D vector 84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 85 PetscCall(DMDAGetCorners(da, &ixs, &iys, &izs, &ixm, &iym, &izm)); 86 if (dim==3) { 87 PetscCall(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 PetscCall(DMDAVecRestoreArray(da, vec_full, &vecdata3d)); 96 } else { 97 PetscCall(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 PetscCall(DMDAVecRestoreArray(da, vec_full, &vecdata2d)); 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 PetscCall(DMDACreatePatchIS(da, &lower, &upper, &selectis, patchis_offproc)); 120 PetscCall(ISView(selectis, PETSC_VIEWER_STDOUT_WORLD)); 121 122 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 123 Use the obtained IS to extract the slice as a subvector 124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 125 PetscCall(VecGetSubVector(vec_full, selectis, &vec_slice)); 126 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)); 133 134 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 135 Restore subvector, destroy data structures and exit. 136 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 137 PetscCall(VecRestoreSubVector(vec_full, selectis, &vec_slice)); 138 139 PetscCall(ISDestroy(&selectis)); 140 PetscCall(DMDestroy(&da)); 141 PetscCall(VecDestroy(&vec_full)); 142 143 PetscCall(PetscFinalize()); 144 return 0; 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