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