Actual source code: gr1.c
1: /*
2: Plots vectors obtained with DMDACreate1d()
3: */
5: #include <petsc/private/dmdaimpl.h>
7: /*@
8: DMDASetUniformCoordinates - Sets a `DMDA` coordinates to be a uniform grid
10: Collective
12: Input Parameters:
13: + da - the `DMDA` object
14: . xmin - min extreme in the x direction
15: . xmax - max extreme in the x direction
16: . ymin - min extreme in the y direction (value ignored for 1 dimensional problems)
17: . ymax - max extreme in the y direction (value ignored for 1 dimensional problems)
18: . zmin - min extreme in the z direction (value ignored for 1 or 2 dimensional problems)
19: - zmax - max extreme in the z direction (value ignored for 1 or 2 dimensional problems)
21: Level: beginner
23: .seealso: [](sec_struct), `DM`, `DMDA`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMStagSetUniformCoordinates()`
24: @*/
25: PetscErrorCode DMDASetUniformCoordinates(DM da, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax, PetscReal zmin, PetscReal zmax)
26: {
27: MPI_Comm comm;
28: DM cda;
29: DM_DA *dd = (DM_DA *)da->data;
30: DMBoundaryType bx, by, bz;
31: Vec xcoor;
32: PetscScalar *coors;
33: PetscReal hx = 0., hy = 0., hz_ = 0.;
34: PetscInt i, j, k, M, N, P, istart, isize, jstart, jsize, kstart, ksize, dim, cnt;
36: PetscFunctionBegin;
38: PetscCheck(dd->gtol, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "Cannot set coordinates until after DMDA has been setup");
39: PetscCall(DMDAGetInfo(da, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, &bx, &by, &bz, NULL));
40: PetscCheck(xmax >= xmin, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_INCOMP, "xmax must be larger than xmin %g %g", (double)xmin, (double)xmax);
41: PetscCheck(!(dim > 1) || !(ymax < ymin), PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_INCOMP, "ymax must be larger than ymin %g %g", (double)ymin, (double)ymax);
42: PetscCheck(!(dim > 2) || !(zmax < zmin), PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_INCOMP, "zmax must be larger than zmin %g %g", (double)zmin, (double)zmax);
43: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
44: PetscCall(DMDAGetCorners(da, &istart, &jstart, &kstart, &isize, &jsize, &ksize));
45: PetscCall(DMGetCoordinateDM(da, &cda));
46: PetscCall(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: PetscCall(VecGetArray(xcoor, &coors));
51: for (i = 0; i < isize; i++) coors[i] = xmin + hx * (i + istart);
52: PetscCall(VecRestoreArray(xcoor, &coors));
53: } else if (dim == 2) {
54: if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax - xmin) / M;
55: else hx = (xmax - xmin) / (M - 1);
56: if (by == DM_BOUNDARY_PERIODIC) hy = (ymax - ymin) / N;
57: else hy = (ymax - ymin) / (N - 1);
58: PetscCall(VecGetArray(xcoor, &coors));
59: cnt = 0;
60: for (j = 0; j < jsize; j++) {
61: for (i = 0; i < isize; i++) {
62: coors[cnt++] = xmin + hx * (i + istart);
63: coors[cnt++] = ymin + hy * (j + jstart);
64: }
65: }
66: PetscCall(VecRestoreArray(xcoor, &coors));
67: } else if (dim == 3) {
68: if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax - xmin) / M;
69: else hx = (xmax - xmin) / (M - 1);
70: if (by == DM_BOUNDARY_PERIODIC) hy = (ymax - ymin) / N;
71: else hy = (ymax - ymin) / (N - 1);
72: if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax - zmin) / P;
73: else hz_ = (zmax - zmin) / (P - 1);
74: PetscCall(VecGetArray(xcoor, &coors));
75: cnt = 0;
76: for (k = 0; k < ksize; k++) {
77: for (j = 0; j < jsize; j++) {
78: for (i = 0; i < isize; i++) {
79: coors[cnt++] = xmin + hx * (i + istart);
80: coors[cnt++] = ymin + hy * (j + jstart);
81: coors[cnt++] = zmin + hz_ * (k + kstart);
82: }
83: }
84: }
85: PetscCall(VecRestoreArray(xcoor, &coors));
86: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Cannot create uniform coordinates for this dimension %" PetscInt_FMT, dim);
87: PetscCall(DMSetCoordinates(da, xcoor));
88: PetscCall(VecDestroy(&xcoor));
89: // Handle periodicity
90: if (bx == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_PERIODIC || bz == DM_BOUNDARY_PERIODIC) {
91: PetscReal maxCell[3] = {hx, hy, hz_};
92: PetscReal Lstart[3] = {xmin, ymin, zmin};
93: PetscReal L[3] = {xmax - xmin, ymax - ymin, zmax - zmin};
95: PetscCall(DMSetPeriodicity(da, maxCell, Lstart, L));
96: }
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: /*
101: Allows a user to select a subset of the fields to be drawn by VecView() when the vector comes from a DMDA
102: */
103: PetscErrorCode DMDASelectFields(DM da, PetscInt *outfields, PetscInt **fields)
104: {
105: PetscInt step, ndisplayfields, *displayfields, k, j;
106: PetscBool flg;
108: PetscFunctionBegin;
109: PetscCall(DMDAGetInfo(da, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &step, NULL, NULL, NULL, NULL, NULL));
110: PetscCall(PetscMalloc1(step, &displayfields));
111: for (k = 0; k < step; k++) displayfields[k] = k;
112: ndisplayfields = step;
113: PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-draw_fields", displayfields, &ndisplayfields, &flg));
114: if (!ndisplayfields) ndisplayfields = step;
115: if (!flg) {
116: char **fields;
117: const char *fieldname;
118: PetscInt nfields = step;
119: PetscCall(PetscMalloc1(step, &fields));
120: PetscCall(PetscOptionsGetStringArray(NULL, NULL, "-draw_fields_by_name", fields, &nfields, &flg));
121: if (flg) {
122: ndisplayfields = 0;
123: for (k = 0; k < nfields; k++) {
124: for (j = 0; j < step; j++) {
125: PetscCall(DMDAGetFieldName(da, j, &fieldname));
126: PetscCall(PetscStrcmp(fieldname, fields[k], &flg));
127: if (flg) goto found;
128: }
129: SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Unknown fieldname %s", fields[k]);
130: found:
131: displayfields[ndisplayfields++] = j;
132: }
133: }
134: for (k = 0; k < nfields; k++) PetscCall(PetscFree(fields[k]));
135: PetscCall(PetscFree(fields));
136: }
137: *fields = displayfields;
138: *outfields = ndisplayfields;
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: #include <petscdraw.h>
144: #if defined(__GNUC__) && !defined(__clang__)
145: #pragma GCC diagnostic push
146: #pragma GCC diagnostic ignored "-Wclobbered"
147: #endif
148: PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin, PetscViewer v)
149: {
150: DM da;
151: PetscMPIInt rank, size, tag;
152: PetscInt i, n, N, dof, istart, isize, j, nbounds;
153: MPI_Status status;
154: PetscReal min, max, xmin = 0.0, xmax = 0.0, tmp = 0.0, xgtmp = 0.0;
155: const PetscScalar *array, *xg;
156: PetscDraw draw;
157: PetscBool isnull, useports = PETSC_FALSE, showmarkers = PETSC_FALSE;
158: MPI_Comm comm;
159: PetscDrawAxis axis;
160: Vec xcoor;
161: DMBoundaryType bx;
162: const char *tlabel = NULL, *xlabel = NULL;
163: const PetscReal *bounds;
164: PetscInt *displayfields;
165: PetscInt ndisplayfields;
166: PetscBool hold;
167: PetscDrawViewPorts *ports = NULL;
168: PetscViewerFormat format;
170: PetscFunctionBegin;
171: PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
172: PetscCall(PetscDrawIsNull(draw, &isnull));
173: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
174: PetscCall(PetscViewerDrawGetBounds(v, &nbounds, &bounds));
176: PetscCall(VecGetDM(xin, &da));
177: PetscCheck(da, PetscObjectComm((PetscObject)xin), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMDA");
178: PetscCall(PetscObjectGetComm((PetscObject)xin, &comm));
179: PetscCallMPI(MPI_Comm_size(comm, &size));
180: PetscCallMPI(MPI_Comm_rank(comm, &rank));
182: PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_vec_use_markers", &showmarkers, NULL));
184: PetscCall(DMDAGetInfo(da, NULL, &N, NULL, NULL, NULL, NULL, NULL, &dof, NULL, &bx, NULL, NULL, NULL));
185: PetscCall(DMDAGetCorners(da, &istart, NULL, NULL, &isize, NULL, NULL));
186: PetscCall(VecGetArrayRead(xin, &array));
187: PetscCall(VecGetLocalSize(xin, &n));
188: n = n / dof;
190: /* Get coordinates of nodes */
191: PetscCall(DMGetCoordinates(da, &xcoor));
192: if (!xcoor) {
193: PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0));
194: PetscCall(DMGetCoordinates(da, &xcoor));
195: }
196: PetscCall(VecGetArrayRead(xcoor, &xg));
197: PetscCall(DMDAGetCoordinateName(da, 0, &xlabel));
199: /* Determine the min and max coordinate in plot */
200: if (rank == 0) xmin = PetscRealPart(xg[0]);
201: if (rank == size - 1) xmax = PetscRealPart(xg[n - 1]);
202: PetscCallMPI(MPI_Bcast(&xmin, 1, MPIU_REAL, 0, comm));
203: PetscCallMPI(MPI_Bcast(&xmax, 1, MPIU_REAL, size - 1, comm));
205: PetscCall(DMDASelectFields(da, &ndisplayfields, &displayfields));
206: PetscCall(PetscViewerGetFormat(v, &format));
207: PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_ports", &useports, NULL));
208: if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
209: if (useports) {
210: PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
211: PetscCall(PetscViewerDrawGetDrawAxis(v, 0, &axis));
212: PetscCall(PetscDrawCheckResizedWindow(draw));
213: PetscCall(PetscDrawClear(draw));
214: PetscCall(PetscDrawViewPortsCreate(draw, ndisplayfields, &ports));
215: }
217: /* Loop over each field; drawing each in a different window */
218: for (PetscInt k = 0; k < ndisplayfields; k++) {
219: j = displayfields[k];
221: /* determine the min and max value in plot */
222: PetscCall(VecStrideMin(xin, j, NULL, &min));
223: PetscCall(VecStrideMax(xin, j, NULL, &max));
224: if (j < nbounds) {
225: min = PetscMin(min, bounds[2 * j]);
226: max = PetscMax(max, bounds[2 * j + 1]);
227: }
228: if (min == max) {
229: min -= 1.e-5;
230: max += 1.e-5;
231: }
233: if (useports) {
234: PetscCall(PetscDrawViewPortsSet(ports, k));
235: PetscCall(DMDAGetFieldName(da, j, &tlabel));
236: } else {
237: const char *title;
238: PetscCall(PetscViewerDrawGetHold(v, &hold));
239: PetscCall(PetscViewerDrawGetDraw(v, k, &draw));
240: PetscCall(PetscViewerDrawGetDrawAxis(v, k, &axis));
241: PetscCall(DMDAGetFieldName(da, j, &title));
242: if (title) PetscCall(PetscDrawSetTitle(draw, title));
243: PetscCall(PetscDrawCheckResizedWindow(draw));
244: if (!hold) PetscCall(PetscDrawClear(draw));
245: }
246: PetscCall(PetscDrawAxisSetLabels(axis, tlabel, xlabel, NULL));
247: PetscCall(PetscDrawAxisSetLimits(axis, xmin, xmax, min, max));
248: PetscCall(PetscDrawAxisDraw(axis));
250: /* draw local part of vector */
251: PetscCall(PetscObjectGetNewTag((PetscObject)xin, &tag));
252: if (rank < size - 1) { /*send value to right */
253: PetscCallMPI(MPI_Send((void *)&xg[n - 1], 1, MPIU_REAL, rank + 1, tag, comm));
254: PetscCallMPI(MPI_Send((void *)&array[j + (n - 1) * dof], 1, MPIU_REAL, rank + 1, tag, comm));
255: }
256: if (rank) { /* receive value from left */
257: PetscCallMPI(MPI_Recv(&xgtmp, 1, MPIU_REAL, rank - 1, tag, comm, &status));
258: PetscCallMPI(MPI_Recv(&tmp, 1, MPIU_REAL, rank - 1, tag, comm, &status));
259: }
260: PetscDrawCollectiveBegin(draw);
261: if (rank) {
262: PetscCall(PetscDrawLine(draw, xgtmp, tmp, PetscRealPart(xg[0]), PetscRealPart(array[j]), PETSC_DRAW_RED));
263: if (showmarkers) PetscCall(PetscDrawPoint(draw, xgtmp, tmp, PETSC_DRAW_BLACK));
264: }
265: for (i = 1; i < n; i++) {
266: PetscCall(PetscDrawLine(draw, PetscRealPart(xg[i - 1]), PetscRealPart(array[j + dof * (i - 1)]), PetscRealPart(xg[i]), PetscRealPart(array[j + dof * i]), PETSC_DRAW_RED));
267: if (showmarkers) PetscCall(PetscDrawMarker(draw, PetscRealPart(xg[i - 1]), PetscRealPart(array[j + dof * (i - 1)]), PETSC_DRAW_BLACK));
268: }
269: if (rank == size - 1) {
270: if (showmarkers) PetscCall(PetscDrawMarker(draw, PetscRealPart(xg[n - 1]), PetscRealPart(array[j + dof * (n - 1)]), PETSC_DRAW_BLACK));
271: }
272: PetscDrawCollectiveEnd(draw);
273: PetscCall(PetscDrawFlush(draw));
274: PetscCall(PetscDrawPause(draw));
275: if (!useports) PetscCall(PetscDrawSave(draw));
276: }
277: if (useports) {
278: PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
279: PetscCall(PetscDrawSave(draw));
280: }
282: PetscCall(PetscDrawViewPortsDestroy(ports));
283: PetscCall(PetscFree(displayfields));
284: PetscCall(VecRestoreArrayRead(xcoor, &xg));
285: PetscCall(VecRestoreArrayRead(xin, &array));
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
288: #if defined(__GNUC__) && !defined(__clang__)
289: #pragma GCC diagnostic pop
290: #endif