1 /*
2 Plots vectors obtained with DMDACreate1d()
3 */
4
5 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
6
7 /*@
8 DMDASetUniformCoordinates - Sets a `DMDA` coordinates to be a uniform grid
9
10 Collective
11
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)
20
21 Level: beginner
22
23 .seealso: [](sec_struct), `DM`, `DMDA`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMStagSetUniformCoordinates()`
24 @*/
DMDASetUniformCoordinates(DM da,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)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;
35
36 PetscFunctionBegin;
37 PetscValidHeaderSpecificType(da, DM_CLASSID, 1, DMDA);
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};
94
95 PetscCall(DMSetPeriodicity(da, maxCell, Lstart, L));
96 }
97 PetscFunctionReturn(PETSC_SUCCESS);
98 }
99
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 */
DMDASelectFields(DM da,PetscInt * outfields,PetscInt ** fields)103 PetscErrorCode DMDASelectFields(DM da, PetscInt *outfields, PetscInt **fields)
104 {
105 PetscInt step, ndisplayfields, *displayfields, k, j;
106 PetscBool flg;
107
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 }
141
142 #include <petscdraw.h>
143
VecView_MPI_Draw_DA1d(Vec xin,PetscViewer v)144 PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin, PetscViewer v)
145 {
146 DM da;
147 PetscMPIInt rank, size, tag;
148 PetscInt i, n, N, dof, istart, isize, j, nbounds;
149 MPI_Status status;
150 PetscReal min, max, xmin = 0.0, xmax = 0.0, tmp = 0.0, xgtmp = 0.0;
151 const PetscScalar *array, *xg;
152 PetscDraw draw;
153 PetscBool isnull, useports = PETSC_FALSE, showmarkers = PETSC_FALSE;
154 MPI_Comm comm;
155 PetscDrawAxis axis;
156 Vec xcoor;
157 DMBoundaryType bx;
158 const char *tlabel = NULL, *xlabel = NULL;
159 const PetscReal *bounds;
160 PetscInt *displayfields;
161 PetscInt k, ndisplayfields;
162 PetscBool hold;
163 PetscDrawViewPorts *ports = NULL;
164 PetscViewerFormat format;
165
166 PetscFunctionBegin;
167 PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
168 PetscCall(PetscDrawIsNull(draw, &isnull));
169 if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
170 PetscCall(PetscViewerDrawGetBounds(v, &nbounds, &bounds));
171
172 PetscCall(VecGetDM(xin, &da));
173 PetscCheck(da, PetscObjectComm((PetscObject)xin), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMDA");
174 PetscCall(PetscObjectGetComm((PetscObject)xin, &comm));
175 PetscCallMPI(MPI_Comm_size(comm, &size));
176 PetscCallMPI(MPI_Comm_rank(comm, &rank));
177
178 PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_vec_use_markers", &showmarkers, NULL));
179
180 PetscCall(DMDAGetInfo(da, NULL, &N, NULL, NULL, NULL, NULL, NULL, &dof, NULL, &bx, NULL, NULL, NULL));
181 PetscCall(DMDAGetCorners(da, &istart, NULL, NULL, &isize, NULL, NULL));
182 PetscCall(VecGetArrayRead(xin, &array));
183 PetscCall(VecGetLocalSize(xin, &n));
184 n = n / dof;
185
186 /* Get coordinates of nodes */
187 PetscCall(DMGetCoordinates(da, &xcoor));
188 if (!xcoor) {
189 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0));
190 PetscCall(DMGetCoordinates(da, &xcoor));
191 }
192 PetscCall(VecGetArrayRead(xcoor, &xg));
193 PetscCall(DMDAGetCoordinateName(da, 0, &xlabel));
194
195 /* Determine the min and max coordinate in plot */
196 if (rank == 0) xmin = PetscRealPart(xg[0]);
197 if (rank == size - 1) xmax = PetscRealPart(xg[n - 1]);
198 PetscCallMPI(MPI_Bcast(&xmin, 1, MPIU_REAL, 0, comm));
199 PetscCallMPI(MPI_Bcast(&xmax, 1, MPIU_REAL, size - 1, comm));
200
201 PetscCall(DMDASelectFields(da, &ndisplayfields, &displayfields));
202 PetscCall(PetscViewerGetFormat(v, &format));
203 PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_ports", &useports, NULL));
204 if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
205 if (useports) {
206 PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
207 PetscCall(PetscViewerDrawGetDrawAxis(v, 0, &axis));
208 PetscCall(PetscDrawCheckResizedWindow(draw));
209 PetscCall(PetscDrawClear(draw));
210 PetscCall(PetscDrawViewPortsCreate(draw, ndisplayfields, &ports));
211 }
212
213 /* Loop over each field; drawing each in a different window */
214 for (k = 0; k < ndisplayfields; k++) {
215 j = displayfields[k];
216
217 /* determine the min and max value in plot */
218 PetscCall(VecStrideMin(xin, j, NULL, &min));
219 PetscCall(VecStrideMax(xin, j, NULL, &max));
220 if (j < nbounds) {
221 min = PetscMin(min, bounds[2 * j]);
222 max = PetscMax(max, bounds[2 * j + 1]);
223 }
224 if (min == max) {
225 min -= 1.e-5;
226 max += 1.e-5;
227 }
228
229 if (useports) {
230 PetscCall(PetscDrawViewPortsSet(ports, k));
231 PetscCall(DMDAGetFieldName(da, j, &tlabel));
232 } else {
233 const char *title;
234 PetscCall(PetscViewerDrawGetHold(v, &hold));
235 PetscCall(PetscViewerDrawGetDraw(v, k, &draw));
236 PetscCall(PetscViewerDrawGetDrawAxis(v, k, &axis));
237 PetscCall(DMDAGetFieldName(da, j, &title));
238 if (title) PetscCall(PetscDrawSetTitle(draw, title));
239 PetscCall(PetscDrawCheckResizedWindow(draw));
240 if (!hold) PetscCall(PetscDrawClear(draw));
241 }
242 PetscCall(PetscDrawAxisSetLabels(axis, tlabel, xlabel, NULL));
243 PetscCall(PetscDrawAxisSetLimits(axis, xmin, xmax, min, max));
244 PetscCall(PetscDrawAxisDraw(axis));
245
246 /* draw local part of vector */
247 PetscCall(PetscObjectGetNewTag((PetscObject)xin, &tag));
248 if (rank < size - 1) { /*send value to right */
249 PetscCallMPI(MPI_Send((void *)&xg[n - 1], 1, MPIU_REAL, rank + 1, tag, comm));
250 PetscCallMPI(MPI_Send((void *)&array[j + (n - 1) * dof], 1, MPIU_REAL, rank + 1, tag, comm));
251 }
252 if (rank) { /* receive value from left */
253 PetscCallMPI(MPI_Recv(&xgtmp, 1, MPIU_REAL, rank - 1, tag, comm, &status));
254 PetscCallMPI(MPI_Recv(&tmp, 1, MPIU_REAL, rank - 1, tag, comm, &status));
255 }
256 PetscDrawCollectiveBegin(draw);
257 if (rank) {
258 PetscCall(PetscDrawLine(draw, xgtmp, tmp, PetscRealPart(xg[0]), PetscRealPart(array[j]), PETSC_DRAW_RED));
259 if (showmarkers) PetscCall(PetscDrawPoint(draw, xgtmp, tmp, PETSC_DRAW_BLACK));
260 }
261 for (i = 1; i < n; i++) {
262 PetscCall(PetscDrawLine(draw, PetscRealPart(xg[i - 1]), PetscRealPart(array[j + dof * (i - 1)]), PetscRealPart(xg[i]), PetscRealPart(array[j + dof * i]), PETSC_DRAW_RED));
263 if (showmarkers) PetscCall(PetscDrawMarker(draw, PetscRealPart(xg[i - 1]), PetscRealPart(array[j + dof * (i - 1)]), PETSC_DRAW_BLACK));
264 }
265 if (rank == size - 1) {
266 if (showmarkers) PetscCall(PetscDrawMarker(draw, PetscRealPart(xg[n - 1]), PetscRealPart(array[j + dof * (n - 1)]), PETSC_DRAW_BLACK));
267 }
268 PetscDrawCollectiveEnd(draw);
269 PetscCall(PetscDrawFlush(draw));
270 PetscCall(PetscDrawPause(draw));
271 if (!useports) PetscCall(PetscDrawSave(draw));
272 }
273 if (useports) {
274 PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
275 PetscCall(PetscDrawSave(draw));
276 }
277
278 PetscCall(PetscDrawViewPortsDestroy(ports));
279 PetscCall(PetscFree(displayfields));
280 PetscCall(VecRestoreArrayRead(xcoor, &xg));
281 PetscCall(VecRestoreArrayRead(xin, &array));
282 PetscFunctionReturn(PETSC_SUCCESS);
283 }
284