Actual source code: gr1.c


  2: /*
  3:    Plots vectors obtained with DMDACreate1d()
  4: */

  6: #include <petsc/private/dmdaimpl.h>

  8: /*@
  9:     DMDASetUniformCoordinates - Sets a DMDA coordinates to be a uniform grid

 11:   Collective on da

 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)

 19:   Level: beginner

 21: .seealso: DMSetCoordinates(), DMGetCoordinates(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMStagSetUniformCoordinates()

 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;

 37:   DMDAGetInfo(da,&dim,&M,&N,&P,NULL,NULL,NULL,NULL,NULL,&bx,&by,&bz,NULL);
 41:   PetscObjectGetComm((PetscObject)da,&comm);
 42:   DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);
 43:   DMGetCoordinateDM(da, &cda);
 44:   DMCreateGlobalVector(cda, &xcoor);
 45:   if (dim == 1) {
 46:     if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/M;
 47:     else hx = (xmax-xmin)/(M-1);
 48:     VecGetArray(xcoor,&coors);
 49:     for (i=0; i<isize; i++) {
 50:       coors[i] = xmin + hx*(i+istart);
 51:     }
 52:     VecRestoreArray(xcoor,&coors);
 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:     VecGetArray(xcoor,&coors);
 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:     VecRestoreArray(xcoor,&coors);
 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:     VecGetArray(xcoor,&coors);
 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:     VecRestoreArray(xcoor,&coors);
 86:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D",dim);
 87:   DMSetCoordinates(da,xcoor);
 88:   PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);
 89:   VecDestroy(&xcoor);
 90:   return 0;
 91: }

 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:   PetscInt       step,ndisplayfields,*displayfields,k,j;
 99:   PetscBool      flg;

101:   DMDAGetInfo(da,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&step,NULL,NULL,NULL,NULL,NULL);
102:   PetscMalloc1(step,&displayfields);
103:   for (k=0; k<step; k++) displayfields[k] = k;
104:   ndisplayfields = step;
105:   PetscOptionsGetIntArray(NULL,NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);
106:   if (!ndisplayfields) ndisplayfields = step;
107:   if (!flg) {
108:     char       **fields;
109:     const char *fieldname;
110:     PetscInt   nfields = step;
111:     PetscMalloc1(step,&fields);
112:     PetscOptionsGetStringArray(NULL,NULL,"-draw_fields_by_name",fields,&nfields,&flg);
113:     if (flg) {
114:       ndisplayfields = 0;
115:       for (k=0; k<nfields;k++) {
116:         for (j=0; j<step; j++) {
117:           DMDAGetFieldName(da,j,&fieldname);
118:           PetscStrcmp(fieldname,fields[k],&flg);
119:           if (flg) {
120:             goto found;
121:           }
122:         }
123:         SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Unknown fieldname %s",fields[k]);
124: found:  displayfields[ndisplayfields++] = j;
125:       }
126:     }
127:     for (k=0; k<nfields; k++) {
128:       PetscFree(fields[k]);
129:     }
130:     PetscFree(fields);
131:   }
132:   *fields    = displayfields;
133:   *outfields = ndisplayfields;
134:   return 0;
135: }

137: #include <petscdraw.h>

139: PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v)
140: {
141:   DM                da;
142:   PetscMPIInt       rank,size,tag;
143:   PetscInt          i,n,N,dof,istart,isize,j,nbounds;
144:   MPI_Status        status;
145:   PetscReal         min,max,xmin = 0.0,xmax = 0.0,tmp = 0.0,xgtmp = 0.0;
146:   const PetscScalar *array,*xg;
147:   PetscDraw         draw;
148:   PetscBool         isnull,useports = PETSC_FALSE,showmarkers = PETSC_FALSE;
149:   MPI_Comm          comm;
150:   PetscDrawAxis     axis;
151:   Vec               xcoor;
152:   DMBoundaryType    bx;
153:   const char        *tlabel = NULL,*xlabel = NULL;
154:   const PetscReal   *bounds;
155:   PetscInt          *displayfields;
156:   PetscInt          k,ndisplayfields;
157:   PetscBool         hold;
158:   PetscDrawViewPorts *ports = NULL;
159:   PetscViewerFormat  format;

161:   PetscViewerDrawGetDraw(v,0,&draw);
162:   PetscDrawIsNull(draw,&isnull);
163:   if (isnull) return 0;
164:   PetscViewerDrawGetBounds(v,&nbounds,&bounds);

166:   VecGetDM(xin,&da);
168:   PetscObjectGetComm((PetscObject)xin,&comm);
169:   MPI_Comm_size(comm,&size);
170:   MPI_Comm_rank(comm,&rank);

172:   PetscOptionsGetBool(NULL,NULL,"-draw_vec_use_markers",&showmarkers,NULL);

174:   DMDAGetInfo(da,NULL,&N,NULL,NULL,NULL,NULL,NULL,&dof,NULL,&bx,NULL,NULL,NULL);
175:   DMDAGetCorners(da,&istart,NULL,NULL,&isize,NULL,NULL);
176:   VecGetArrayRead(xin,&array);
177:   VecGetLocalSize(xin,&n);
178:   n    = n/dof;

180:   /* Get coordinates of nodes */
181:   DMGetCoordinates(da,&xcoor);
182:   if (!xcoor) {
183:     DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
184:     DMGetCoordinates(da,&xcoor);
185:   }
186:   VecGetArrayRead(xcoor,&xg);
187:   DMDAGetCoordinateName(da,0,&xlabel);

