Actual source code: gr1.c
petsc-3.5.4 2015-05-23
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 (values ignored for 1 dimensional problems)
19: - zmin,zmax - extremes in the z direction (values 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: PetscErrorCode DMDASelectFields(DM da,PetscInt *outfields,PetscInt **fields)
175: {
177: PetscInt step,ndisplayfields,*displayfields,k,j;
178: PetscBool flg;
181: DMDAGetInfo(da,0,0,0,0,0,0,0,&step,0,0,0,0,0);
182: PetscMalloc1(step,&displayfields);
183: for (k=0; k<step; k++) displayfields[k] = k;
184: ndisplayfields = step;
185: PetscOptionsGetIntArray(NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);
186: if (!ndisplayfields) ndisplayfields = step;
187: if (!flg) {
188: char **fields;
189: const char *fieldname;
190: PetscInt nfields = step;
191: PetscMalloc1(step,&fields);
192: PetscOptionsGetStringArray(NULL,"-draw_fields_by_name",fields,&nfields,&flg);
193: if (flg) {
194: ndisplayfields = 0;
195: for (k=0; k<nfields;k++) {
196: for (j=0; j<step; j++) {
197: DMDAGetFieldName(da,j,&fieldname);
198: PetscStrcmp(fieldname,fields[k],&flg);
199: if (flg) {
200: goto found;
201: }
202: }
203: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Unknown fieldname %s",fields[k]);
204: found: displayfields[ndisplayfields++] = j;
205: }
206: }
207: for (k=0; k<nfields; k++) {
208: PetscFree(fields[k]);
209: }
210: PetscFree(fields);
211: }
212: *fields = displayfields;
213: *outfields = ndisplayfields;
214: return(0);
215: }
217: #include <petscdraw.h>
220: PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v)
221: {
222: DM da;
223: PetscErrorCode ierr;
224: PetscMPIInt rank,size,tag1,tag2;
225: PetscInt i,n,N,step,istart,isize,j,nbounds;
226: MPI_Status status;
227: PetscReal coors[4],ymin,ymax,min,max,xmin = 0.0,xmax = 0.0,tmp = 0.0,xgtmp = 0.0;
228: const PetscScalar *array,*xg;
229: PetscDraw draw;
230: PetscBool isnull,showpoints = PETSC_FALSE;
231: MPI_Comm comm;
232: PetscDrawAxis axis;
233: Vec xcoor;
234: DMBoundaryType bx;
235: const PetscReal *bounds;
236: PetscInt *displayfields;
237: PetscInt k,ndisplayfields;
238: PetscBool hold;
241: PetscViewerDrawGetDraw(v,0,&draw);
242: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
243: PetscViewerDrawGetBounds(v,&nbounds,&bounds);
245: VecGetDM(xin,&da);
246: if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
248: PetscOptionsGetBool(NULL,"-draw_vec_mark_points",&showpoints,NULL);
250: DMDAGetInfo(da,0,&N,0,0,0,0,0,&step,0,&bx,0,0,0);
251: DMDAGetCorners(da,&istart,0,0,&isize,0,0);
252: VecGetArrayRead(xin,&array);
253: VecGetLocalSize(xin,&n);
254: n = n/step;
256: /* get coordinates of nodes */
257: DMGetCoordinates(da,&xcoor);
258: if (!xcoor) {
259: DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
260: DMGetCoordinates(da,&xcoor);
261: }
262: VecGetArrayRead(xcoor,&xg);
264: PetscObjectGetComm((PetscObject)xin,&comm);
265: MPI_Comm_size(comm,&size);
266: MPI_Comm_rank(comm,&rank);
268: /*
269: Determine the min and max x coordinate in plot
270: */
271: if (!rank) {
272: xmin = PetscRealPart(xg[0]);
273: }
274: if (rank == size-1) {
275: xmax = PetscRealPart(xg[n-1]);
276: }
277: MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
278: MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);
280: DMDASelectFields(da,&ndisplayfields,&displayfields);
281: for (k=0; k<ndisplayfields; k++) {
282: j = displayfields[k];
283: PetscViewerDrawGetDraw(v,k,&draw);
284: PetscDrawCheckResizedWindow(draw);
286: /*
287: Determine the min and max y coordinate in plot
288: */
289: min = 1.e20; max = -1.e20;
290: for (i=0; i<n; i++) {
291: if (PetscRealPart(array[j+i*step]) < min) min = PetscRealPart(array[j+i*step]);
292: if (PetscRealPart(array[j+i*step]) > max) max = PetscRealPart(array[j+i*step]);
293: }
294: if (min + 1.e-10 > max) {
295: min -= 1.e-5;
296: max += 1.e-5;
297: }
298: if (j < nbounds) {
299: min = PetscMin(min,bounds[2*j]);
300: max = PetscMax(max,bounds[2*j+1]);
301: }
303: MPI_Reduce(&min,&ymin,1,MPIU_REAL,MPIU_MIN,0,comm);
304: MPI_Reduce(&max,&ymax,1,MPIU_REAL,MPIU_MAX,0,comm);
306: PetscViewerDrawGetHold(v,&hold);
307: if (!hold) {
308: PetscDrawSynchronizedClear(draw);
309: }
310: PetscViewerDrawGetDrawAxis(v,k,&axis);
311: PetscLogObjectParent((PetscObject)draw,(PetscObject)axis);
312: if (!rank) {
313: const char *title;
315: PetscDrawAxisSetLimits(axis,xmin,xmax,ymin,ymax);
316: PetscDrawAxisDraw(axis);
317: PetscDrawGetCoordinates(draw,coors,coors+1,coors+2,coors+3);
318: DMDAGetFieldName(da,j,&title);
319: if (title) {PetscDrawSetTitle(draw,title);}
320: }
321: MPI_Bcast(coors,4,MPIU_REAL,0,comm);
322: if (rank) {
323: PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);
324: }
326: /* draw local part of vector */
327: PetscObjectGetNewTag((PetscObject)xin,&tag1);
328: PetscObjectGetNewTag((PetscObject)xin,&tag2);
329: if (rank < size-1) { /*send value to right */
330: MPI_Send((void*)&array[j+(n-1)*step],1,MPIU_REAL,rank+1,tag1,comm);
331: MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag1,comm);
332: }
333: if (!rank && bx == DM_BOUNDARY_PERIODIC && size > 1) { /* first processor sends first value to last */
334: MPI_Send((void*)&array[j],1,MPIU_REAL,size-1,tag2,comm);
335: }
337: for (i=1; i<n; i++) {
338: PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+step*i]),PETSC_DRAW_RED);
339: if (showpoints) {
340: PetscDrawPoint(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PETSC_DRAW_BLACK);
341: }
342: }
343: if (rank) { /* receive value from left */
344: MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag1,comm,&status);
345: MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag1,comm,&status);
346: PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);
347: if (showpoints) {
348: PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);
349: }
350: }
351: if (rank == size-1 && bx == DM_BOUNDARY_PERIODIC && size > 1) {
352: MPI_Recv(&tmp,1,MPIU_REAL,0,tag2,comm,&status);
353: /* 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 */
354: 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);
355: if (showpoints) {
356: PetscDrawPoint(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),PETSC_DRAW_BLACK);
357: }
358: }
359: PetscDrawSynchronizedFlush(draw);
360: PetscDrawPause(draw);
361: }
362: PetscFree(displayfields);
363: VecRestoreArrayRead(xcoor,&xg);
364: VecRestoreArrayRead(xin,&array);
365: return(0);
366: }