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 PetscSection section; 30 DM cda; 31 DMDABoundaryType bx,by,bz; 32 Vec xcoor; 33 PetscScalar *coors; 34 PetscReal hx,hy,hz_; 35 PetscInt i,j,k,M,N,P,istart,isize,jstart,jsize,kstart,ksize,dim,cnt; 36 PetscErrorCode ierr; 37 38 PetscFunctionBegin; 39 PetscValidHeaderSpecific(da,DM_CLASSID,1); 40 ierr = DMDAGetDimension(da, &dim);CHKERRQ(ierr); 41 if (xmax <= xmin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %G %G",xmin,xmax); 42 if ((ymax <= ymin) && (dim > 1)) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %G %G",ymin,ymax); 43 if ((zmax <= zmin) && (dim > 2)) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"zmax must be larger than zmin %G %G",zmin,zmax); 44 45 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 46 ierr = DMGetDefaultSection(da,§ion);CHKERRQ(ierr); 47 ierr = DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 48 ierr = DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);CHKERRQ(ierr); 49 ierr = DMDAGetCoordinateDA(da, &cda);CHKERRQ(ierr); 50 if (section) { 51 /* This would be better as a vector, but this is compatible */ 52 PetscInt numComp[3] = {1, 1, 1}; 53 PetscInt numVertexDof[3] = {1, 1, 1}; 54 55 ierr = DMDASetFieldName(cda, 0, "x");CHKERRQ(ierr); 56 if (dim > 1) {ierr = DMDASetFieldName(cda, 1, "y");CHKERRQ(ierr);} 57 if (dim > 2) {ierr = DMDASetFieldName(cda, 2, "z");CHKERRQ(ierr);} 58 ierr = DMDACreateSection(cda, numComp, numVertexDof, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 59 } 60 ierr = DMCreateGlobalVector(cda, &xcoor);CHKERRQ(ierr); 61 if (section) { 62 PetscSection csection; 63 PetscInt vStart, vEnd; 64 65 ierr = DMGetDefaultGlobalSection(cda,&csection);CHKERRQ(ierr); 66 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 67 ierr = DMDAGetHeightStratum(da, dim, &vStart, &vEnd);CHKERRQ(ierr); 68 if (bx == DMDA_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M+1); 69 else hx = (xmax-xmin)/(M ? M : 1); 70 if (by == DMDA_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N+1); 71 else hy = (ymax-ymin)/(N ? N : 1); 72 if (bz == DMDA_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P+1); 73 else hz_ = (zmax-zmin)/(P ? P : 1); 74 switch (dim) { 75 case 1: 76 for(i = 0; i < isize+1; ++i) { 77 PetscInt v = i+vStart, dof, off; 78 79 ierr = PetscSectionGetDof(csection, v, &dof);CHKERRQ(ierr); 80 ierr = PetscSectionGetOffset(csection, v, &off);CHKERRQ(ierr); 81 if (off >= 0) { 82 coors[off] = xmin + hx*(i+istart); 83 } 84 } 85 break; 86 case 2: 87 for(j = 0; j < jsize+1; ++j) { 88 for(i = 0; i < isize+1; ++i) { 89 PetscInt v = j*(isize+1)+i+vStart, dof, off; 90 91 ierr = PetscSectionGetDof(csection, v, &dof);CHKERRQ(ierr); 92 ierr = PetscSectionGetOffset(csection, v, &off);CHKERRQ(ierr); 93 if (off >= 0) { 94 coors[off+0] = xmin + hx*(i+istart); 95 coors[off+1] = ymin + hy*(j+jstart); 96 } 97 } 98 } 99 break; 100 case 3: 101 for(k = 0; k < ksize+1; ++k) { 102 for(j = 0; j < jsize+1; ++j) { 103 for(i = 0; i < isize+1; ++i) { 104 PetscInt v = (k*(jsize+1)+j)*(isize+1)+i+vStart, dof, off; 105 106 ierr = PetscSectionGetDof(csection, v, &dof);CHKERRQ(ierr); 107 ierr = PetscSectionGetOffset(csection, v, &off);CHKERRQ(ierr); 108 if (off >= 0) { 109 coors[off+0] = xmin + hx*(i+istart); 110 coors[off+1] = ymin + hy*(j+jstart); 111 coors[off+2] = zmin + hz_*(k+kstart); 112 } 113 } 114 } 115 } 116 break; 117 default: 118 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim); 119 } 120 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 121 ierr = DMDASetCoordinates(da,xcoor);CHKERRQ(ierr); 122 ierr = PetscLogObjectParent(da,xcoor);CHKERRQ(ierr); 123 ierr = VecDestroy(&xcoor);CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 if (dim == 1) { 127 if (bx == DMDA_BOUNDARY_PERIODIC) hx = (xmax-xmin)/M; 128 else hx = (xmax-xmin)/(M-1); 129 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 130 for (i=0; i<isize; i++) { 131 coors[i] = xmin + hx*(i+istart); 132 } 133 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 134 } else if (dim == 2) { 135 if (bx == DMDA_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M); 136 else hx = (xmax-xmin)/(M-1); 137 if (by == DMDA_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N); 138 else hy = (ymax-ymin)/(N-1); 139 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 140 cnt = 0; 141 for (j=0; j<jsize; j++) { 142 for (i=0; i<isize; i++) { 143 coors[cnt++] = xmin + hx*(i+istart); 144 coors[cnt++] = ymin + hy*(j+jstart); 145 } 146 } 147 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 148 } else if (dim == 3) { 149 if (bx == DMDA_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M); 150 else hx = (xmax-xmin)/(M-1); 151 if (by == DMDA_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N); 152 else hy = (ymax-ymin)/(N-1); 153 if (bz == DMDA_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P); 154 else hz_ = (zmax-zmin)/(P-1); 155 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 156 cnt = 0; 157 for (k=0; k<ksize; k++) { 158 for (j=0; j<jsize; j++) { 159 for (i=0; i<isize; i++) { 160 coors[cnt++] = xmin + hx*(i+istart); 161 coors[cnt++] = ymin + hy*(j+jstart); 162 coors[cnt++] = zmin + hz_*(k+kstart); 163 } 164 } 165 } 166 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 167 } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim); 168 ierr = DMDASetCoordinates(da,xcoor);CHKERRQ(ierr); 169 ierr = PetscLogObjectParent(da,xcoor);CHKERRQ(ierr); 170 ierr = VecDestroy(&xcoor);CHKERRQ(ierr); 171 PetscFunctionReturn(0); 172 } 173 174 #undef __FUNCT__ 175 #define __FUNCT__ "VecView_MPI_Draw_DA1d" 176 PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v) 177 { 178 DM da; 179 PetscErrorCode ierr; 180 PetscMPIInt rank,size,tag1,tag2; 181 PetscInt i,n,N,step,istart,isize,j,nbounds; 182 MPI_Status status; 183 PetscReal coors[4],ymin,ymax,min,max,xmin,xmax,tmp,xgtmp; 184 const PetscScalar *array,*xg; 185 PetscDraw draw; 186 PetscBool isnull,showpoints = PETSC_FALSE; 187 MPI_Comm comm; 188 PetscDrawAxis axis; 189 Vec xcoor; 190 DMDABoundaryType bx; 191 const PetscReal *bounds; 192 PetscInt *displayfields; 193 PetscInt k,ndisplayfields; 194 PetscBool flg,hold; 195 196 PetscFunctionBegin; 197 ierr = PetscViewerDrawGetDraw(v,0,&draw);CHKERRQ(ierr); 198 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 199 ierr = PetscViewerDrawGetBounds(v,&nbounds,&bounds);CHKERRQ(ierr); 200 201 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 202 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 203 204 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_vec_mark_points",&showpoints,PETSC_NULL);CHKERRQ(ierr); 205 206 ierr = DMDAGetInfo(da,0,&N,0,0,0,0,0,&step,0,&bx,0,0,0);CHKERRQ(ierr); 207 ierr = DMDAGetCorners(da,&istart,0,0,&isize,0,0);CHKERRQ(ierr); 208 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 209 ierr = VecGetLocalSize(xin,&n);CHKERRQ(ierr); 210 n = n/step; 211 212 /* get coordinates of nodes */ 213 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 214 if (!xcoor) { 215 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);CHKERRQ(ierr); 216 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 217 } 218 ierr = VecGetArrayRead(xcoor,&xg);CHKERRQ(ierr); 219 220 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 221 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 222 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 223 224 /* 225 Determine the min and max x coordinate in plot 226 */ 227 if (!rank) { 228 xmin = PetscRealPart(xg[0]); 229 } 230 if (rank == size-1) { 231 xmax = PetscRealPart(xg[n-1]); 232 } 233 ierr = MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);CHKERRQ(ierr); 234 ierr = MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);CHKERRQ(ierr); 235 236 ierr = PetscMalloc(step*sizeof(PetscInt),&displayfields);CHKERRQ(ierr); 237 for (i=0; i<step; i++) displayfields[i] = i; 238 ndisplayfields = step; 239 ierr = PetscOptionsGetIntArray(PETSC_NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);CHKERRQ(ierr); 240 if (!flg) ndisplayfields = step; 241 for (k=0; k<ndisplayfields; k++) { 242 j = displayfields[k]; 243 ierr = PetscViewerDrawGetDraw(v,k,&draw);CHKERRQ(ierr); 244 ierr = PetscDrawCheckResizedWindow(draw);CHKERRQ(ierr); 245 246 /* 247 Determine the min and max y coordinate in plot 248 */ 249 min = 1.e20; max = -1.e20; 250 for (i=0; i<n; i++) { 251 if (PetscRealPart(array[j+i*step]) < min) min = PetscRealPart(array[j+i*step]); 252 if (PetscRealPart(array[j+i*step]) > max) max = PetscRealPart(array[j+i*step]); 253 } 254 if (min + 1.e-10 > max) { 255 min -= 1.e-5; 256 max += 1.e-5; 257 } 258 if (j < nbounds) { 259 min = PetscMin(min,bounds[2*j]); 260 max = PetscMax(max,bounds[2*j+1]); 261 } 262 263 ierr = MPI_Reduce(&min,&ymin,1,MPIU_REAL,MPIU_MIN,0,comm);CHKERRQ(ierr); 264 ierr = MPI_Reduce(&max,&ymax,1,MPIU_REAL,MPIU_MAX,0,comm);CHKERRQ(ierr); 265 266 ierr = PetscViewerDrawGetHold(v,&hold);CHKERRQ(ierr); 267 if (!hold) { 268 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 269 } 270 ierr = PetscViewerDrawGetDrawAxis(v,k,&axis);CHKERRQ(ierr); 271 ierr = PetscLogObjectParent(draw,axis);CHKERRQ(ierr); 272 if (!rank) { 273 const char *title; 274 275 ierr = PetscDrawAxisSetLimits(axis,xmin,xmax,ymin,ymax);CHKERRQ(ierr); 276 ierr = PetscDrawAxisDraw(axis);CHKERRQ(ierr); 277 ierr = PetscDrawGetCoordinates(draw,coors,coors+1,coors+2,coors+3);CHKERRQ(ierr); 278 ierr = DMDAGetFieldName(da,j,&title);CHKERRQ(ierr); 279 if (title) {ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr);} 280 } 281 ierr = MPI_Bcast(coors,4,MPIU_REAL,0,comm);CHKERRQ(ierr); 282 if (rank) { 283 ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 284 } 285 286 /* draw local part of vector */ 287 ierr = PetscObjectGetNewTag((PetscObject)xin,&tag1);CHKERRQ(ierr); 288 ierr = PetscObjectGetNewTag((PetscObject)xin,&tag2);CHKERRQ(ierr); 289 if (rank < size-1) { /*send value to right */ 290 ierr = MPI_Send((void*)&array[j+(n-1)*step],1,MPIU_REAL,rank+1,tag1,comm);CHKERRQ(ierr); 291 ierr = MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag1,comm);CHKERRQ(ierr); 292 } 293 if (!rank && bx == DMDA_BOUNDARY_PERIODIC && size > 1) { /* first processor sends first value to last */ 294 ierr = MPI_Send((void*)&array[j],1,MPIU_REAL,size-1,tag2,comm);CHKERRQ(ierr); 295 } 296 297 for (i=1; i<n; i++) { 298 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); 299 if (showpoints) { 300 ierr = PetscDrawPoint(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PETSC_DRAW_BLACK);CHKERRQ(ierr); 301 } 302 } 303 if (rank) { /* receive value from left */ 304 ierr = MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag1,comm,&status);CHKERRQ(ierr); 305 ierr = MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag1,comm,&status);CHKERRQ(ierr); 306 ierr = PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);CHKERRQ(ierr); 307 if (showpoints) { 308 ierr = PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);CHKERRQ(ierr); 309 } 310 } 311 if (rank == size-1 && bx == DMDA_BOUNDARY_PERIODIC && size > 1) { 312 ierr = MPI_Recv(&tmp,1,MPIU_REAL,0,tag2,comm,&status);CHKERRQ(ierr); 313 ierr = PetscDrawLine(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),PetscRealPart(xg[n-1]),tmp,PETSC_DRAW_RED);CHKERRQ(ierr); 314 if (showpoints) { 315 ierr = PetscDrawPoint(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),PETSC_DRAW_BLACK);CHKERRQ(ierr); 316 } 317 } 318 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 319 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 320 } 321 ierr = PetscFree(displayfields);CHKERRQ(ierr); 322 ierr = VecRestoreArrayRead(xcoor,&xg);CHKERRQ(ierr); 323 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 324 PetscFunctionReturn(0); 325 } 326 327