189:   /* Determine the min and max coordinate in plot */
190:   if (rank == 0) xmin = PetscRealPart(xg[0]);
191:   if (rank == size-1) xmax = PetscRealPart(xg[n-1]);
192:   MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
193:   MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);

195:   DMDASelectFields(da,&ndisplayfields,&displayfields);
196:   PetscViewerGetFormat(v,&format);
197:   PetscOptionsGetBool(NULL,NULL,"-draw_ports",&useports,NULL);
198:   if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
199:   if (useports) {
200:     PetscViewerDrawGetDraw(v,0,&draw);
201:     PetscViewerDrawGetDrawAxis(v,0,&axis);
202:     PetscDrawCheckResizedWindow(draw);
203:     PetscDrawClear(draw);
204:     PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);
205:   }

207:   /* Loop over each field; drawing each in a different window */
208:   for (k=0; k<ndisplayfields; k++) {
210:     j = displayfields[k];

212:     /* determine the min and max value in plot */
213:     VecStrideMin(xin,j,NULL,&min);
214:     VecStrideMax(xin,j,NULL,&max);
215:     if (j < nbounds) {
216:       min = PetscMin(min,bounds[2*j]);
217:       max = PetscMax(max,bounds[2*j+1]);
218:     }
219:     if (min == max) {
220:       min -= 1.e-5;
221:       max += 1.e-5;
222:     }

224:     if (useports) {
225:       PetscDrawViewPortsSet(ports,k);
226:       DMDAGetFieldName(da,j,&tlabel);
227:     } else {
228:       const char *title;
229:       PetscViewerDrawGetHold(v,&hold);
230:       PetscViewerDrawGetDraw(v,k,&draw);
231:       PetscViewerDrawGetDrawAxis(v,k,&axis);
232:       DMDAGetFieldName(da,j,&title);
233:       if (title) PetscDrawSetTitle(draw,title);
234:       PetscDrawCheckResizedWindow(draw);
235:       if (!hold) PetscDrawClear(draw);
236:     }
237:     PetscDrawAxisSetLabels(axis,tlabel,xlabel,NULL);
238:     PetscDrawAxisSetLimits(axis,xmin,xmax,min,max);
239:     PetscDrawAxisDraw(axis);

241:     /* draw local part of vector */
242:     PetscObjectGetNewTag((PetscObject)xin,&tag);
243:     if (rank < size-1) { /*send value to right */
244:       MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag,comm);
245:       MPI_Send((void*)&array[j+(n-1)*dof],1,MPIU_REAL,rank+1,tag,comm);
246:     }
247:     if (rank) { /* receive value from left */
248:       MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag,comm,&status);
249:       MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag,comm,&status);
250:     }
251:     PetscDrawCollectiveBegin(draw);
252:     if (rank) {
253:       PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);
254:       if (showmarkers) PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);
255:     }
256:     for (i=1; i<n; i++) {
257:       PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+dof*i]),PETSC_DRAW_RED);
258:       if (showmarkers) PetscDrawMarker(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PETSC_DRAW_BLACK);
259:     }
260:     if (rank == size-1) {
261:       if (showmarkers) PetscDrawMarker(draw,PetscRealPart(xg[n-1]),PetscRealPart(array[j+dof*(n-1)]),PETSC_DRAW_BLACK);
262:     }
263:     PetscDrawCollectiveEnd(draw);
264:     PetscDrawFlush(draw);
265:     PetscDrawPause(draw);
266:     if (!useports) PetscDrawSave(draw);
267:   }
268:   if (useports) {
269:     PetscViewerDrawGetDraw(v,0,&draw);
270:     PetscDrawSave(draw);
271:   }

273:   PetscDrawViewPortsDestroy(ports);
274:   PetscFree(displayfields);
275:   VecRestoreArrayRead(xcoor,&xg);
276:   VecRestoreArrayRead(xin,&array);
277:   return 0;
278: }