Actual source code: gr1.c

petsc-3.9.4 2018-09-11
Report Typos and Errors

  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 DMDA

 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()

 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;

 38:   DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&bx,&by,&bz,0);
 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:   PetscObjectGetComm((PetscObject)da,&comm);
 43:   DMGetDefaultSection(da,&section);
 44:   DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);
 45:   DMGetCoordinateDM(da, &cda);
 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};

 51:     DMDASetFieldName(cda, 0, "x");
 52:     if (dim > 1) {DMDASetFieldName(cda, 1, "y");}
 53:     if (dim > 2) {DMDASetFieldName(cda, 2, "z");}
 54:     DMDACreateSection(cda, numComp, numVertexDof, NULL, NULL);
 55:   }
 56:   DMCreateGlobalVector(cda, &xcoor);
 57:   if (section) {
 58:     PetscSection csection;
 59:     PetscInt     vStart, vEnd;

 61:     DMGetDefaultGlobalSection(cda,&csection);
 62:     VecGetArray(xcoor,&coors);
 63:     DMDAGetHeightStratum(da, dim, &vStart, &vEnd);
 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;

 75:         PetscSectionGetDof(csection, v, &dof);
 76:         PetscSectionGetOffset(csection, v, &off);
 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;

 87:           PetscSectionGetDof(csection, v, &dof);
 88:           PetscSectionGetOffset(csection, v, &off);
 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;

102:             PetscSectionGetDof(csection, v, &dof);
103:             PetscSectionGetOffset(csection, v, &off);
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:     VecRestoreArray(xcoor,&coors);
117:     DMSetCoordinates(da,xcoor);
118:     PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);
119:     VecDestroy(&xcoor);
120:     return(0);
121:   }
122:   if (dim == 1) {
123:     if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/M;
124:     else hx = (xmax-xmin)/(M-1);
125:     VecGetArray(xcoor,&coors);
126:     for (i=0; i<isize; i++) {
127:       coors[i] = xmin + hx*(i+istart);
128:     }
129:     VecRestoreArray(xcoor,&coors);
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:     VecGetArray(xcoor,&coors);
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:     VecRestoreArray(xcoor,&coors);
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:     VecGetArray(xcoor,&coors);
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:     VecRestoreArray(xcoor,&coors);
163:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim);
164:   DMSetCoordinates(da,xcoor);
165:   PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);
166:   VecDestroy(&xcoor);
167:   return(0);
168: }

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

180:   DMDAGetInfo(da,0,0,0,0,0,0,0,&step,0,0,0,0,0);
181:   PetscMalloc1(step,&displayfields);
182:   for (k=0; k<step; k++) displayfields[k] = k;
183:   ndisplayfields = step;
184:   PetscOptionsGetIntArray(NULL,NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);
185:   if (!ndisplayfields) ndisplayfields = step;
186:   if (!flg) {
187:     char       **fields;
188:     const char *fieldname;
189:     PetscInt   nfields = step;
190:     PetscMalloc1(step,&fields);
191:     PetscOptionsGetStringArray(NULL,NULL,"-draw_fields_by_name",fields,&nfields,&flg);
192:     if (flg) {
193:       ndisplayfields = 0;
194:       for (k=0; k<nfields;k++) {
195:         for (j=0; j<step; j++) {
196:           DMDAGetFieldName(da,j,&fieldname);
197:           PetscStrcmp(fieldname,fields[k],&flg);
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:       PetscFree(fields[k]);
208:     }
209:     PetscFree(fields);
210:   }
211:   *fields    = displayfields;
212:   *outfields = ndisplayfields;
213:   return(0);
214: }

216:  #include <petscdraw.h>

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;

242:   PetscViewerDrawGetDraw(v,0,&draw);
243:   PetscDrawIsNull(draw,&isnull);
244:   if (isnull) return(0);
245:   PetscViewerDrawGetBounds(v,&nbounds,&bounds);

247:   VecGetDM(xin,&da);
248:   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
249:   PetscObjectGetComm((PetscObject)xin,&comm);
250:   MPI_Comm_size(comm,&size);
251:   MPI_Comm_rank(comm,&rank);

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

255:   DMDAGetInfo(da,NULL,&N,NULL,NULL,NULL,NULL,NULL,&dof,NULL,&bx,NULL,NULL,NULL);
256:   DMDAGetCorners(da,&istart,NULL,NULL,&isize,NULL,NULL);
257:   VecGetArrayRead(xin,&array);
258:   VecGetLocalSize(xin,&n);
259:   n    = n/dof;

261:   /* Get coordinates of nodes */
262:   DMGetCoordinates(da,&xcoor);
263:   if (!xcoor) {
264:     DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
265:     DMGetCoordinates(da,&xcoor);
266:   }
267:   VecGetArrayRead(xcoor,&xg);
268:   DMDAGetCoordinateName(da,0,&xlabel);

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:   MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
274:   MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);

