1 2 /* 3 Plots vectors obtained with DMDACreate1d() 4 */ 5 6 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 7 8 /*@ 9 DMDASetUniformCoordinates - Sets a DMDA coordinates to be a uniform grid 10 11 Collective on da 12 13 Input Parameters: 14 + da - the distributed array object 15 . xmin,xmax - extremes in the x direction 16 . ymin,ymax - extremes in the y direction (value ignored for 1 dimensional problems) 17 - zmin,zmax - extremes in the z direction (value ignored for 1 or 2 dimensional problems) 18 19 Level: beginner 20 21 .seealso: DMSetCoordinates(), DMGetCoordinates(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMStagSetUniformCoordinates() 22 23 @*/ 24 PetscErrorCode DMDASetUniformCoordinates(DM da,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax) 25 { 26 MPI_Comm comm; 27 DM cda; 28 DM_DA *dd = (DM_DA*)da->data; 29 DMBoundaryType bx,by,bz; 30 Vec xcoor; 31 PetscScalar *coors; 32 PetscReal hx,hy,hz_; 33 PetscInt i,j,k,M,N,P,istart,isize,jstart,jsize,kstart,ksize,dim,cnt; 34 PetscErrorCode ierr; 35 36 PetscFunctionBegin; 37 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 38 if (!dd->gtol) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set coordinates until after DMDA has been setup"); 39 ierr = DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 40 if (xmax < xmin) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %g %g",(double)xmin,(double)xmax); 41 if ((dim > 1) && (ymax < ymin)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %g %g",(double)ymin,(double)ymax); 42 if ((dim > 2) && (zmax < zmin)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"zmax must be larger than zmin %g %g",(double)zmin,(double)zmax); 43 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 44 ierr = DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);CHKERRQ(ierr); 45 ierr = DMGetCoordinateDM(da, &cda);CHKERRQ(ierr); 46 ierr = DMCreateGlobalVector(cda, &xcoor);CHKERRQ(ierr); 47 if (dim == 1) { 48 if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/M; 49 else hx = (xmax-xmin)/(M-1); 50 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 51 for (i=0; i<isize; i++) { 52 coors[i] = xmin + hx*(i+istart); 53 } 54 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 55 } else if (dim == 2) { 56 if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M); 57 else hx = (xmax-xmin)/(M-1); 58 if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N); 59 else hy = (ymax-ymin)/(N-1); 60 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 61 cnt = 0; 62 for (j=0; j<jsize; j++) { 63 for (i=0; i<isize; i++) { 64 coors[cnt++] = xmin + hx*(i+istart); 65 coors[cnt++] = ymin + hy*(j+jstart); 66 } 67 } 68 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 69 } else if (dim == 3) { 70 if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M); 71 else hx = (xmax-xmin)/(M-1); 72 if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N); 73 else hy = (ymax-ymin)/(N-1); 74 if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P); 75 else hz_ = (zmax-zmin)/(P-1); 76 ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr); 77 cnt = 0; 78 for (k=0; k<ksize; k++) { 79 for (j=0; j<jsize; j++) { 80 for (i=0; i<isize; i++) { 81 coors[cnt++] = xmin + hx*(i+istart); 82 coors[cnt++] = ymin + hy*(j+jstart); 83 coors[cnt++] = zmin + hz_*(k+kstart); 84 } 85 } 86 } 87 ierr = VecRestoreArray(xcoor,&coors);CHKERRQ(ierr); 88 } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim); 89 ierr = DMSetCoordinates(da,xcoor);CHKERRQ(ierr); 90 ierr = PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);CHKERRQ(ierr); 91 ierr = VecDestroy(&xcoor);CHKERRQ(ierr); 92 PetscFunctionReturn(0); 93 } 94 95 /* 96 Allows a user to select a subset of the fields to be drawn by VecView() when the vector comes from a DMDA 97 */ 98 PetscErrorCode DMDASelectFields(DM da,PetscInt *outfields,PetscInt **fields) 99 { 100 PetscErrorCode ierr; 101 PetscInt step,ndisplayfields,*displayfields,k,j; 102 PetscBool flg; 103 104 PetscFunctionBegin; 105 ierr = DMDAGetInfo(da,0,0,0,0,0,0,0,&step,0,0,0,0,0);CHKERRQ(ierr); 106 ierr = PetscMalloc1(step,&displayfields);CHKERRQ(ierr); 107 for (k=0; k<step; k++) displayfields[k] = k; 108 ndisplayfields = step; 109 ierr = PetscOptionsGetIntArray(NULL,NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);CHKERRQ(ierr); 110 if (!ndisplayfields) ndisplayfields = step; 111 if (!flg) { 112 char **fields; 113 const char *fieldname; 114 PetscInt nfields = step; 115 ierr = PetscMalloc1(step,&fields);CHKERRQ(ierr); 116 ierr = PetscOptionsGetStringArray(NULL,NULL,"-draw_fields_by_name",fields,&nfields,&flg);CHKERRQ(ierr); 117 if (flg) { 118 ndisplayfields = 0; 119 for (k=0; k<nfields;k++) { 120 for (j=0; j<step; j++) { 121 ierr = DMDAGetFieldName(da,j,&fieldname);CHKERRQ(ierr); 122 ierr = PetscStrcmp(fieldname,fields[k],&flg);CHKERRQ(ierr); 123 if (flg) { 124 goto found; 125 } 126 } 127 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Unknown fieldname %s",fields[k]); 128 found: displayfields[ndisplayfields++] = j; 129 } 130 } 131 for (k=0; k<nfields; k++) { 132 ierr = PetscFree(fields[k]);CHKERRQ(ierr); 133 } 134 ierr = PetscFree(fields);CHKERRQ(ierr); 135 } 136 *fields = displayfields; 137 *outfields = ndisplayfields; 138 PetscFunctionReturn(0); 139 } 140 141 #include <petscdraw.h> 142 143 PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v) 144 { 145 DM da; 146 PetscErrorCode ierr; 147 PetscMPIInt rank,size,tag; 148 PetscInt i,n,N,dof,istart,isize,j,nbounds; 149 MPI_Status status; 150 PetscReal min,max,xmin = 0.0,xmax = 0.0,tmp = 0.0,xgtmp = 0.0; 151 const PetscScalar *array,*xg; 152 PetscDraw draw; 153 PetscBool isnull,useports = PETSC_FALSE,showmarkers = PETSC_FALSE; 154 MPI_Comm comm; 155 PetscDrawAxis axis; 156 Vec xcoor; 157 DMBoundaryType bx; 158 const char *tlabel = NULL,*xlabel = NULL; 159 const PetscReal *bounds; 160 PetscInt *displayfields; 161 PetscInt k,ndisplayfields; 162 PetscBool hold; 163 PetscDrawViewPorts *ports = NULL; 164 PetscViewerFormat format; 165 166 PetscFunctionBegin; 167 ierr = PetscViewerDrawGetDraw(v,0,&draw);CHKERRQ(ierr); 168 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 169 if (isnull) PetscFunctionReturn(0); 170 ierr = PetscViewerDrawGetBounds(v,&nbounds,&bounds);CHKERRQ(ierr); 171 172 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 173 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 174 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 175 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 176 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 177 178 ierr = PetscOptionsGetBool(NULL,NULL,"-draw_vec_use_markers",&showmarkers,NULL);CHKERRQ(ierr); 179 180 ierr = DMDAGetInfo(da,NULL,&N,NULL,NULL,NULL,NULL,NULL,&dof,NULL,&bx,NULL,NULL,NULL);CHKERRQ(ierr); 181 ierr = DMDAGetCorners(da,&istart,NULL,NULL,&isize,NULL,NULL);CHKERRQ(ierr); 182 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 183 ierr = VecGetLocalSize(xin,&n);CHKERRQ(ierr); 184 n = n/dof; 185 186 /* Get coordinates of nodes */ 187 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 188 if (!xcoor) { 189 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);CHKERRQ(ierr); 190 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 191 } 192 ierr = VecGetArrayRead(xcoor,&xg);CHKERRQ(ierr); 193 ierr = DMDAGetCoordinateName(da,0,&xlabel);CHKERRQ(ierr); 194 195 /* Determine the min and max coordinate in plot */ 196 if (!rank) xmin = PetscRealPart(xg[0]); 197 if (rank == size-1) xmax = PetscRealPart(xg[n-1]); 198 ierr = MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);CHKERRQ(ierr); 199 ierr = MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);CHKERRQ(ierr); 200 201 ierr = DMDASelectFields(da,&ndisplayfields,&displayfields);CHKERRQ(ierr); 202 ierr = PetscViewerGetFormat(v,&format);CHKERRQ(ierr); 203 ierr = PetscOptionsGetBool(NULL,NULL,"-draw_ports",&useports,NULL);CHKERRQ(ierr); 204 if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE; 205 if (useports) { 206 ierr = PetscViewerDrawGetDraw(v,0,&draw);CHKERRQ(ierr); 207 ierr = PetscViewerDrawGetDrawAxis(v,0,&axis);CHKERRQ(ierr); 208 ierr = PetscDrawCheckResizedWindow(draw);CHKERRQ(ierr); 209 ierr = PetscDrawClear(draw);CHKERRQ(ierr); 210 ierr = PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);CHKERRQ(ierr); 211 } 212 213 /* Loop over each field; drawing each in a different window */ 214 for (k=0; k<ndisplayfields; k++) { 215 j = displayfields[k]; 216 217 /* determine the min and max value in plot */ 218 ierr = VecStrideMin(xin,j,NULL,&min);CHKERRQ(ierr); 219 ierr = VecStrideMax(xin,j,NULL,&max);CHKERRQ(ierr); 220 if (j < nbounds) { 221 min = PetscMin(min,bounds[2*j]); 222 max = PetscMax(max,bounds[2*j+1]); 223 } 224 if (min == max) { 225 min -= 1.e-5; 226 max += 1.e-5; 227 } 228 229 if (useports) { 230 ierr = PetscDrawViewPortsSet(ports,k);CHKERRQ(ierr); 231 ierr = DMDAGetFieldName(da,j,&tlabel);CHKERRQ(ierr); 232 } else { 233 const char *title; 234 ierr = PetscViewerDrawGetHold(v,&hold);CHKERRQ(ierr); 235 ierr = PetscViewerDrawGetDraw(v,k,&draw);CHKERRQ(ierr); 236 ierr = PetscViewerDrawGetDrawAxis(v,k,&axis);CHKERRQ(ierr); 237 ierr = DMDAGetFieldName(da,j,&title);CHKERRQ(ierr); 238 if (title) {ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr);} 239 ierr = PetscDrawCheckResizedWindow(draw);CHKERRQ(ierr); 240 if (!hold) {ierr = PetscDrawClear(draw);CHKERRQ(ierr);} 241 } 242 ierr = PetscDrawAxisSetLabels(axis,tlabel,xlabel,NULL);CHKERRQ(ierr); 243 ierr = PetscDrawAxisSetLimits(axis,xmin,xmax,min,max);CHKERRQ(ierr); 244 ierr = PetscDrawAxisDraw(axis);CHKERRQ(ierr); 245 246 /* draw local part of vector */ 247 ierr = PetscObjectGetNewTag((PetscObject)xin,&tag);CHKERRQ(ierr); 248 if (rank < size-1) { /*send value to right */ 249 ierr = MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag,comm);CHKERRQ(ierr); 250 ierr = MPI_Send((void*)&array[j+(n-1)*dof],1,MPIU_REAL,rank+1,tag,comm);CHKERRQ(ierr); 251 } 252 if (rank) { /* receive value from left */ 253 ierr = MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag,comm,&status);CHKERRQ(ierr); 254 ierr = MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag,comm,&status);CHKERRQ(ierr); 255 } 256 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 257 if (rank) { 258 ierr = PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);CHKERRQ(ierr); 259 if (showmarkers) {ierr = PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);CHKERRQ(ierr);} 260 } 261 for (i=1; i<n; i++) { 262 ierr = PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+dof*i]),PETSC_DRAW_RED);CHKERRQ(ierr); 263 if (showmarkers) {ierr = PetscDrawMarker(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PETSC_DRAW_BLACK);CHKERRQ(ierr);} 264 } 265 if (rank == size-1) { 266 if (showmarkers) {ierr = PetscDrawMarker(draw,PetscRealPart(xg[n-1]),PetscRealPart(array[j+dof*(n-1)]),PETSC_DRAW_BLACK);CHKERRQ(ierr);} 267 } 268 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 269 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 270 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 271 if (!useports) {ierr = PetscDrawSave(draw);CHKERRQ(ierr);} 272 } 273 if (useports) { 274 ierr = PetscViewerDrawGetDraw(v,0,&draw);CHKERRQ(ierr); 275 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 276 } 277 278 ierr = PetscDrawViewPortsDestroy(ports);CHKERRQ(ierr); 279 ierr = PetscFree(displayfields);CHKERRQ(ierr); 280 ierr = VecRestoreArrayRead(xcoor,&xg);CHKERRQ(ierr); 281 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 282 PetscFunctionReturn(0); 283 } 284 285