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