276:   DMDASelectFields(da,&ndisplayfields,&displayfields);
277:   PetscViewerGetFormat(v,&format);
278:   PetscOptionsGetBool(NULL,NULL,"-draw_ports",&useports,NULL);
279:   if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
280:   if (useports) {
281:     PetscViewerDrawGetDraw(v,0,&draw);
282:     PetscViewerDrawGetDrawAxis(v,0,&axis);
283:     PetscDrawCheckResizedWindow(draw);
284:     PetscDrawClear(draw);
285:     PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);
286:   }

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

292:     /* determine the min and max value in plot */
293:     VecStrideMin(xin,j,NULL,&min);
294:     VecStrideMax(xin,j,NULL,&max);
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:     }

304:     if (useports) {
305:       PetscDrawViewPortsSet(ports,k);
306:       DMDAGetFieldName(da,j,&tlabel);
307:     } else {
308:       const char *title;
309:       PetscViewerDrawGetHold(v,&hold);
310:       PetscViewerDrawGetDraw(v,k,&draw);
311:       PetscViewerDrawGetDrawAxis(v,k,&axis);
312:       DMDAGetFieldName(da,j,&title);
313:       if (title) {PetscDrawSetTitle(draw,title);}
314:       PetscDrawCheckResizedWindow(draw);
315:       if (!hold) {PetscDrawClear(draw);}
316:     }
317:     PetscDrawAxisSetLabels(axis,tlabel,xlabel,NULL);
318:     PetscDrawAxisSetLimits(axis,xmin,xmax,min,max);
319:     PetscDrawAxisDraw(axis);

321:     /* draw local part of vector */
322:     PetscObjectGetNewTag((PetscObject)xin,&tag);
323:     if (rank < size-1) { /*send value to right */
324:       MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag,comm);
325:       MPI_Send((void*)&array[j+(n-1)*dof],1,MPIU_REAL,rank+1,tag,comm);
326:     }
327:     if (rank) { /* receive value from left */
328:       MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag,comm,&status);
329:       MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag,comm,&status);
330:     }
331:     PetscDrawCollectiveBegin(draw);
332:     if (rank) {
333:       PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);
334:       if (showmarkers) {PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);}
335:     }
336:     for (i=1; i<n; i++) {
337:       PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+dof*i]),PETSC_DRAW_RED);
338:       if (showmarkers) {PetscDrawMarker(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PETSC_DRAW_BLACK);}
339:     }
340:     if (rank == size-1) {
341:       if (showmarkers) {PetscDrawMarker(draw,PetscRealPart(xg[n-1]),PetscRealPart(array[j+dof*(n-1)]),PETSC_DRAW_BLACK);}
342:     }
343:     PetscDrawCollectiveEnd(draw);
344:     PetscDrawFlush(draw);
345:     PetscDrawPause(draw);
346:     if (!useports) {PetscDrawSave(draw);}
347:   }
348:   if (useports) {
349:     PetscViewerDrawGetDraw(v,0,&draw);
350:     PetscDrawSave(draw);
351:   }

353:   PetscDrawViewPortsDestroy(ports);
354:   PetscFree(displayfields);
355:   VecRestoreArrayRead(xcoor,&xg);
356:   VecRestoreArrayRead(xin,&array);
357:   return(0);
358: }