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