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 34 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 35 Initialize program and set problem parameters 36 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 37 PetscFunctionBeginUser; 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 PetscCall(DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,mx, my, mz, 54 PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,1, 1,NULL, NULL, NULL,&da)); 55 PetscCall(DMSetFromOptions(da)); 56 PetscCall(DMSetUp(da)); 57 } else { 58 PetscCall(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,mx, my,PETSC_DECIDE, PETSC_DECIDE, 59 1, 1,NULL, NULL,&da)); 60 PetscCall(DMSetFromOptions(da)); 61 PetscCall(DMSetUp(da)); 62 } 63 64 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 65 Create the parent vector 66 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 67 PetscCall(DMCreateGlobalVector(da, &vec_full)); 68 PetscCall(PetscObjectSetName((PetscObject) vec_full, "full_vector")); 69 70 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 71 Populate the 3D vector 72 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 73 PetscCall(DMDAGetCorners(da, &ixs, &iys, &izs, &ixm, &iym, &izm)); 74 if (dim==3) { 75 PetscCall(DMDAVecGetArray(da, vec_full, &vecdata3d)); 76 for (k=izs; k<izs+izm; k++) { 77 for (j=iys; j<iys+iym; j++) { 78 for (i=ixs; i<ixs+ixm; i++) { 79 vecdata3d[k][j][i] = ((i-mx/2)*(j+mx/2))+k*100; 80 } 81 } 82 } 83 PetscCall(DMDAVecRestoreArray(da, vec_full, &vecdata3d)); 84 } else { 85 PetscCall(DMDAVecGetArray(da, vec_full, &vecdata2d)); 86 for (j=iys; j<iys+iym; j++) { 87 for (i=ixs; i<ixs+ixm; i++) { 88 vecdata2d[j][i] = ((i-mx/2)*(j+mx/2)); 89 } 90 } 91 PetscCall(DMDAVecRestoreArray(da, vec_full, &vecdata2d)); 92 } 93 94 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 95 Get an IS corresponding to a 2D slice 96 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 97 if (sliceaxis==0) { 98 lower.i = sliceid; lower.j = 0; lower.k = 0; lower.c = 1; 99 upper.i = sliceid; upper.j = my; upper.k = mz; upper.c = 1; 100 } else if (sliceaxis==1) { 101 lower.i = 0; lower.j = sliceid; lower.k = 0; lower.c = 1; 102 upper.i = mx; upper.j = sliceid; upper.k = mz; upper.c = 1; 103 } else if (sliceaxis==2 && dim==3) { 104 lower.i = 0; lower.j = 0; lower.k = sliceid; lower.c = 1; 105 upper.i = mx; upper.j = my; upper.k = sliceid; upper.c = 1; 106 } 107 PetscCall(DMDACreatePatchIS(da, &lower, &upper, &selectis, patchis_offproc)); 108 PetscCall(ISView(selectis, PETSC_VIEWER_STDOUT_WORLD)); 109 110 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 111 Use the obtained IS to extract the slice as a subvector 112 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 113 PetscCall(VecGetSubVector(vec_full, selectis, &vec_slice)); 114 115 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 116 View the extracted subvector 117 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 118 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_DENSE)); 119 PetscCall(VecView(vec_slice, PETSC_VIEWER_STDOUT_WORLD)); 120 PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD)); 121 122 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 123 Restore subvector, destroy data structures and exit. 124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 125 PetscCall(VecRestoreSubVector(vec_full, selectis, &vec_slice)); 126 127 PetscCall(ISDestroy(&selectis)); 128 PetscCall(DMDestroy(&da)); 129 PetscCall(VecDestroy(&vec_full)); 130 131 PetscCall(PetscFinalize()); 132 return 0; 133 } 134 135 /*TEST 136 137 test: 138 nsize: 1 139 args: -sliceaxis 0 140 141 test: 142 suffix: 2 143 nsize: 2 144 args: -sliceaxis 1 145 146 test: 147 suffix: 3 148 nsize: 4 149 args: -sliceaxis 2 150 151 test: 152 suffix: 4 153 nsize: 2 154 args: -sliceaxis 1 -dim 2 155 156 test: 157 suffix: 5 158 nsize: 3 159 args: -sliceaxis 0 -dim 2 160 161 TEST*/ 162