Actual source code: gr1.c
petsc-3.6.1 2015-08-06
2: /*
3: Plots vectors obtained with DMDACreate1d()
4: */
6: #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
10: /*@
11: DMDASetUniformCoordinates - Sets a DMDA coordinates to be a uniform grid
13: Collective on DMDA
15: Input Parameters:
16: + da - the distributed array object
17: . xmin,xmax - extremes in the x direction
18: . ymin,ymax - extremes in the y direction (value ignored for 1 dimensional problems)
19: - zmin,zmax - extremes in the z direction (value ignored for 1 or 2 dimensional problems)
21: Level: beginner
23: .seealso: DMSetCoordinates(), DMGetCoordinates(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
25: @*/
26: PetscErrorCode DMDASetUniformCoordinates(DM da,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
27: {
28: MPI_Comm comm;
29: PetscSection section;
30: DM cda;
31: DMBoundaryType bx,by,bz;
32: Vec xcoor;
33: PetscScalar *coors;
34: PetscReal hx,hy,hz_;
35: PetscInt i,j,k,M,N,P,istart,isize,jstart,jsize,kstart,ksize,dim,cnt;
36: PetscErrorCode ierr;
40: DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&bx,&by,&bz,0);
41: if (xmax < xmin) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %g %g",(double)xmin,(double)xmax);
42: 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);
43: 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);
44: PetscObjectGetComm((PetscObject)da,&comm);
45: DMGetDefaultSection(da,§ion);
46: DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);
47: DMGetCoordinateDM(da, &cda);
48: if (section) {
49: /* This would be better as a vector, but this is compatible */
50: PetscInt numComp[3] = {1, 1, 1};
51: PetscInt numVertexDof[3] = {1, 1, 1};
53: DMDASetFieldName(cda, 0, "x");
54: if (dim > 1) {DMDASetFieldName(cda, 1, "y");}
55: if (dim > 2) {DMDASetFieldName(cda, 2, "z");}
56: DMDACreateSection(cda, numComp, numVertexDof, NULL, NULL);
57: }
58: DMCreateGlobalVector(cda, &xcoor);
59: if (section) {
60: PetscSection csection;
61: PetscInt vStart, vEnd;
63: DMGetDefaultGlobalSection(cda,&csection);
64: VecGetArray(xcoor,&coors);
65: DMDAGetHeightStratum(da, dim, &vStart, &vEnd);
66: if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M+1);
67: else hx = (xmax-xmin)/(M ? M : 1);
68: if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N+1);
69: else hy = (ymax-ymin)/(N ? N : 1);
70: if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P+1);
71: else hz_ = (zmax-zmin)/(P ? P : 1);
72: switch (dim) {
73: case 1:
74: for (i = 0; i < isize+1; ++i) {
75: PetscInt v = i+vStart, dof, off;
77: PetscSectionGetDof(csection, v, &dof);
78: PetscSectionGetOffset(csection, v, &off);
79: if (off >= 0) {
80: coors[off] = xmin + hx*(i+istart);
81: }
82: }
83: break;
84: case 2:
85: for (j = 0; j < jsize+1; ++j) {
86: for (i = 0; i < isize+1; ++i) {
87: PetscInt v = j*(isize+1)+i+vStart, dof, off;
89: PetscSectionGetDof(csection, v, &dof);
90: PetscSectionGetOffset(csection, v, &off);
91: if (off >= 0) {
92: coors[off+0] = xmin + hx*(i+istart);
93: coors[off+1] = ymin + hy*(j+jstart);
94: }
95: }
96: }
97: break;
98: case 3:
99: for (k = 0; k < ksize+1; ++k) {
100: for (j = 0; j < jsize+1; ++j) {
101: for (i = 0; i < isize+1; ++i) {
102: PetscInt v = (k*(jsize+1)+j)*(isize+1)+i+vStart, dof, off;
104: PetscSectionGetDof(csection, v, &dof);
105: PetscSectionGetOffset(csection, v, &off);
106: if (off >= 0) {
107: coors[off+0] = xmin + hx*(i+istart);
108: coors[off+1] = ymin + hy*(j+jstart);
109: coors[off+2] = zmin + hz_*(k+kstart);
110: }
111: }
112: }
113: }
114: break;
115: default:
116: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim);
117: }
118: VecRestoreArray(xcoor,&coors);
119: DMSetCoordinates(da,xcoor);
120: PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);
121: VecDestroy(&xcoor);
122: return(0);
123: }
124: if (dim == 1) {
125: if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/M;
126: else hx = (xmax-xmin)/(M-1);
127: VecGetArray(xcoor,&coors);
128: for (i=0; i<isize; i++) {
129: coors[i] = xmin + hx*(i+istart);
130: }
131: VecRestoreArray(xcoor,&coors);
132: } else if (dim == 2) {
133: if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
134: else hx = (xmax-xmin)/(M-1);
135: if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N);
136: else hy = (ymax-ymin)/(N-1);
137: VecGetArray(xcoor,&coors);
138: cnt = 0;
139: for (j=0; j<jsize; j++) {
140: for (i=0; i<isize; i++) {
141: coors[cnt++] = xmin + hx*(i+istart);
142: coors[cnt++] = ymin + hy*(j+jstart);
143: }
144: }
145: VecRestoreArray(xcoor,&coors);
146: } else if (dim == 3) {
147: if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
148: else hx = (xmax-xmin)/(M-1);
149: if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N);
150: else hy = (ymax-ymin)/(N-1);
151: if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P);
152: else hz_ = (zmax-zmin)/(P-1);
153: VecGetArray(xcoor,&coors);
154: cnt = 0;
155: for (k=0; k<ksize; k++) {
156: for (j=0; j<jsize; j++) {
157: for (i=0; i<isize; i++) {
158: coors[cnt++] = xmin + hx*(i+istart);
159: coors[cnt++] = ymin + hy*(j+jstart);
160: coors[cnt++] = zmin + hz_*(k+kstart);
161: }
162: }
163: }
164: VecRestoreArray(xcoor,&coors);
165: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim);
166: DMSetCoordinates(da,xcoor);
167: PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);
168: VecDestroy(&xcoor);
169: return(0);
170: }
174: /*
175: Allows a user to select a subset of the fields to be drawn by VecView() when the vector comes from a DMDA
176: */
177: PetscErrorCode DMDASelectFields(DM da,PetscInt *outfields,PetscInt **fields)
178: {
180: PetscInt step,ndisplayfields,*displayfields,k,j;
181: PetscBool flg;
184: DMDAGetInfo(da,0,0,0,0,0,0,0,&step,0,0,0,0,0);
185: PetscMalloc1(step,&displayfields);
186: for (k=0; k<step; k++) displayfields[k] = k;
187: ndisplayfields = step;
188: PetscOptionsGetIntArray(NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);
189: if (!ndisplayfields) ndisplayfields = step;
190: if (!flg) {
191: char **fields;
192: const char *fieldname;
193: PetscInt nfields = step;
194: PetscMalloc1(step,&fields);
195: PetscOptionsGetStringArray(NULL,"-draw_fields_by_name",fields,&nfields,&flg);
196: if (flg) {
197: ndisplayfields = 0;
198: for (k=0; k<nfields;k++) {
199: for (j=0; j<step; j++) {
200: DMDAGetFieldName(da,j,&fieldname);
201: PetscStrcmp(fieldname,fields[k],&flg);
202: if (flg) {
203: goto found;
204: }
205: }
206: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Unknown fieldname %s",fields[k]);
207: found: displayfields[ndisplayfields++] = j;
208: }
209: }
210: for (k=0; k<nfields; k++) {
211: PetscFree(fields[k]);
212: }
213: PetscFree(fields);
214: }
215: *fields = displayfields;
216: *outfields = ndisplayfields;
217: return(0);
218: }
220: #include <petscdraw.h>
221: #if defined(PETSC_HAVE_SETJMP_H) && defined(PETSC_HAVE_X)
222: #include <X11/Xlib.h>
223: #include <X11/Xutil.h>
224: #include <setjmp.h>
225: static jmp_buf PetscXIOErrorJumpBuf;
226: static void PetscXIOHandler(Display *dpy)
227: {
228: longjmp(PetscXIOErrorJumpBuf, 1);
229: }
230: #endif
234: PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v)
235: {
236: DM da;
237: PetscErrorCode ierr;
238: PetscMPIInt rank,size,tag1,tag2;
239: PetscInt i,n,N,step,istart,isize,j,nbounds;
240: MPI_Status status;
241: PetscReal coors[4],ymin,ymax,min,max,xmin = 0.0,xmax = 0.0,tmp = 0.0,xgtmp = 0.0;
242: const PetscScalar *array,*xg;
243: PetscDraw draw;
244: PetscBool isnull,showmarkers = PETSC_FALSE;
245: MPI_Comm comm;
246: PetscDrawAxis axis;
247: Vec xcoor;
248: DMBoundaryType bx;
249: const PetscReal *bounds;
250: PetscInt *displayfields;
251: PetscInt k,ndisplayfields;
252: PetscBool hold;
255: PetscViewerDrawGetDraw(v,0,&draw);
256: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
257: PetscViewerDrawGetBounds(v,&nbounds,&bounds);
259: VecGetDM(xin,&da);
260: if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
262: PetscOptionsGetBool(NULL,"-draw_vec_use_markers",&showmarkers,NULL);
264: DMDAGetInfo(da,0,&N,0,0,0,0,0,&step,0,&bx,0,0,0);
265: DMDAGetCorners(da,&istart,0,0,&isize,0,0);
266: VecGetArrayRead(xin,&array);
267: VecGetLocalSize(xin,&n);
268: n = n/step;
270: /* get coordinates of nodes */
271: DMGetCoordinates(da,&xcoor);
272: if (!xcoor) {
273: DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
274: DMGetCoordinates(da,&xcoor);
275: }
276: VecGetArrayRead(xcoor,&xg);
278: PetscObjectGetComm((PetscObject)xin,&comm);
279: MPI_Comm_size(comm,&size);
280: MPI_Comm_rank(comm,&rank);
282: /*
283: Determine the min and max x coordinate in plot
284: */
285: if (!rank) {
286: xmin = PetscRealPart(xg[0]);
287: }
288: if (rank == size-1) {
289: xmax = PetscRealPart(xg[n-1]);
290: }
291: MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
292: MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);
294: DMDASelectFields(da,&ndisplayfields,&displayfields);
295: #if defined(PETSC_HAVE_SETJMP_H) && defined(PETSC_HAVE_X)
296: if (!setjmp(PetscXIOErrorJumpBuf)) XSetIOErrorHandler((XIOErrorHandler)PetscXIOHandler);
297: else {
298: XSetIOErrorHandler(NULL);
299: PetscDrawSetType(draw,PETSC_DRAW_NULL);
300: return(0);
301: }
302: #endif
303: for (k=0; k<ndisplayfields; k++) {
304: j = displayfields[k];
305: PetscViewerDrawGetDraw(v,k,&draw);
306: PetscDrawCheckResizedWindow(draw);
308: /*
309: Determine the min and max y coordinate in plot
310: */
311: min = 1.e20; max = -1.e20;
312: for (i=0; i<n; i++) {
313: if (PetscRealPart(array[j+i*step]) < min) min = PetscRealPart(array[j+i*step]);
314: if (PetscRealPart(array[j+i*step]) > max) max = PetscRealPart(array[j+i*step]);
315: }
316: if (min + 1.e-10 > max) {
317: min -= 1.e-5;
318: max += 1.e-5;
319: }
320: if (j < nbounds) {
321: min = PetscMin(min,bounds[2*j]);
322: max = PetscMax(max,bounds[2*j+1]);
323: }
325: MPI_Reduce(&min,&ymin,1,MPIU_REAL,MPIU_MIN,0,comm);
326: MPI_Reduce(&max,&ymax,1,MPIU_REAL,MPIU_MAX,0,comm);
328: PetscViewerDrawGetHold(v,&hold);
329: if (!hold) {
330: PetscDrawSynchronizedClear(draw);
331: }
332: PetscViewerDrawGetDrawAxis(v,k,&axis);
333: PetscLogObjectParent((PetscObject)draw,(PetscObject)axis);
334: if (!rank) {
335: const char *title;
337: PetscDrawAxisSetLimits(axis,xmin,xmax,ymin,ymax);
338: PetscDrawAxisDraw(axis);
339: PetscDrawGetCoordinates(draw,coors,coors+1,coors+2,coors+3);
340: DMDAGetFieldName(da,j,&title);
341: if (title) {PetscDrawSetTitle(draw,title);}
342: }
343: MPI_Bcast(coors,4,MPIU_REAL,0,comm);
344: if (rank) {
345: PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);
346: }
348: /* draw local part of vector */
349: PetscObjectGetNewTag((PetscObject)xin,&tag1);
350: PetscObjectGetNewTag((PetscObject)xin,&tag2);
351: if (rank < size-1) { /*send value to right */
352: MPI_Send((void*)&array[j+(n-1)*step],1,MPIU_REAL,rank+1,tag1,comm);
353: MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag1,comm);
354: }
355: if (!rank && bx == DM_BOUNDARY_PERIODIC && size > 1) { /* first processor sends first value to last */
356: MPI_Send((void*)&array[j],1,MPIU_REAL,size-1,tag2,comm);
357: }
359: for (i=1; i<n; i++) {
360: PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+step*i]),PETSC_DRAW_RED);
361: if (showmarkers) {
362: PetscDrawMarker(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PETSC_DRAW_BLACK);
363: }
364: }
365: if (rank) { /* receive value from left */
366: MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag1,comm,&status);
367: MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag1,comm,&status);
368: PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);
369: if (showmarkers) {
370: PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);
371: }
372: }
373: if (rank == size-1 && bx == DM_BOUNDARY_PERIODIC && size > 1) {
374: MPI_Recv(&tmp,1,MPIU_REAL,0,tag2,comm,&status);
375: /* If the mesh is not uniform we do not know the mesh spacing between the last point on the right and the first ghost point */
376: PetscDrawLine(draw,PetscRealPart(xg[n-1]),PetscRealPart(array[j+step*(n-1)]),PetscRealPart(xg[n-1]+(xg[n-1]-xg[n-2])),tmp,PETSC_DRAW_RED);
377: if (showmarkers) {
378: PetscDrawMarker(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),PETSC_DRAW_BLACK);
379: }
380: }
381: PetscDrawSynchronizedFlush(draw);
382: PetscDrawPause(draw);
383: }
384: #if defined(PETSC_HAVE_SETJMP_H) && defined(PETSC_HAVE_X)
385: XSetIOErrorHandler(NULL);
386: #endif
387: PetscFree(displayfields);
388: VecRestoreArrayRead(xcoor,&xg);
389: VecRestoreArrayRead(xin,&array);
390: return(0);
391: }