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