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