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 @*/ 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 */ 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 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