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