Actual source code: gr1.c
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;
37: DMDAGetInfo(da,&dim,&M,&N,&P,NULL,NULL,NULL,NULL,NULL,&bx,&by,&bz,NULL);
41: PetscObjectGetComm((PetscObject)da,&comm);
42: DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);
43: DMGetCoordinateDM(da, &cda);
44: DMCreateGlobalVector(cda, &xcoor);
45: if (dim == 1) {
46: if (bx == DM_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 (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
55: else hx = (xmax-xmin)/(M-1);
56: if (by == DM_BOUNDARY_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 (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
69: else hx = (xmax-xmin)/(M-1);
70: if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N);
71: else hy = (ymax-ymin)/(N-1);
72: if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P);
73: else hz_ = (zmax-zmin)/(P-1);
74: VecGetArray(xcoor,&coors);
75: cnt = 0;
76: for (k=0; k<ksize; k++) {
77: for (j=0; j<jsize; j++) {
78: for (i=0; i<isize; i++) {
79: coors[cnt++] = xmin + hx*(i+istart);
80: coors[cnt++] = ymin + hy*(j+jstart);
81: coors[cnt++] = zmin + hz_*(k+kstart);
82: }
83: }
84: }
85: VecRestoreArray(xcoor,&coors);
86: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D",dim);
87: DMSetCoordinates(da,xcoor);
88: PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);
89: VecDestroy(&xcoor);
90: return 0;
91: }
93: /*
94: Allows a user to select a subset of the fields to be drawn by VecView() when the vector comes from a DMDA
95: */
96: PetscErrorCode DMDASelectFields(DM da,PetscInt *outfields,PetscInt **fields)
97: {
98: PetscInt step,ndisplayfields,*displayfields,k,j;
99: PetscBool flg;
101: DMDAGetInfo(da,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&step,NULL,NULL,NULL,NULL,NULL);
102: PetscMalloc1(step,&displayfields);
103: for (k=0; k<step; k++) displayfields[k] = k;
104: ndisplayfields = step;
105: PetscOptionsGetIntArray(NULL,NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);
106: if (!ndisplayfields) ndisplayfields = step;
107: if (!flg) {
108: char **fields;
109: const char *fieldname;
110: PetscInt nfields = step;
111: PetscMalloc1(step,&fields);
112: PetscOptionsGetStringArray(NULL,NULL,"-draw_fields_by_name",fields,&nfields,&flg);
113: if (flg) {
114: ndisplayfields = 0;
115: for (k=0; k<nfields;k++) {
116: for (j=0; j<step; j++) {
117: DMDAGetFieldName(da,j,&fieldname);
118: PetscStrcmp(fieldname,fields[k],&flg);
119: if (flg) {
120: goto found;
121: }
122: }
123: SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Unknown fieldname %s",fields[k]);
124: found: displayfields[ndisplayfields++] = j;
125: }
126: }
127: for (k=0; k<nfields; k++) {
128: PetscFree(fields[k]);
129: }
130: PetscFree(fields);
131: }
132: *fields = displayfields;
133: *outfields = ndisplayfields;
134: return 0;
135: }
137: #include <petscdraw.h>
139: PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v)
140: {
141: DM da;
142: PetscMPIInt rank,size,tag;
143: PetscInt i,n,N,dof,istart,isize,j,nbounds;
144: MPI_Status status;
145: PetscReal min,max,xmin = 0.0,xmax = 0.0,tmp = 0.0,xgtmp = 0.0;
146: const PetscScalar *array,*xg;
147: PetscDraw draw;
148: PetscBool isnull,useports = PETSC_FALSE,showmarkers = PETSC_FALSE;
149: MPI_Comm comm;
150: PetscDrawAxis axis;
151: Vec xcoor;
152: DMBoundaryType bx;
153: const char *tlabel = NULL,*xlabel = NULL;
154: const PetscReal *bounds;
155: PetscInt *displayfields;
156: PetscInt k,ndisplayfields;
157: PetscBool hold;
158: PetscDrawViewPorts *ports = NULL;
159: PetscViewerFormat format;
161: PetscViewerDrawGetDraw(v,0,&draw);
162: PetscDrawIsNull(draw,&isnull);
163: if (isnull) return 0;
164: PetscViewerDrawGetBounds(v,&nbounds,&bounds);
166: VecGetDM(xin,&da);
168: PetscObjectGetComm((PetscObject)xin,&comm);
169: MPI_Comm_size(comm,&size);
170: MPI_Comm_rank(comm,&rank);
172: PetscOptionsGetBool(NULL,NULL,"-draw_vec_use_markers",&showmarkers,NULL);
174: DMDAGetInfo(da,NULL,&N,NULL,NULL,NULL,NULL,NULL,&dof,NULL,&bx,NULL,NULL,NULL);
175: DMDAGetCorners(da,&istart,NULL,NULL,&isize,NULL,NULL);
176: VecGetArrayRead(xin,&array);
177: VecGetLocalSize(xin,&n);
178: n = n/dof;
180: /* Get coordinates of nodes */
181: DMGetCoordinates(da,&xcoor);
182: if (!xcoor) {
183: DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
184: DMGetCoordinates(da,&xcoor);
185: }
186: VecGetArrayRead(xcoor,&xg);
187: DMDAGetCoordinateName(da,0,&xlabel);
189: /* Determine the min and max coordinate in plot */
190: if (rank == 0) xmin = PetscRealPart(xg[0]);
191: if (rank == size-1) xmax = PetscRealPart(xg[n-1]);
192: MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
193: MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);
195: DMDASelectFields(da,&ndisplayfields,&displayfields);
196: PetscViewerGetFormat(v,&format);
197: PetscOptionsGetBool(NULL,NULL,"-draw_ports",&useports,NULL);
198: if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
199: if (useports) {
200: PetscViewerDrawGetDraw(v,0,&draw);
201: PetscViewerDrawGetDrawAxis(v,0,&axis);
202: PetscDrawCheckResizedWindow(draw);
203: PetscDrawClear(draw);
204: PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);
205: }
207: /* Loop over each field; drawing each in a different window */
208: for (k=0; k<ndisplayfields; k++) {
210: j = displayfields[k];
212: /* determine the min and max value in plot */
213: VecStrideMin(xin,j,NULL,&min);
214: VecStrideMax(xin,j,NULL,&max);
215: if (j < nbounds) {
216: min = PetscMin(min,bounds[2*j]);
217: max = PetscMax(max,bounds[2*j+1]);
218: }
219: if (min == max) {
220: min -= 1.e-5;
221: max += 1.e-5;
222: }
224: if (useports) {
225: PetscDrawViewPortsSet(ports,k);
226: DMDAGetFieldName(da,j,&tlabel);
227: } else {
228: const char *title;
229: PetscViewerDrawGetHold(v,&hold);
230: PetscViewerDrawGetDraw(v,k,&draw);
231: PetscViewerDrawGetDrawAxis(v,k,&axis);
232: DMDAGetFieldName(da,j,&title);
233: if (title) PetscDrawSetTitle(draw,title);
234: PetscDrawCheckResizedWindow(draw);
235: if (!hold) PetscDrawClear(draw);
236: }
237: PetscDrawAxisSetLabels(axis,tlabel,xlabel,NULL);
238: PetscDrawAxisSetLimits(axis,xmin,xmax,min,max);
239: PetscDrawAxisDraw(axis);
241: /* draw local part of vector */
242: PetscObjectGetNewTag((PetscObject)xin,&tag);
243: if (rank < size-1) { /*send value to right */
244: MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag,comm);
245: MPI_Send((void*)&array[j+(n-1)*dof],1,MPIU_REAL,rank+1,tag,comm);
246: }
247: if (rank) { /* receive value from left */
248: MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag,comm,&status);
249: MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag,comm,&status);
250: }
251: PetscDrawCollectiveBegin(draw);
252: if (rank) {
253: PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);
254: if (showmarkers) PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);
255: }
256: for (i=1; i<n; i++) {
257: PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+dof*i]),PETSC_DRAW_RED);
258: if (showmarkers) PetscDrawMarker(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PETSC_DRAW_BLACK);
259: }
260: if (rank == size-1) {
261: if (showmarkers) PetscDrawMarker(draw,PetscRealPart(xg[n-1]),PetscRealPart(array[j+dof*(n-1)]),PETSC_DRAW_BLACK);
262: }
263: PetscDrawCollectiveEnd(draw);
264: PetscDrawFlush(draw);
265: PetscDrawPause(draw);
266: if (!useports) PetscDrawSave(draw);
267: }
268: if (useports) {
269: PetscViewerDrawGetDraw(v,0,&draw);
270: PetscDrawSave(draw);
271: }
273: PetscDrawViewPortsDestroy(ports);
274: PetscFree(displayfields);
275: VecRestoreArrayRead(xcoor,&xg);
276: VecRestoreArrayRead(xin,&array);
277: return 0;
278: }