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