1 2 /* 3 Plots vectors obtained with DMDACreate1d() 4 */ 5 6 #include <petscdmda.h> /*I "petscdmda.h" I*/ 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "DMDASetUniformCoordinates" 10 /*@ 11 DMDASetUniformCoordinates - Sets a DMDA coordinates to be a uniform grid 12 13 Collective on DMDA 14 15 Input Parameters: 16 + da - the distributed array object 17 . xmin,xmax - extremes in the x direction 18 . ymin,ymax - extremes in the y direction (use PETSC_NULL for 1 dimensional problems) 19 - zmin,zmax - extremes in the z direction (use PETSC_NULL for 1 or 2 dimensional problems) 20 21 Level: beginner 22 23 .seealso: DMDASetCoordinates(), DMDAGetCoordinates(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 24 25 @*/ 26 PetscErrorCode DMDASetUniformCoordinates(DM da,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax) 27 { 28 MPI_Comm comm; 29 DM cda; 30 DMDABoundaryType bx,by,bz; 31 Vec xcoor; 32 PetscScalar *coors; 33 PetscReal hx,hy,hz_; 34 PetscInt i,j,k,M,N,P,istart,isize,jstart,jsize,kstart,ksize,dim,cnt; 35 PetscErrorCode ierr; 36 37 PetscFunctionBegin; 38 if (xmax <= xmin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %G %G",xmin,xmax); 39 40 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 41 ierr = DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 42 ierr = DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);CHKERRQ(ierr); 43 ierr = DMDAGetCoordinateDA(da, &cda);CHKERRQ(ierr); 44 ierr = DMCreateGlobalVector(cda, &xcoor);CHKERRQ(ierr); 45 if (dim == 1) { 46 if (bx == DMDA_BOUNDARY_PERIODIC) hx = (xmax-xmin)/M; 47 else hx = (xmax-xmin)/(M-1); 48 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 49 for (i=0; i<isize; i++) { 50 coors[i] = xmin + hx*(i+istart); 51 } 52 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 53 } else if (dim == 2) { 54 if (ymax <= ymin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %G %G",ymin,ymax); 55 if (bx == DMDA_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M); 56 else hx = (xmax-xmin)/(M-1); 57 if (by == DMDA_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N); 58 else hy = (ymax-ymin)/(N-1); 59 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 60 cnt = 0; 61 for (j=0; j<jsize; j++) { 62 for (i=0; i<isize; i++) { 63 coors[cnt++] = xmin + hx*(i+istart); 64 coors[cnt++] = ymin + hy*(j+jstart); 65 } 66 } 67 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 68 } else if (dim == 3) { 69 if (ymax <= ymin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %G %G",ymin,ymax); 70 if (zmax <= zmin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"zmax must be larger than zmin %G %G",zmin,zmax); 71 if (bx == DMDA_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M); 72 else hx = (xmax-xmin)/(M-1); 73 if (by == DMDA_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N); 74 else hy = (ymax-ymin)/(N-1); 75 if (bz == DMDA_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P); 76 else hz_ = (zmax-zmin)/(P-1); 77 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 78 cnt = 0; 79 for (k=0; k<ksize; k++) { 80 for (j=0; j<jsize; j++) { 81 for (i=0; i<isize; i++) { 82 coors[cnt++] = xmin + hx*(i+istart); 83 coors[cnt++] = ymin + hy*(j+jstart); 84 coors[cnt++] = zmin + hz_*(k+kstart); 85 } 86 } 87 } 88 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 89 } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim); 90 ierr = DMDASetCoordinates(da,xcoor);CHKERRQ(ierr); 91 ierr = PetscLogObjectParent(da,xcoor);CHKERRQ(ierr); 92 ierr = VecDestroy(&xcoor);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 #undef __FUNCT__ 97 #define __FUNCT__ "VecView_MPI_Draw_DA1d" 98 PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v) 99 { 100 DM da; 101 PetscErrorCode ierr; 102 PetscMPIInt rank,size,tag1,tag2; 103 PetscInt i,n,N,step,istart,isize,j,nbounds; 104 MPI_Status status; 105 PetscReal coors[4],ymin,ymax,min,max,xmin,xmax,tmp,xgtmp; 106 const PetscScalar *array,*xg; 107 PetscDraw draw; 108 PetscBool isnull,showpoints = PETSC_FALSE; 109 MPI_Comm comm; 110 PetscDrawAxis axis; 111 Vec xcoor; 112 DMDABoundaryType bx; 113 const PetscReal *bounds; 114 PetscInt *displayfields; 115 PetscInt k,ndisplayfields; 116 PetscBool flg; 117 118 PetscFunctionBegin; 119 ierr = PetscViewerDrawGetDraw(v,0,&draw);CHKERRQ(ierr); 120 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 121 ierr = PetscViewerDrawGetBounds(v,&nbounds,&bounds);CHKERRQ(ierr); 122 123 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 124 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 125 126 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_vec_mark_points",&showpoints,PETSC_NULL);CHKERRQ(ierr); 127 128 ierr = DMDAGetInfo(da,0,&N,0,0,0,0,0,&step,0,&bx,0,0,0);CHKERRQ(ierr); 129 ierr = DMDAGetCorners(da,&istart,0,0,&isize,0,0);CHKERRQ(ierr); 130 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 131 ierr = VecGetLocalSize(xin,&n);CHKERRQ(ierr); 132 n = n/step; 133 134 /* get coordinates of nodes */ 135 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 136 if (!xcoor) { 137 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);CHKERRQ(ierr); 138 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 139 } 140 ierr = VecGetArrayRead(xcoor,&xg);CHKERRQ(ierr); 141 142 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 143 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 144 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 145 146 /* 147 Determine the min and max x coordinate in plot 148 */ 149 if (!rank) { 150 xmin = PetscRealPart(xg[0]); 151 } 152 if (rank == size-1) { 153 xmax = PetscRealPart(xg[n-1]); 154 } 155 ierr = MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);CHKERRQ(ierr); 156 ierr = MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);CHKERRQ(ierr); 157 158 ierr = PetscMalloc(step*sizeof(PetscInt),&displayfields);CHKERRQ(ierr); 159 for (i=0; i<step; i++) displayfields[i] = i; 160 ndisplayfields = step; 161 ierr = PetscOptionsGetIntArray(PETSC_NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);CHKERRQ(ierr); 162 if (!flg) ndisplayfields = step; 163 for (k=0; k<ndisplayfields; k++) { 164 j = displayfields[k]; 165 ierr = PetscViewerDrawGetDraw(v,k,&draw);CHKERRQ(ierr); 166 ierr = PetscDrawCheckResizedWindow(draw);CHKERRQ(ierr); 167 168 /* 169 Determine the min and max y coordinate in plot 170 */ 171 min = 1.e20; max = -1.e20; 172 for (i=0; i<n; i++) { 173 #if defined(PETSC_USE_COMPLEX) 174 if (PetscRealPart(array[j+i*step]) < min) min = PetscRealPart(array[j+i*step]); 175 if (PetscRealPart(array[j+i*step]) > max) max = PetscRealPart(array[j+i*step]); 176 #else 177 if (array[j+i*step] < min) min = array[j+i*step]; 178 if (array[j+i*step] > max) max = array[j+i*step]; 179 #endif 180 } 181 if (min + 1.e-10 > max) { 182 min -= 1.e-5; 183 max += 1.e-5; 184 } 185 if (j < nbounds) { 186 min = PetscMin(min,bounds[2*j]); 187 max = PetscMax(max,bounds[2*j+1]); 188 } 189 190 ierr = MPI_Reduce(&min,&ymin,1,MPIU_REAL,MPIU_MIN,0,comm);CHKERRQ(ierr); 191 ierr = MPI_Reduce(&max,&ymax,1,MPIU_REAL,MPIU_MAX,0,comm);CHKERRQ(ierr); 192 193 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 194 ierr = PetscViewerDrawGetDrawAxis(v,k,&axis);CHKERRQ(ierr); 195 ierr = PetscLogObjectParent(draw,axis);CHKERRQ(ierr); 196 if (!rank) { 197 const char *title; 198 199 ierr = PetscDrawAxisSetLimits(axis,xmin,xmax,ymin,ymax);CHKERRQ(ierr); 200 ierr = PetscDrawAxisDraw(axis);CHKERRQ(ierr); 201 ierr = PetscDrawGetCoordinates(draw,coors,coors+1,coors+2,coors+3);CHKERRQ(ierr); 202 ierr = DMDAGetFieldName(da,j,&title);CHKERRQ(ierr); 203 if (title) {ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr);} 204 } 205 ierr = MPI_Bcast(coors,4,MPIU_REAL,0,comm);CHKERRQ(ierr); 206 if (rank) { 207 ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 208 } 209 210 /* draw local part of vector */ 211 PetscObjectGetNewTag((PetscObject)xin,&tag1);CHKERRQ(ierr); 212 PetscObjectGetNewTag((PetscObject)xin,&tag2);CHKERRQ(ierr); 213 if (rank < size-1) { /*send value to right */ 214 ierr = MPI_Send((void*)&array[j+(n-1)*step],1,MPIU_REAL,rank+1,tag1,comm);CHKERRQ(ierr); 215 ierr = MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag1,comm);CHKERRQ(ierr); 216 } 217 if (!rank && bx == DMDA_BOUNDARY_PERIODIC && size > 1) { /* first processor sends first value to last */ 218 ierr = MPI_Send((void*)&array[j],1,MPIU_REAL,size-1,tag2,comm);CHKERRQ(ierr); 219 } 220 221 for (i=1; i<n; i++) { 222 #if !defined(PETSC_USE_COMPLEX) 223 ierr = PetscDrawLine(draw,xg[i-1],array[j+step*(i-1)],xg[i],array[j+step*i],PETSC_DRAW_RED);CHKERRQ(ierr); 224 #else 225 ierr = PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+step*i]),PETSC_DRAW_RED);CHKERRQ(ierr); 226 #endif 227 if (showpoints) { 228 ierr = PetscDrawPoint(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PETSC_DRAW_BLACK);CHKERRQ(ierr); 229 } 230 } 231 if (rank) { /* receive value from left */ 232 ierr = MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag1,comm,&status);CHKERRQ(ierr); 233 ierr = MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag1,comm,&status);CHKERRQ(ierr); 234 #if !defined(PETSC_USE_COMPLEX) 235 ierr = PetscDrawLine(draw,xgtmp,tmp,xg[0],array[j],PETSC_DRAW_RED);CHKERRQ(ierr); 236 #else 237 ierr = PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);CHKERRQ(ierr); 238 #endif 239 if (showpoints) { 240 ierr = PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);CHKERRQ(ierr); 241 } 242 } 243 if (rank == size-1 && bx == DMDA_BOUNDARY_PERIODIC && size > 1) { 244 ierr = MPI_Recv(&tmp,1,MPIU_REAL,0,tag2,comm,&status);CHKERRQ(ierr); 245 #if !defined(PETSC_USE_COMPLEX) 246 ierr = PetscDrawLine(draw,xg[n-2],array[j+step*(n-1)],xg[n-1],tmp,PETSC_DRAW_RED);CHKERRQ(ierr); 247 #else 248 ierr = PetscDrawLine(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]), 249 PetscRealPart(xg[n-1]),tmp,PETSC_DRAW_RED);CHKERRQ(ierr); 250 #endif 251 if (showpoints) { 252 ierr = PetscDrawPoint(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),PETSC_DRAW_BLACK);CHKERRQ(ierr); 253 } 254 } 255 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 256 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 257 } 258 ierr = VecRestoreArrayRead(xcoor,&xg);CHKERRQ(ierr); 259 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 260 PetscFunctionReturn(0); 261 } 262 263