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