Actual source code: gr1.c
petsc-3.10.5 2019-03-28
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 DMDA
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()
23: @*/
24: PetscErrorCode DMDASetUniformCoordinates(DM da,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
25: {
26: MPI_Comm comm;
27: PetscSection section;
28: DM cda;
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: DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&bx,&by,&bz,0);
39: if (xmax < xmin) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %g %g",(double)xmin,(double)xmax);
40: 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);
41: 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);
42: PetscObjectGetComm((PetscObject)da,&comm);
43: DMGetSection(da,§ion);
44: DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);
45: DMGetCoordinateDM(da, &cda);
46: if (section) {
47: /* This would be better as a vector, but this is compatible */
48: PetscInt numComp[3] = {1, 1, 1};
49: PetscInt numVertexDof[3] = {1, 1, 1};
51: DMDASetFieldName(cda, 0, "x");
52: if (dim > 1) {DMDASetFieldName(cda, 1, "y");}
53: if (dim > 2) {DMDASetFieldName(cda, 2, "z");}
54: DMDACreateSection(cda, numComp, numVertexDof, NULL, NULL);
55: }
56: DMCreateGlobalVector(cda, &xcoor);
57: if (section) {
58: PetscSection csection;
59: PetscInt vStart, vEnd;
61: DMGetGlobalSection(cda,&csection);
62: VecGetArray(xcoor,&coors);
63: DMDAGetHeightStratum(da, dim, &vStart, &vEnd);
64: if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M+1);
65: else hx = (xmax-xmin)/(M ? M : 1);
66: if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N+1);
67: else hy = (ymax-ymin)/(N ? N : 1);
68: if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P+1);
69: else hz_ = (zmax-zmin)/(P ? P : 1);
70: switch (dim) {
71: case 1:
72: for (i = 0; i < isize+1; ++i) {
73: PetscInt v = i+vStart, dof, off;
75: PetscSectionGetDof(csection, v, &dof);
76: PetscSectionGetOffset(csection, v, &off);
77: if (off >= 0) {
78: coors[off] = xmin + hx*(i+istart);
79: }
80: }
81: break;
82: case 2:
83: for (j = 0; j < jsize+1; ++j) {
84: for (i = 0; i < isize+1; ++i) {
85: PetscInt v = j*(isize+1)+i+vStart, dof, off;
87: PetscSectionGetDof(csection, v, &dof);
88: PetscSectionGetOffset(csection, v, &off);
89: if (off >= 0) {
90: coors[off+0] = xmin + hx*(i+istart);
91: coors[off+1] = ymin + hy*(j+jstart);
92: }
93: }
94: }
95: break;
96: case 3:
97: for (k = 0; k < ksize+1; ++k) {
98: for (j = 0; j < jsize+1; ++j) {
99: for (i = 0; i < isize+1; ++i) {
100: PetscInt v = (k*(jsize+1)+j)*(isize+1)+i+vStart, dof, off;
102: PetscSectionGetDof(csection, v, &dof);
103: PetscSectionGetOffset(csection, v, &off);
104: if (off >= 0) {
105: coors[off+0] = xmin + hx*(i+istart);
106: coors[off+1] = ymin + hy*(j+jstart);
107: coors[off+2] = zmin + hz_*(k+kstart);
108: }
109: }
110: }
111: }
112: break;
113: default:
114: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim);
115: }
116: VecRestoreArray(xcoor,&coors);
117: DMSetCoordinates(da,xcoor);
118: PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);
119: VecDestroy(&xcoor);
120: return(0);
121: }
122: if (dim == 1) {
123: if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/M;
124: else hx = (xmax-xmin)/(M-1);
125: VecGetArray(xcoor,&coors);
126: for (i=0; i<isize; i++) {
127: coors[i] = xmin + hx*(i+istart);
128: }
129: VecRestoreArray(xcoor,&coors);
130: } else if (dim == 2) {
131: if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
132: else hx = (xmax-xmin)/(M-1);
133: if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N);
134: else hy = (ymax-ymin)/(N-1);
135: VecGetArray(xcoor,&coors);
136: cnt = 0;
137: for (j=0; j<jsize; j++) {
138: for (i=0; i<isize; i++) {
139: coors[cnt++] = xmin + hx*(i+istart);
140: coors[cnt++] = ymin + hy*(j+jstart);
141: }
142: }
143: VecRestoreArray(xcoor,&coors);
144: } else if (dim == 3) {
145: if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax-xmin)/(M);
146: else hx = (xmax-xmin)/(M-1);
147: if (by == DM_BOUNDARY_PERIODIC) hy = (ymax-ymin)/(N);
148: else hy = (ymax-ymin)/(N-1);
149: if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P);
150: else hz_ = (zmax-zmin)/(P-1);
151: VecGetArray(xcoor,&coors);
152: cnt = 0;
153: for (k=0; k<ksize; k++) {
154: for (j=0; j<jsize; j++) {
155: for (i=0; i<isize; i++) {
156: coors[cnt++] = xmin + hx*(i+istart);
157: coors[cnt++] = ymin + hy*(j+jstart);
158: coors[cnt++] = zmin + hz_*(k+kstart);
159: }
160: }
161: }
162: VecRestoreArray(xcoor,&coors);
163: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot create uniform coordinates for this dimension %D\n",dim);
164: DMSetCoordinates(da,xcoor);
165: PetscLogObjectParent((PetscObject)da,(PetscObject)xcoor);
166: VecDestroy(&xcoor);
167: return(0);
168: }
170: /*
171: Allows a user to select a subset of the fields to be drawn by VecView() when the vector comes from a DMDA
172: */
173: PetscErrorCode DMDASelectFields(DM da,PetscInt *outfields,PetscInt **fields)
174: {
176: PetscInt step,ndisplayfields,*displayfields,k,j;
177: PetscBool flg;
180: DMDAGetInfo(da,0,0,0,0,0,0,0,&step,0,0,0,0,0);
181: PetscMalloc1(step,&displayfields);
182: for (k=0; k<step; k++) displayfields[k] = k;
183: ndisplayfields = step;
184: PetscOptionsGetIntArray(NULL,NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);
185: if (!ndisplayfields) ndisplayfields = step;
186: if (!flg) {
187: char **fields;
188: const char *fieldname;
189: PetscInt nfields = step;
190: PetscMalloc1(step,&fields);
191: PetscOptionsGetStringArray(NULL,NULL,"-draw_fields_by_name",fields,&nfields,&flg);
192: if (flg) {
193: ndisplayfields = 0;
194: for (k=0; k<nfields;k++) {
195: for (j=0; j<step; j++) {
196: DMDAGetFieldName(da,j,&fieldname);
197: PetscStrcmp(fieldname,fields[k],&flg);
198: if (flg) {
199: goto found;
200: }
201: }
202: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Unknown fieldname %s",fields[k]);
203: found: displayfields[ndisplayfields++] = j;
204: }
205: }
206: for (k=0; k<nfields; k++) {
207: PetscFree(fields[k]);
208: }
209: PetscFree(fields);
210: }
211: *fields = displayfields;
212: *outfields = ndisplayfields;
213: return(0);
214: }
216: #include <petscdraw.h>
218: PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v)
219: {
220: DM da;
221: PetscErrorCode ierr;
222: PetscMPIInt rank,size,tag;
223: PetscInt i,n,N,dof,istart,isize,j,nbounds;
224: MPI_Status status;
225: PetscReal min,max,xmin = 0.0,xmax = 0.0,tmp = 0.0,xgtmp = 0.0;
226: const PetscScalar *array,*xg;
227: PetscDraw draw;
228: PetscBool isnull,useports = PETSC_FALSE,showmarkers = PETSC_FALSE;
229: MPI_Comm comm;
230: PetscDrawAxis axis;
231: Vec xcoor;
232: DMBoundaryType bx;
233: const char *tlabel = NULL,*xlabel = NULL;
234: const PetscReal *bounds;
235: PetscInt *displayfields;
236: PetscInt k,ndisplayfields;
237: PetscBool hold;
238: PetscDrawViewPorts *ports = NULL;
239: PetscViewerFormat format;
242: PetscViewerDrawGetDraw(v,0,&draw);
243: PetscDrawIsNull(draw,&isnull);
244: if (isnull) return(0);
245: PetscViewerDrawGetBounds(v,&nbounds,&bounds);
247: VecGetDM(xin,&da);
248: if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
249: PetscObjectGetComm((PetscObject)xin,&comm);
250: MPI_Comm_size(comm,&size);
251: MPI_Comm_rank(comm,&rank);
253: PetscOptionsGetBool(NULL,NULL,"-draw_vec_use_markers",&showmarkers,NULL);
255: DMDAGetInfo(da,NULL,&N,NULL,NULL,NULL,NULL,NULL,&dof,NULL,&bx,NULL,NULL,NULL);
256: DMDAGetCorners(da,&istart,NULL,NULL,&isize,NULL,NULL);
257: VecGetArrayRead(xin,&array);
258: VecGetLocalSize(xin,&n);
259: n = n/dof;
261: /* Get coordinates of nodes */
262: DMGetCoordinates(da,&xcoor);
263: if (!xcoor) {
264: DMDASetUniformCoordinates(da,0.0,1.0,0.0,0.0,0.0,0.0);
265: DMGetCoordinates(da,&xcoor);
266: }
267: VecGetArrayRead(xcoor,&xg);
268: DMDAGetCoordinateName(da,0,&xlabel);
270: /* Determine the min and max coordinate in plot */
271: if (!rank) xmin = PetscRealPart(xg[0]);
272: if (rank == size-1) xmax = PetscRealPart(xg[n-1]);
273: MPI_Bcast(&xmin,1,MPIU_REAL,0,comm);
274: MPI_Bcast(&xmax,1,MPIU_REAL,size-1,comm);
276: DMDASelectFields(da,&ndisplayfields,&displayfields);
277: PetscViewerGetFormat(v,&format);
278: PetscOptionsGetBool(NULL,NULL,"-draw_ports",&useports,NULL);
279: if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
280: if (useports) {
281: PetscViewerDrawGetDraw(v,0,&draw);
282: PetscViewerDrawGetDrawAxis(v,0,&axis);
283: PetscDrawCheckResizedWindow(draw);
284: PetscDrawClear(draw);
285: PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);
286: }
288: /* Loop over each field; drawing each in a different window */
289: for (k=0; k<ndisplayfields; k++) {
290: j = displayfields[k];
292: /* determine the min and max value in plot */
293: VecStrideMin(xin,j,NULL,&min);
294: VecStrideMax(xin,j,NULL,&max);
295: if (j < nbounds) {
296: min = PetscMin(min,bounds[2*j]);
297: max = PetscMax(max,bounds[2*j+1]);
298: }
299: if (min == max) {
300: min -= 1.e-5;
301: max += 1.e-5;
302: }
304: if (useports) {
305: PetscDrawViewPortsSet(ports,k);
306: DMDAGetFieldName(da,j,&tlabel);
307: } else {
308: const char *title;
309: PetscViewerDrawGetHold(v,&hold);
310: PetscViewerDrawGetDraw(v,k,&draw);
311: PetscViewerDrawGetDrawAxis(v,k,&axis);
312: DMDAGetFieldName(da,j,&title);
313: if (title) {PetscDrawSetTitle(draw,title);}
314: PetscDrawCheckResizedWindow(draw);
315: if (!hold) {PetscDrawClear(draw);}
316: }
317: PetscDrawAxisSetLabels(axis,tlabel,xlabel,NULL);
318: PetscDrawAxisSetLimits(axis,xmin,xmax,min,max);
319: PetscDrawAxisDraw(axis);
321: /* draw local part of vector */
322: PetscObjectGetNewTag((PetscObject)xin,&tag);
323: if (rank < size-1) { /*send value to right */
324: MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag,comm);
325: MPI_Send((void*)&array[j+(n-1)*dof],1,MPIU_REAL,rank+1,tag,comm);
326: }
327: if (rank) { /* receive value from left */
328: MPI_Recv(&xgtmp,1,MPIU_REAL,rank-1,tag,comm,&status);
329: MPI_Recv(&tmp,1,MPIU_REAL,rank-1,tag,comm,&status);
330: }
331: PetscDrawCollectiveBegin(draw);
332: if (rank) {
333: PetscDrawLine(draw,xgtmp,tmp,PetscRealPart(xg[0]),PetscRealPart(array[j]),PETSC_DRAW_RED);
334: if (showmarkers) {PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);}
335: }
336: for (i=1; i<n; i++) {
337: PetscDrawLine(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PetscRealPart(xg[i]),PetscRealPart(array[j+dof*i]),PETSC_DRAW_RED);
338: if (showmarkers) {PetscDrawMarker(draw,PetscRealPart(xg[i-1]),PetscRealPart(array[j+dof*(i-1)]),PETSC_DRAW_BLACK);}
339: }
340: if (rank == size-1) {
341: if (showmarkers) {PetscDrawMarker(draw,PetscRealPart(xg[n-1]),PetscRealPart(array[j+dof*(n-1)]),PETSC_DRAW_BLACK);}
342: }
343: PetscDrawCollectiveEnd(draw);
344: PetscDrawFlush(draw);
345: PetscDrawPause(draw);
346: if (!useports) {PetscDrawSave(draw);}
347: }
348: if (useports) {
349: PetscViewerDrawGetDraw(v,0,&draw);
350: PetscDrawSave(draw);
351: }
353: PetscDrawViewPortsDestroy(ports);
354: PetscFree(displayfields);
355: VecRestoreArrayRead(xcoor,&xg);
356: VecRestoreArrayRead(xin,&array);
357: return(0);
358: }