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