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: }