Actual source code: gr1.c
petsc-3.3-p7 2013-05-11
2: /*
3: Plots vectors obtained with DMDACreate1d()
4: */
6: #include <petscdmda.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 (use PETSC_NULL for 1 dimensional problems)
19: - zmin,zmax - extremes in the z direction (use PETSC_NULL for 1 or 2 dimensional problems)
21: Level: beginner
23: .seealso: DMDASetCoordinates(), DMDAGetCoordinates(), 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: DM cda;
30: DMDABoundaryType bx,by,bz;
31: Vec xcoor;
32: PetscScalar *coors;
33: PetscReal hx,hy,hz_;
34: PetscInt i,j,k,M,N,P,istart,isize,jstart,jsize,kstart,ksize,dim,cnt;
35: PetscErrorCode ierr;
38: if (xmax <= xmin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %G %G",xmin,xmax);
40: PetscObjectGetComm((PetscObject)da,&comm);
41: DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&bx,&by,&bz,0);
42: DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);
43: DMDAGetCoordinateDA(da, &cda);
44: DMCreateGlobalVector(cda, &xcoor);
45: if (dim == 1) {
46: if (bx == DMDA_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 (ymax <= ymin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %G %G",ymin,ymax);
55: if (bx == DMDA_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
56: else hx = (xmax-xmin)/(M-1);
57: if (by == DMDA_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N);
58: else hy = (ymax-ymin)/(N-1);
59: VecGetArray(xcoor,&coors);
60: cnt = 0;
61: for (j=0; j<jsize; j++) {
62: for (i=0; i<isize; i++) {
63: coors[cnt++] = xmin + hx*(i+istart);
64: coors[cnt++] = ymin + hy*(j+jstart);
65: }
66: }
67: VecRestoreArray(xcoor,&coors);
68: } else if (dim == 3) {
69: if (ymax <= ymin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %G %G",ymin,ymax);
70: if (zmax <= zmin) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"zmax must be larger than zmin %G %G",zmin,zmax);
71: if (bx == DMDA_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
72: else hx = (xmax-xmin)/(M-1);
73: if (by == DMDA_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N);
74: else hy = (ymax-ymin)/(N-1);
75: if (bz == DMDA_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P);
76: else hz_ = (zmax-zmin)/(P-1);
77: VecGetArray(xcoor,&coors);
78: cnt = 0;
79: for (k=0; k<ksize; k++) {
80: for (j=0; j<jsize; j++) {
81: for (i=0; i<isize; i++) {
82: coors[cnt++] = xmin + hx*(i+istart);
83: coors[cnt++] = ymin + hy*(j+jstart);
84: coors[cnt++] = zmin + hz_*(k+kstart);
85: }
86: }
87: }
88: VecRestoreArray(xcoor,&coors);
89: } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim);
90: DMDASetCoordinates(da,xcoor);
91: PetscLogObjectParent(da,xcoor);
92: VecDestroy(&xcoor);
93: return(0);
94: }
98: PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v)
99: {
100: DM da;
101: PetscErrorCode ierr;
102: PetscMPIInt rank,size,tag1,tag2;
103: PetscInt i,n,N,step,istart,isize,j,nbounds;
104: MPI_Status status;
105: PetscReal coors[4],ymin,ymax,min,max,xmin,xmax,tmp,xgtmp;
106: const PetscScalar *array,*xg;
107: PetscDraw draw;
108: PetscBool isnull,showpoints = PETSC_FALSE;
109: MPI_Comm comm;
110: PetscDrawAxis axis;
111: Vec xcoor;
112: DMDABoundaryType bx;
113: const PetscReal *bounds;
114: PetscInt *displayfields;
115: PetscInt k,ndisplayfields;
116: PetscBool flg,hold;
119: PetscViewerDrawGetDraw(v,0,&draw);
120: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
121: PetscViewerDrawGetBounds(v,&nbounds,&bounds);
123: PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);
124: if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
126: PetscOptionsGetBool(PETSC_NULL,"-draw_vec_mark_points",&showpoints,PETSC_NULL);
128: DMDAGetInfo(da,0,&N,0,0,0,0,0,&step,0,&bx,0,0,0);
129: DMDAGetCorners(da,&istart,0,0,&isize,0,0);
130: VecGetArrayRead(xin,&array);
131: VecGetLocalSize(xin,&n);
132: n = n/step;
134: /* get coordinates of nodes */
135: DMDAGetCoordinates(da,&xcoor);
136: if (!xcoor) {
137: DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
138: DMDAGetCoordinates(da,&xcoor);
139: }
140: VecGetArrayRead(xcoor,&xg);
142: PetscObjectGetComm((PetscObject)xin,&comm);
143: MPI_Comm_size(comm,&size);
144: MPI_Comm_rank(comm,&rank);
146: /*
147: Determine the min and max x coordinate in plot
148: */
149: if (!rank) {
150: xmin = PetscRealPart(xg[0]);
151: }
152: if (rank == size-1) {
153: xmax = PetscRealPart(xg[n-1]);
154: }
155: MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
156: MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);
158: PetscMalloc(step*sizeof(PetscInt),&displayfields);
159: for (i=0; i<step; i++) displayfields[i] = i;
160: ndisplayfields = step;
161: PetscOptionsGetIntArray(PETSC_NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);
162: if (!flg) ndisplayfields = step;
163: for (k=0; k<ndisplayfields; k++) {
164: j = displayfields[k];
165: PetscViewerDrawGetDraw(v,k,&draw);
166: PetscDrawCheckResizedWindow(draw);
168: /*
169: Determine the min and max y coordinate in plot
170: */
171: min = 1.e20; max = -1.e20;
172: for (i=0; i<n; i++) {
173: if (PetscRealPart(array[j+i*step]) < min) min = PetscRealPart(array[j+i*step]);
174: if (PetscRealPart(array[j+i*step]) > max) max = PetscRealPart(array[j+i*step]);
175: }
176: if (min + 1.e-10 > max) {
177: min -= 1.e-5;
178: max += 1.e-5;
179: }
180: if (j < nbounds) {
181: min = PetscMin(min,bounds[2*j]);
182: max = PetscMax(max,bounds[2*j+1]);
183: }
185: MPI_Reduce(&min,&ymin,1,MPIU_REAL,MPIU_MIN,0,comm);
186: MPI_Reduce(&max,&ymax,1,MPIU_REAL,MPIU_MAX,0,comm);
188: PetscViewerDrawGetHold(v,&hold);
189: if (!hold) {
190: PetscDrawSynchronizedClear(draw);
191: }
192: PetscViewerDrawGetDrawAxis(v,k,&axis);
193: PetscLogObjectParent(draw,axis);
194: if (!rank) {
195: const char *title;
197: PetscDrawAxisSetLimits(axis,xmin,xmax,ymin,ymax);
198: PetscDrawAxisDraw(axis);
199: PetscDrawGetCoordinates(draw,coors,coors+1,coors+2,coors+3);
200: DMDAGetFieldName(da,j,&title);
201: if (title) {PetscDrawSetTitle(draw,title);}
202: }
203: MPI_Bcast(coors,4,MPIU_REAL,0,comm);
204: if (rank) {
205: PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);
206: }
208: /* draw local part of vector */
209: PetscObjectGetNewTag((PetscObject)xin,&tag1);
210: PetscObjectGetNewTag((PetscObject)xin,&tag2);
211: if (rank < size-1) { /*send value to right */
212: MPI_Send((void*)&array[j+(n-1)*step],1,MPIU_REAL,rank+1,tag1,comm);
213: MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag1,comm);
214: }
215: if (!rank && bx == DMDA_BOUNDARY_PERIODIC && size > 1) { /* first processor sends first value to last */
216: MPI_Send((void*)&array[j],1,MPIU_REAL,size-1,tag2,comm);
217: }
219: for (i=1; i<n; i++) {
220: PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+step*i]),PETSC_DRAW_RED);
221: if (showpoints) {
222: PetscDrawPoint(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PETSC_DRAW_BLACK);
223: }
224: }
225: if (rank) { /* receive value from left */
226: MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag1,comm,&status);
227: MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag1,comm,&status);
228: PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);
229: if (showpoints) {
230: PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);
231: }
232: }
233: if (rank == size-1 && bx == DMDA_BOUNDARY_PERIODIC && size > 1) {
234: MPI_Recv(&tmp,1,MPIU_REAL,0,tag2,comm,&status);
235: PetscDrawLine(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),PetscRealPart(xg[n-1]),tmp,PETSC_DRAW_RED);
236: if (showpoints) {
237: PetscDrawPoint(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),PETSC_DRAW_BLACK);
238: }
239: }
240: PetscDrawSynchronizedFlush(draw);
241: PetscDrawPause(draw);
242: }
243: PetscFree(displayfields);
244: VecRestoreArrayRead(xcoor,&xg);
245: VecRestoreArrayRead(xin,&array);
246: return(0);
247: }