1 2 /* 3 Plots vectors obtained with DMDACreate1d() 4 */ 5 6 #include <petsc/private/dmdaimpl.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 (value ignored for 1 dimensional problems) 19 - zmin,zmax - extremes in the z direction (value ignored for 1 or 2 dimensional problems) 20 21 Level: beginner 22 23 .seealso: DMSetCoordinates(), DMGetCoordinates(), 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 DMBoundaryType 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 = DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 41 if (xmax < xmin) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %g %g",(double)xmin,(double)xmax); 42 if ((ymax < ymin) && (dim > 1)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %g %g",(double)ymin,(double)ymax); 43 if ((zmax < zmin) && (dim > 2)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"zmax must be larger than zmin %g %g",(double)zmin,(double)zmax); 44 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 45 ierr = DMGetDefaultSection(da,§ion);CHKERRQ(ierr); 46 ierr = DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);CHKERRQ(ierr); 47 ierr = DMGetCoordinateDM(da, &cda);CHKERRQ(ierr); 48 if (section) { 49 /* This would be better as a vector, but this is compatible */ 50 PetscInt numComp[3] = {1, 1, 1}; 51 PetscInt numVertexDof[3] = {1, 1, 1}; 52 53 ierr = DMDASetFieldName(cda, 0, "x");CHKERRQ(ierr); 54 if (dim > 1) {ierr = DMDASetFieldName(cda, 1, "y");CHKERRQ(ierr);} 55 if (dim > 2) {ierr = DMDASetFieldName(cda, 2, "z");CHKERRQ(ierr);} 56 ierr = DMDACreateSection(cda, numComp, numVertexDof, NULL, NULL);CHKERRQ(ierr); 57 } 58 ierr = DMCreateGlobalVector(cda, &xcoor);CHKERRQ(ierr); 59 if (section) { 60 PetscSection csection; 61 PetscInt vStart, vEnd; 62 63 ierr = DMGetDefaultGlobalSection(cda,&csection);CHKERRQ(ierr); 64 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 65 ierr = DMDAGetHeightStratum(da, dim, &vStart, &vEnd);CHKERRQ(ierr); 66 if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M+1); 67 else hx = (xmax-xmin)/(M ? M : 1); 68 if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N+1); 69 else hy = (ymax-ymin)/(N ? N : 1); 70 if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P+1); 71 else hz_ = (zmax-zmin)/(P ? P : 1); 72 switch (dim) { 73 case 1: 74 for (i = 0; i < isize+1; ++i) { 75 PetscInt v = i+vStart, dof, off; 76 77 ierr = PetscSectionGetDof(csection, v, &dof);CHKERRQ(ierr); 78 ierr = PetscSectionGetOffset(csection, v, &off);CHKERRQ(ierr); 79 if (off >= 0) { 80 coors[off] = xmin + hx*(i+istart); 81 } 82 } 83 break; 84 case 2: 85 for (j = 0; j < jsize+1; ++j) { 86 for (i = 0; i < isize+1; ++i) { 87 PetscInt v = j*(isize+1)+i+vStart, dof, off; 88 89 ierr = PetscSectionGetDof(csection, v, &dof);CHKERRQ(ierr); 90 ierr = PetscSectionGetOffset(csection, v, &off);CHKERRQ(ierr); 91 if (off >= 0) { 92 coors[off+0] = xmin + hx*(i+istart); 93 coors[off+1] = ymin + hy*(j+jstart); 94 } 95 } 96 } 97 break; 98 case 3: 99 for (k = 0; k < ksize+1; ++k) { 100 for (j = 0; j < jsize+1; ++j) { 101 for (i = 0; i < isize+1; ++i) { 102 PetscInt v = (k*(jsize+1)+j)*(isize+1)+i+vStart, dof, off; 103 104 ierr = PetscSectionGetDof(csection, v, &dof);CHKERRQ(ierr); 105 ierr = PetscSectionGetOffset(csection, v, &off);CHKERRQ(ierr); 106 if (off >= 0) { 107 coors[off+0] = xmin + hx*(i+istart); 108 coors[off+1] = ymin + hy*(j+jstart); 109 coors[off+2] = zmin + hz_*(k+kstart); 110 } 111 } 112 } 113 } 114 break; 115 default: 116 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim); 117 } 118 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 119 ierr = DMSetCoordinates(da,xcoor);CHKERRQ(ierr); 120 ierr = PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);CHKERRQ(ierr); 121 ierr = VecDestroy(&xcoor);CHKERRQ(ierr); 122 PetscFunctionReturn(0); 123 } 124 if (dim == 1) { 125 if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/M; 126 else hx = (xmax-xmin)/(M-1); 127 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 128 for (i=0; i<isize; i++) { 129 coors[i] = xmin + hx*(i+istart); 130 } 131 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 132 } else if (dim == 2) { 133 if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M); 134 else hx = (xmax-xmin)/(M-1); 135 if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N); 136 else hy = (ymax-ymin)/(N-1); 137 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 138 cnt = 0; 139 for (j=0; j<jsize; j++) { 140 for (i=0; i<isize; i++) { 141 coors[cnt++] = xmin + hx*(i+istart); 142 coors[cnt++] = ymin + hy*(j+jstart); 143 } 144 } 145 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 146 } else if (dim == 3) { 147 if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M); 148 else hx = (xmax-xmin)/(M-1); 149 if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N); 150 else hy = (ymax-ymin)/(N-1); 151 if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P); 152 else hz_ = (zmax-zmin)/(P-1); 153 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 154 cnt = 0; 155 for (k=0; k<ksize; k++) { 156 for (j=0; j<jsize; j++) { 157 for (i=0; i<isize; i++) { 158 coors[cnt++] = xmin + hx*(i+istart); 159 coors[cnt++] = ymin + hy*(j+jstart); 160 coors[cnt++] = zmin + hz_*(k+kstart); 161 } 162 } 163 } 164 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 165 } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim); 166 ierr = DMSetCoordinates(da,xcoor);CHKERRQ(ierr); 167 ierr = PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);CHKERRQ(ierr); 168 ierr = VecDestroy(&xcoor);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "DMDASelectFields" 174 /* 175 Allows a user to select a subset of the fields to be drawn by VecView() when the vector comes from a DMDA 176 */ 177 PetscErrorCode DMDASelectFields(DM da,PetscInt *outfields,PetscInt **fields) 178 { 179 PetscErrorCode ierr; 180 PetscInt step,ndisplayfields,*displayfields,k,j; 181 PetscBool flg; 182 183 PetscFunctionBegin; 184 ierr = DMDAGetInfo(da,0,0,0,0,0,0,0,&step,0,0,0,0,0);CHKERRQ(ierr); 185 ierr = PetscMalloc1(step,&displayfields);CHKERRQ(ierr); 186 for (k=0; k<step; k++) displayfields[k] = k; 187 ndisplayfields = step; 188 ierr = PetscOptionsGetIntArray(NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);CHKERRQ(ierr); 189 if (!ndisplayfields) ndisplayfields = step; 190 if (!flg) { 191 char **fields; 192 const char *fieldname; 193 PetscInt nfields = step; 194 ierr = PetscMalloc1(step,&fields);CHKERRQ(ierr); 195 ierr = PetscOptionsGetStringArray(NULL,"-draw_fields_by_name",fields,&nfields,&flg);CHKERRQ(ierr); 196 if (flg) { 197 ndisplayfields = 0; 198 for (k=0; k<nfields;k++) { 199 for (j=0; j<step; j++) { 200 ierr = DMDAGetFieldName(da,j,&fieldname);CHKERRQ(ierr); 201 ierr = PetscStrcmp(fieldname,fields[k],&flg);CHKERRQ(ierr); 202 if (flg) { 203 goto found; 204 } 205 } 206 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Unknown fieldname %s",fields[k]); 207 found: displayfields[ndisplayfields++] = j; 208 } 209 } 210 for (k=0; k<nfields; k++) { 211 ierr = PetscFree(fields[k]);CHKERRQ(ierr); 212 } 213 ierr = PetscFree(fields);CHKERRQ(ierr); 214 } 215 *fields = displayfields; 216 *outfields = ndisplayfields; 217 PetscFunctionReturn(0); 218 } 219 220 #include <petscdraw.h> 221 #if defined(PETSC_HAVE_SETJMP_H) && defined(PETSC_HAVE_X) 222 #include <X11/Xlib.h> 223 #include <X11/Xutil.h> 224 #include <setjmp.h> 225 static jmp_buf PetscXIOErrorJumpBuf; 226 static void PetscXIOHandler(Display *dpy) 227 { 228 longjmp(PetscXIOErrorJumpBuf, 1); 229 } 230 #endif 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "VecView_MPI_Draw_DA1d" 234 PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v) 235 { 236 DM da; 237 PetscErrorCode ierr; 238 PetscMPIInt rank,size,tag1,tag2; 239 PetscInt i,n,N,step,istart,isize,j,nbounds; 240 MPI_Status status; 241 PetscReal coors[4],ymin,ymax,min,max,xmin = 0.0,xmax = 0.0,tmp = 0.0,xgtmp = 0.0; 242 const PetscScalar *array,*xg; 243 PetscDraw draw; 244 PetscBool isnull,showmarkers = PETSC_FALSE; 245 MPI_Comm comm; 246 PetscDrawAxis axis; 247 Vec xcoor; 248 DMBoundaryType bx; 249 const PetscReal *bounds; 250 PetscInt *displayfields; 251 PetscInt k,ndisplayfields; 252 PetscBool hold; 253 254 PetscFunctionBegin; 255 ierr = PetscViewerDrawGetDraw(v,0,&draw);CHKERRQ(ierr); 256 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 257 ierr = PetscViewerDrawGetBounds(v,&nbounds,&bounds);CHKERRQ(ierr); 258 259 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 260 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 261 262 ierr = PetscOptionsGetBool(NULL,"-draw_vec_use_markers",&showmarkers,NULL);CHKERRQ(ierr); 263 264 ierr = DMDAGetInfo(da,0,&N,0,0,0,0,0,&step,0,&bx,0,0,0);CHKERRQ(ierr); 265 ierr = DMDAGetCorners(da,&istart,0,0,&isize,0,0);CHKERRQ(ierr); 266 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 267 ierr = VecGetLocalSize(xin,&n);CHKERRQ(ierr); 268 n = n/step; 269 270 /* get coordinates of nodes */ 271 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 272 if (!xcoor) { 273 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);CHKERRQ(ierr); 274 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 275 } 276 ierr = VecGetArrayRead(xcoor,&xg);CHKERRQ(ierr); 277 278 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 279 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 280 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 281 282 /* 283 Determine the min and max x coordinate in plot 284 */ 285 if (!rank) { 286 xmin = PetscRealPart(xg[0]); 287 } 288 if (rank == size-1) { 289 xmax = PetscRealPart(xg[n-1]); 290 } 291 ierr = MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);CHKERRQ(ierr); 292 ierr = MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);CHKERRQ(ierr); 293 294 ierr = DMDASelectFields(da,&ndisplayfields,&displayfields);CHKERRQ(ierr); 295 #if defined(PETSC_HAVE_SETJMP_H) && defined(PETSC_HAVE_X) 296 if (!setjmp(PetscXIOErrorJumpBuf)) XSetIOErrorHandler((XIOErrorHandler)PetscXIOHandler); 297 else { 298 XSetIOErrorHandler(NULL); 299 ierr = PetscDrawSetType(draw,PETSC_DRAW_NULL);CHKERRQ(ierr); 300 PetscFunctionReturn(0); 301 } 302 #endif 303 for (k=0; k<ndisplayfields; k++) { 304 j = displayfields[k]; 305 ierr = PetscViewerDrawGetDraw(v,k,&draw);CHKERRQ(ierr); 306 ierr = PetscDrawCheckResizedWindow(draw);CHKERRQ(ierr); 307 308 /* 309 Determine the min and max y coordinate in plot 310 */ 311 min = 1.e20; max = -1.e20; 312 for (i=0; i<n; i++) { 313 if (PetscRealPart(array[j+i*step]) < min) min = PetscRealPart(array[j+i*step]); 314 if (PetscRealPart(array[j+i*step]) > max) max = PetscRealPart(array[j+i*step]); 315 } 316 if (min + 1.e-10 > max) { 317 min -= 1.e-5; 318 max += 1.e-5; 319 } 320 if (j < nbounds) { 321 min = PetscMin(min,bounds[2*j]); 322 max = PetscMax(max,bounds[2*j+1]); 323 } 324 325 ierr = MPI_Reduce(&min,&ymin,1,MPIU_REAL,MPIU_MIN,0,comm);CHKERRQ(ierr); 326 ierr = MPI_Reduce(&max,&ymax,1,MPIU_REAL,MPIU_MAX,0,comm);CHKERRQ(ierr); 327 328 ierr = PetscViewerDrawGetHold(v,&hold);CHKERRQ(ierr); 329 if (!hold) { 330 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 331 } 332 ierr = PetscViewerDrawGetDrawAxis(v,k,&axis);CHKERRQ(ierr); 333 ierr = PetscLogObjectParent((PetscObject)draw,(PetscObject)axis);CHKERRQ(ierr); 334 if (!rank) { 335 const char *title; 336 337 ierr = PetscDrawAxisSetLimits(axis,xmin,xmax,ymin,ymax);CHKERRQ(ierr); 338 ierr = PetscDrawAxisDraw(axis);CHKERRQ(ierr); 339 ierr = PetscDrawGetCoordinates(draw,coors,coors+1,coors+2,coors+3);CHKERRQ(ierr); 340 ierr = DMDAGetFieldName(da,j,&title);CHKERRQ(ierr); 341 if (title) {ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr);} 342 } 343 ierr = MPI_Bcast(coors,4,MPIU_REAL,0,comm);CHKERRQ(ierr); 344 if (rank) { 345 ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 346 } 347 348 /* draw local part of vector */ 349 ierr = PetscObjectGetNewTag((PetscObject)xin,&tag1);CHKERRQ(ierr); 350 ierr = PetscObjectGetNewTag((PetscObject)xin,&tag2);CHKERRQ(ierr); 351 if (rank < size-1) { /*send value to right */ 352 ierr = MPI_Send((void*)&array[j+(n-1)*step],1,MPIU_REAL,rank+1,tag1,comm);CHKERRQ(ierr); 353 ierr = MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag1,comm);CHKERRQ(ierr); 354 } 355 if (!rank && bx == DM_BOUNDARY_PERIODIC && size > 1) { /* first processor sends first value to last */ 356 ierr = MPI_Send((void*)&array[j],1,MPIU_REAL,size-1,tag2,comm);CHKERRQ(ierr); 357 } 358 359 for (i=1; i<n; i++) { 360 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); 361 if (showmarkers) { 362 ierr = PetscDrawMarker(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PETSC_DRAW_BLACK);CHKERRQ(ierr); 363 } 364 } 365 if (rank) { /* receive value from left */ 366 ierr = MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag1,comm,&status);CHKERRQ(ierr); 367 ierr = MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag1,comm,&status);CHKERRQ(ierr); 368 ierr = PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);CHKERRQ(ierr); 369 if (showmarkers) { 370 ierr = PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);CHKERRQ(ierr); 371 } 372 } 373 if (rank == size-1 && bx == DM_BOUNDARY_PERIODIC && size > 1) { 374 ierr = MPI_Recv(&tmp,1,MPIU_REAL,0,tag2,comm,&status);CHKERRQ(ierr); 375 /* If the mesh is not uniform we do not know the mesh spacing between the last point on the right and the first ghost point */ 376 ierr = PetscDrawLine(draw,PetscRealPart(xg[n-1]),PetscRealPart(array[j+step*(n-1)]),PetscRealPart(xg[n-1]+(xg[n-1]-xg[n-2])),tmp,PETSC_DRAW_RED);CHKERRQ(ierr); 377 if (showmarkers) { 378 ierr = PetscDrawMarker(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),PETSC_DRAW_BLACK);CHKERRQ(ierr); 379 } 380 } 381 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 382 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 383 } 384 #if defined(PETSC_HAVE_SETJMP_H) && defined(PETSC_HAVE_X) 385 XSetIOErrorHandler(NULL); 386 #endif 387 ierr = PetscFree(displayfields);CHKERRQ(ierr); 388 ierr = VecRestoreArrayRead(xcoor,&xg);CHKERRQ(ierr); 389 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 390 PetscFunctionReturn(0); 391 } 392 393