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