Actual source code: gr1.c

petsc-3.14.6 2021-03-30
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 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;
 34:   PetscErrorCode   ierr;

 38:   if (!dd->gtol) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set coordinates until after DMDA has been setup");
 39:   DMDAGetInfo(da,&dim,&M,&N,&P,NULL,NULL,NULL,NULL,NULL,&bx,&by,&bz,NULL);
 40:   if (xmax < xmin) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %g %g",(double)xmin,(double)xmax);
 41:   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);
 42:   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);
 43:   PetscObjectGetComm((PetscObject)da,&comm);
 44:   DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);
 45:   DMGetCoordinateDM(da, &cda);
 46:   DMCreateGlobalVector(cda, &xcoor);
 47:   if (dim == 1) {
 48:     if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/M;
 49:     else hx = (xmax-xmin)/(M-1);
 50:     VecGetArray(xcoor,&coors);
 51:     for (i=0; i<isize; i++) {
 52:       coors[i] = xmin + hx*(i+istart);
 53:     }
 54:     VecRestoreArray(xcoor,&coors);
 55:   } else if (dim == 2) {
 56:     if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
 57:     else hx = (xmax-xmin)/(M-1);
 58:     if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N);
 59:     else hy = (ymax-ymin)/(N-1);
 60:     VecGetArray(xcoor,&coors);
 61:     cnt  = 0;
 62:     for (j=0; j<jsize; j++) {
 63:       for (i=0; i<isize; i++) {
 64:         coors[cnt++] = xmin + hx*(i+istart);
 65:         coors[cnt++] = ymin + hy*(j+jstart);
 66:       }
 67:     }
 68:     VecRestoreArray(xcoor,&coors);
 69:   } else if (dim == 3) {
 70:     if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
 71:     else hx = (xmax-xmin)/(M-1);
 72:     if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N);
 73:     else hy = (ymax-ymin)/(N-1);
 74:     if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P);
 75:     else hz_ = (zmax-zmin)/(P-1);
 76:     VecGetArray(xcoor,&coors);
 77:     cnt  = 0;
 78:     for (k=0; k<ksize; k++) {
 79:       for (j=0; j<jsize; j++) {
 80:         for (i=0; i<isize; i++) {
 81:           coors[cnt++] = xmin + hx*(i+istart);
 82:           coors[cnt++] = ymin + hy*(j+jstart);
 83:           coors[cnt++] = zmin + hz_*(k+kstart);
 84:         }
 85:       }
 86:     }
 87:     VecRestoreArray(xcoor,&coors);
 88:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim);
 89:   DMSetCoordinates(da,xcoor);
 90:   PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);
 91:   VecDestroy(&xcoor);
 92:   return(0);
 93: }

 95: /*
 96:     Allows a user to select a subset of the fields to be drawn by VecView() when the vector comes from a DMDA
 97: */
 98: PetscErrorCode DMDASelectFields(DM da,PetscInt *outfields,PetscInt **fields)
 99: {
101:   PetscInt       step,ndisplayfields,*displayfields,k,j;
102:   PetscBool      flg;

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

141: #include <petscdraw.h>

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

167:   PetscViewerDrawGetDraw(v,0,&draw);
168:   PetscDrawIsNull(draw,&isnull);
169:   if (isnull) return(0);
170:   PetscViewerDrawGetBounds(v,&nbounds,&bounds);

172:   VecGetDM(xin,&da);
173:   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
174:   PetscObjectGetComm((PetscObject)xin,&comm);
175:   MPI_Comm_size(comm,&size);
176:   MPI_Comm_rank(comm,&rank);

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

180:   DMDAGetInfo(da,NULL,&N,NULL,NULL,NULL,NULL,NULL,&dof,NULL,&bx,NULL,NULL,NULL);
181:   DMDAGetCorners(da,&istart,NULL,NULL,&isize,NULL,NULL);
182:   VecGetArrayRead(xin,&array);
183:   VecGetLocalSize(xin,&n);
184:   n    = n/dof;

186:   /* Get coordinates of nodes */
187:   DMGetCoordinates(da,&xcoor);
188:   if (!xcoor) {
189:     DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
190:     DMGetCoordinates(da,&xcoor);
191:   }
192:   VecGetArrayRead(xcoor,&xg);
193:   DMDAGetCoordinateName(da,0,&xlabel);

195:   /* Determine the min and max coordinate in plot */
196:   if (!rank) xmin = PetscRealPart(xg[0]);
197:   if (rank == size-1) xmax = PetscRealPart(xg[n-1]);
198:   MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
199:   MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);

201:   DMDASelectFields(da,&ndisplayfields,&displayfields);
202:   PetscViewerGetFormat(v,&format);
203:   PetscOptionsGetBool(NULL,NULL,"-draw_ports",&useports,NULL);
204:   if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
205:   if (useports) {
206:     PetscViewerDrawGetDraw(v,0,&draw);
207:     PetscViewerDrawGetDrawAxis(v,0,&axis);
208:     PetscDrawCheckResizedWindow(draw);
209:     PetscDrawClear(draw);
210:     PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);
211:   }

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

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

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

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

278:   PetscDrawViewPortsDestroy(ports);
279:   PetscFree(displayfields);
280:   VecRestoreArrayRead(xcoor,&xg);
281:   VecRestoreArrayRead(xin,&array);
282:   return(0);
283: }