Actual source code: gr1.c
1: /*$Id: gr1.c,v 1.30 2001/08/10 03:34:34 bsmith Exp $*/
3: /*
4: Plots vectors obtained with DACreate1d()
5: */
7: #include petscda.h
9: #undef __FUNCT__
11: /*@
12: DASetUniformCoordinates - Sets a DA coordinates to be a uniform grid
14: Collective on DA
16: Input Parameters:
17: + da - the distributed array object
18: . xmin,xmax - extremes in the x direction
19: . xmin,xmax - extremes in the y direction
20: - xmin,xmax - extremes in the z direction
22: Level: beginner
24: .seealso: DASetCoordinates(), DAGetCoordinates(), DACreate1d(), DACreate2d(), DACreate3d()
26: @*/
27: int DASetUniformCoordinates(DA da,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
28: {
29: int i,j,k,ierr,M,N,P,istart,isize,jstart,jsize,kstart,ksize,dim,cnt;
30: PetscReal hx,hy,hz_;
31: Vec xcoor;
32: DAPeriodicType periodic;
33: PetscScalar *coors;
36: if (xmax <= xmin) SETERRQ2(1,"Xmax must be larger than xmin %g %g",xmin,xmax);
38: DAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&periodic,0);
39: DAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);
41: if (dim == 1) {
42: VecCreateMPI(PETSC_COMM_WORLD,isize,PETSC_DETERMINE,&xcoor);
43: if (periodic == DA_NONPERIODIC) hx = (xmax-xmin)/(M-1);
44: else hx = (xmax-xmin)/M;
45: VecGetArray(xcoor,&coors);
46: for (i=0; i<isize; i++) {
47: coors[i] = xmin + hx*(i+istart);
48: }
49: VecRestoreArray(xcoor,&coors);
50: } else if (dim == 2) {
51: if (ymax <= ymin) SETERRQ2(1,"Ymax must be larger than ymin %g %g",ymin,ymax);
52: VecCreateMPI(PETSC_COMM_WORLD,2*isize*jsize,PETSC_DETERMINE,&xcoor);
53: VecSetBlockSize(xcoor,2);
54: if (DAXPeriodic(periodic)) hx = (xmax-xmin)/(M);
55: else hx = (xmax-xmin)/(M-1);
56: if (DAYPeriodic(periodic)) hy = (ymax-ymin)/(N);
57: else hy = (ymax-ymin)/(N-1);
58: VecGetArray(xcoor,&coors);
59: cnt = 0;
60: for (j=0; j<jsize; j++) {
61: for (i=0; i<isize; i++) {
62: coors[cnt++] = xmin + hx*(i+istart);
63: coors[cnt++] = ymin + hy*(j+jstart);
64: }
65: }
66: VecRestoreArray(xcoor,&coors);
67: } else if (dim == 3) {
68: if (ymax <= ymin) SETERRQ2(1,"Ymax must be larger than ymin %g %g",ymin,ymax);
69: if (zmax <= zmin) SETERRQ2(1,"Zmax must be larger than zmin %g %g",zmin,zmax);
70: VecCreateMPI(PETSC_COMM_WORLD,3*isize*jsize*ksize,PETSC_DETERMINE,&xcoor);
71: VecSetBlockSize(xcoor,3);
72: if (DAXPeriodic(periodic)) hx = (xmax-xmin)/(M);
73: else hx = (xmax-xmin)/(M-1);
74: if (DAYPeriodic(periodic)) hy = (ymax-ymin)/(N);
75: else hy = (ymax-ymin)/(N-1);
76: if (DAZPeriodic(periodic)) hz_ = (zmax-zmin)/(P);
77: else hz_ = (zmax-zmin)/(P-1);
78: VecGetArray(xcoor,&coors);
79: cnt = 0;
80: for (k=0; k<ksize; k++) {
81: for (j=0; j<jsize; j++) {
82: for (i=0; i<isize; i++) {
83: coors[cnt++] = xmin + hx*(i+istart);
84: coors[cnt++] = ymin + hy*(j+jstart);
85: coors[cnt++] = zmin + hz_*(k+kstart);
86: }
87: }
88: }
89: VecRestoreArray(xcoor,&coors);
90: } else {
91: SETERRQ1(1,"Cannot create uniform coordinates for this dimension %dn",dim);
92: }
93: DASetCoordinates(da,xcoor);
94: PetscLogObjectParent(da,xcoor);
96: return(0);
97: }
99: #undef __FUNCT__
101: int VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v)
102: {
103: DA da;
104: int i,rank,size,ierr,n,tag1,tag2,N,step;
105: int istart,isize,j;
106: MPI_Status status;
107: PetscReal coors[4],ymin,ymax,min,max,xmin,xmax,tmp,xgtmp;
108: PetscScalar *array,*xg;
109: PetscDraw draw;
110: PetscTruth isnull,showpoints;
111: MPI_Comm comm;
112: PetscDrawAxis axis;
113: Vec xcoor;
114: DAPeriodicType periodic;
117: PetscViewerDrawGetDraw(v,0,&draw);
118: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
120: PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);
121: if (!da) SETERRQ(1,"Vector not generated from a DA");
123: PetscOptionsHasName(PETSC_NULL,"-draw_vec_mark_points",&showpoints);
125: DAGetInfo(da,0,&N,0,0,0,0,0,&step,0,&periodic,0);
126: DAGetCorners(da,&istart,0,0,&isize,0,0);
127: VecGetArray(xin,&array);
128: VecGetLocalSize(xin,&n);
129: n = n/step;
131: /* get coordinates of nodes */
132: DAGetCoordinates(da,&xcoor);
133: if (!xcoor) {
134: DASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
135: DAGetCoordinates(da,&xcoor);
136: }
137: VecGetArray(xcoor,&xg);
139: PetscObjectGetComm((PetscObject)xin,&comm);
140: MPI_Comm_size(comm,&size);
141: MPI_Comm_rank(comm,&rank);
143: /*
144: Determine the min and max x coordinate in plot
145: */
146: if (!rank) {
147: xmin = PetscRealPart(xg[0]);
148: }
149: if (rank == size-1) {
150: xmax = PetscRealPart(xg[n-1]);
151: }
152: MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
153: MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);
155: for (j=0; j<step; j++) {
156: PetscViewerDrawGetDraw(v,j,&draw);
157: PetscDrawCheckResizedWindow(draw);
159: /*
160: Determine the min and max y coordinate in plot
161: */
162: min = 1.e20; max = -1.e20;
163: for (i=0; i<n; i++) {
164: #if defined(PETSC_USE_COMPLEX)
165: if (PetscRealPart(array[j+i*step]) < min) min = PetscRealPart(array[j+i*step]);
166: if (PetscRealPart(array[j+i*step]) > max) max = PetscRealPart(array[j+i*step]);
167: #else
168: if (array[j+i*step] < min) min = array[j+i*step];
169: if (array[j+i*step] > max) max = array[j+i*step];
170: #endif
171: }
172: if (min + 1.e-10 > max) {
173: min -= 1.e-5;
174: max += 1.e-5;
175: }
176: MPI_Reduce(&min,&ymin,1,MPIU_REAL,MPI_MIN,0,comm);
177: MPI_Reduce(&max,&ymax,1,MPIU_REAL,MPI_MAX,0,comm);
179: PetscDrawSynchronizedClear(draw);
180: PetscViewerDrawGetDrawAxis(v,j,&axis);
181: PetscLogObjectParent(draw,axis);
182: if (!rank) {
183: char *title;
185: PetscDrawAxisSetLimits(axis,xmin,xmax,ymin,ymax);
186: PetscDrawAxisDraw(axis);
187: PetscDrawGetCoordinates(draw,coors,coors+1,coors+2,coors+3);
188: DAGetFieldName(da,j,&title);
189: if (title) {PetscDrawSetTitle(draw,title);}
190: }
191: MPI_Bcast(coors,4,MPIU_REAL,0,comm);
192: if (rank) {
193: PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);
194: }
196: /* draw local part of vector */
197: PetscObjectGetNewTag((PetscObject)xin,&tag1);
198: PetscObjectGetNewTag((PetscObject)xin,&tag2);
199: if (rank < size-1) { /*send value to right */
200: MPI_Send(&array[j+(n-1)*step],1,MPIU_REAL,rank+1,tag1,comm);
201: MPI_Send(&xg[n-1],1,MPIU_REAL,rank+1,tag1,comm);
202: }
203: if (!rank && periodic) { /* first processor sends first value to last */
204: MPI_Send(&array[j],1,MPIU_REAL,size-1,tag2,comm);
205: }
207: for (i=1; i<n; i++) {
208: #if !defined(PETSC_USE_COMPLEX)
209: PetscDrawLine(draw,xg[i-1],array[j+step*(i-1)],xg[i],array[j+step*i],
210: PETSC_DRAW_RED);
211: #else
212: PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),
213: PetscRealPart(xg[i]),PetscRealPart(array[j+step*i]),PETSC_DRAW_RED);
214: #endif
215: if (showpoints) {
216: PetscDrawPoint(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+step*(i-1)]),PETSC_DRAW_BLACK);
217: }
218: }
219: if (rank) { /* receive value from left */
220: MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag1,comm,&status);
221: MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag1,comm,&status);
222: #if !defined(PETSC_USE_COMPLEX)
223: PetscDrawLine(draw,xgtmp,tmp,xg[0],array[j],PETSC_DRAW_RED);
224: #else
225: PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),
226: PETSC_DRAW_RED);
227: #endif
228: if (showpoints) {
229: PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);
230: }
231: }
232: if (rank == size-1 && periodic) {
233: MPI_Recv(&tmp,1,MPIU_REAL,0,tag2,comm,&status);
234: #if !defined(PETSC_USE_COMPLEX)
235: PetscDrawLine(draw,xg[n-2],array[j+step*(n-1)],xg[n-1],tmp,PETSC_DRAW_RED);
236: #else
237: PetscDrawLine(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),
238: PetscRealPart(xg[n-1]),tmp,PETSC_DRAW_RED);
239: #endif
240: if (showpoints) {
241: PetscDrawPoint(draw,PetscRealPart(xg[n-2]),PetscRealPart(array[j+step*(n-1)]),PETSC_DRAW_BLACK);
242: }
243: }
244: PetscDrawSynchronizedFlush(draw);
245: PetscDrawPause(draw);
246: }
247: VecRestoreArray(xcoor,&xg);
248: VecRestoreArray(xin,&array);
249: return(0);
250: }