Actual source code: gr1.c
petsc-3.12.5 2020-03-29
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: DMBoundaryType bx,by,bz;
29: Vec xcoor;
30: PetscScalar *coors;
31: PetscReal hx,hy,hz_;
32: PetscInt i,j,k,M,N,P,istart,isize,jstart,jsize,kstart,ksize,dim,cnt;
33: PetscErrorCode ierr;
37: DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&bx,&by,&bz,0);
38: if (xmax < xmin) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %g %g",(double)xmin,(double)xmax);
39: if ((dim > 1) && (ymax < ymin)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %g %g",(double)ymin,(double)ymax);
40: if ((dim > 2) && (zmax < zmin)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"zmax must be larger than zmin %g %g",(double)zmin,(double)zmax);
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 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",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: {
99: PetscInt step,ndisplayfields,*displayfields,k,j;
100: PetscBool flg;
103: DMDAGetInfo(da,0,0,0,0,0,0,0,&step,0,0,0,0,0);
104: PetscMalloc1(step,&displayfields);
105: for (k=0; k<step; k++) displayfields[k] = k;
106: ndisplayfields = step;
107: PetscOptionsGetIntArray(NULL,NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);
108: if (!ndisplayfields) ndisplayfields = step;
109: if (!flg) {
110: char **fields;
111: const char *fieldname;
112: PetscInt nfields = step;
113: PetscMalloc1(step,&fields);
114: PetscOptionsGetStringArray(NULL,NULL,"-draw_fields_by_name",fields,&nfields,&flg);
115: if (flg) {
116: ndisplayfields = 0;
117: for (k=0; k<nfields;k++) {
118: for (j=0; j<step; j++) {
119: DMDAGetFieldName(da,j,&fieldname);
120: PetscStrcmp(fieldname,fields[k],&flg);
121: if (flg) {
122: goto found;
123: }
124: }
125: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Unknown fieldname %s",fields[k]);
126: found: displayfields[ndisplayfields++] = j;
127: }
128: }
129: for (k=0; k<nfields; k++) {
130: PetscFree(fields[k]);
131: }
132: PetscFree(fields);
133: }
134: *fields = displayfields;
135: *outfields = ndisplayfields;
136: return(0);
137: }
139: #include <petscdraw.h>
141: PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v)
142: {
143: DM da;
144: PetscErrorCode ierr;
145: PetscMPIInt rank,size,tag;
146: PetscInt i,n,N,dof,istart,isize,j,nbounds;
147: MPI_Status status;
148: PetscReal min,max,xmin = 0.0,xmax = 0.0,tmp = 0.0,xgtmp = 0.0;
149: const PetscScalar *array,*xg;
150: PetscDraw draw;
151: PetscBool isnull,useports = PETSC_FALSE,showmarkers = PETSC_FALSE;
152: MPI_Comm comm;
153: PetscDrawAxis axis;
154: Vec xcoor;
155: DMBoundaryType bx;
156: const char *tlabel = NULL,*xlabel = NULL;
157: const PetscReal *bounds;
158: PetscInt *displayfields;
159: PetscInt k,ndisplayfields;
160: PetscBool hold;
161: PetscDrawViewPorts *ports = NULL;
162: PetscViewerFormat format;
165: PetscViewerDrawGetDraw(v,0,&draw);
166: PetscDrawIsNull(draw,&isnull);
167: if (isnull) return(0);
168: PetscViewerDrawGetBounds(v,&nbounds,&bounds);
170: VecGetDM(xin,&da);
171: if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
172: PetscObjectGetComm((PetscObject)xin,&comm);
173: MPI_Comm_size(comm,&size);
174: MPI_Comm_rank(comm,&rank);
176: PetscOptionsGetBool(NULL,NULL,"-draw_vec_use_markers",&showmarkers,NULL);
178: DMDAGetInfo(da,NULL,&N,NULL,NULL,NULL,NULL,NULL,&dof,NULL,&bx,NULL,NULL,NULL);
179: DMDAGetCorners(da,&istart,NULL,NULL,&isize,NULL,NULL);
180: VecGetArrayRead(xin,&array);
181: VecGetLocalSize(xin,&n);
182: n = n/dof;
184: /* Get coordinates of nodes */
185: DMGetCoordinates(da,&xcoor);
186: if (!xcoor) {
187: DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
188: DMGetCoordinates(da,&xcoor);
189: }
190: VecGetArrayRead(xcoor,&xg);
191: DMDAGetCoordinateName(da,0,&xlabel);
193: /* Determine the min and max coordinate in plot */
194: if (!rank) xmin = PetscRealPart(xg[0]);
195: if (rank == size-1) xmax = PetscRealPart(xg[n-1]);
196: MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
197: MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);
199: DMDASelectFields(da,&ndisplayfields,&displayfields);
200: PetscViewerGetFormat(v,&format);
201: PetscOptionsGetBool(NULL,NULL,"-draw_ports",&useports,NULL);
202: if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
203: if (useports) {
204: PetscViewerDrawGetDraw(v,0,&draw);
205: PetscViewerDrawGetDrawAxis(v,0,&axis);
206: PetscDrawCheckResizedWindow(draw);
207: PetscDrawClear(draw);
208: PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);
209: }
211: /* Loop over each field; drawing each in a different window */
212: for (k=0; k<ndisplayfields; k++) {
213: j = displayfields[k];
215: /* determine the min and max value in plot */
216: VecStrideMin(xin,j,NULL,&min);
217: VecStrideMax(xin,j,NULL,&max);
218: if (j < nbounds) {
219: min = PetscMin(min,bounds[2*j]);
220: max = PetscMax(max,bounds[2*j+1]);
221: }
222: if (min == max) {
223: min -= 1.e-5;
224: max += 1.e-5;
225: }
227: if (useports) {
228: PetscDrawViewPortsSet(ports,k);
229: DMDAGetFieldName(da,j,&tlabel);
230: } else {
231: const char *title;
232: PetscViewerDrawGetHold(v,&hold);
233: PetscViewerDrawGetDraw(v,k,&draw);
234: PetscViewerDrawGetDrawAxis(v,k,&axis);
235: DMDAGetFieldName(da,j,&title);
236: if (title) {PetscDrawSetTitle(draw,title);}
237: PetscDrawCheckResizedWindow(draw);
238: if (!hold) {PetscDrawClear(draw);}
239: }
240: PetscDrawAxisSetLabels(axis,tlabel,xlabel,NULL);
241: PetscDrawAxisSetLimits(axis,xmin,xmax,min,max);
242: PetscDrawAxisDraw(axis);
244: /* draw local part of vector */
245: PetscObjectGetNewTag((PetscObject)xin,&tag);
246: if (rank < size-1) { /*send value to right */
247: MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag,comm);
248: MPI_Send((void*)&array[j+(n-1)*dof],1,MPIU_REAL,rank+1,tag,comm);
249: }
250: if (rank) { /* receive value from left */
251: MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag,comm,&status);
252: MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag,comm,&status);
253: }
254: PetscDrawCollectiveBegin(draw);
255: if (rank) {
256: PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);
257: if (showmarkers) {PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);}
258: }
259: for (i=1; i<n; i++) {
260: PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+dof*i]),PETSC_DRAW_RED);
261: if (showmarkers) {PetscDrawMarker(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PETSC_DRAW_BLACK);}
262: }
263: if (rank == size-1) {
264: if (showmarkers) {PetscDrawMarker(draw,PetscRealPart(xg[n-1]),PetscRealPart(array[j+dof*(n-1)]),PETSC_DRAW_BLACK);}
265: }
266: PetscDrawCollectiveEnd(draw);
267: PetscDrawFlush(draw);
268: PetscDrawPause(draw);
269: if (!useports) {PetscDrawSave(draw);}
270: }
271: if (useports) {
272: PetscViewerDrawGetDraw(v,0,&draw);
273: PetscDrawSave(draw);
274: }
276: PetscDrawViewPortsDestroy(ports);
277: PetscFree(displayfields);
278: VecRestoreArrayRead(xcoor,&xg);
279: VecRestoreArrayRead(xin,&array);
280: return(0);
281: }