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 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; 33 34 PetscFunctionBegin; 35 PetscValidHeaderSpecificType(da, DM_CLASSID, 1, DMDA); 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(PetscLogObjectParent((PetscObject)da, (PetscObject)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 PetscInt step, ndisplayfields, *displayfields, k, j; 96 PetscBool flg; 97 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(0); 130 } 131 132 #include <petscdraw.h> 133 134 PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin, PetscViewer v) { 135 DM da; 136 PetscMPIInt rank, size, tag; 137 PetscInt i, n, N, dof, istart, isize, j, nbounds; 138 MPI_Status status; 139 PetscReal min, max, xmin = 0.0, xmax = 0.0, tmp = 0.0, xgtmp = 0.0; 140 const PetscScalar *array, *xg; 141 PetscDraw draw; 142 PetscBool isnull, useports = PETSC_FALSE, showmarkers = PETSC_FALSE; 143 MPI_Comm comm; 144 PetscDrawAxis axis; 145 Vec xcoor; 146 DMBoundaryType bx; 147 const char *tlabel = NULL, *xlabel = NULL; 148 const PetscReal *bounds; 149 PetscInt *displayfields; 150 PetscInt k, ndisplayfields; 151 PetscBool hold; 152 PetscDrawViewPorts *ports = NULL; 153 PetscViewerFormat format; 154 155 PetscFunctionBegin; 156 PetscCall(PetscViewerDrawGetDraw(v, 0, &draw)); 157 PetscCall(PetscDrawIsNull(draw, &isnull)); 158 if (isnull) PetscFunctionReturn(0); 159 PetscCall(PetscViewerDrawGetBounds(v, &nbounds, &bounds)); 160 161 PetscCall(VecGetDM(xin, &da)); 162 PetscCheck(da, PetscObjectComm((PetscObject)xin), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMDA"); 163 PetscCall(PetscObjectGetComm((PetscObject)xin, &comm)); 164 PetscCallMPI(MPI_Comm_size(comm, &size)); 165 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 166 167 PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_vec_use_markers", &showmarkers, NULL)); 168 169 PetscCall(DMDAGetInfo(da, NULL, &N, NULL, NULL, NULL, NULL, NULL, &dof, NULL, &bx, NULL, NULL, NULL)); 170 PetscCall(DMDAGetCorners(da, &istart, NULL, NULL, &isize, NULL, NULL)); 171 PetscCall(VecGetArrayRead(xin, &array)); 172 PetscCall(VecGetLocalSize(xin, &n)); 173 n = n / dof; 174 175 /* Get coordinates of nodes */ 176 PetscCall(DMGetCoordinates(da, &xcoor)); 177 if (!xcoor) { 178 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0)); 179 PetscCall(DMGetCoordinates(da, &xcoor)); 180 } 181 PetscCall(VecGetArrayRead(xcoor, &xg)); 182 PetscCall(DMDAGetCoordinateName(da, 0, &xlabel)); 183 184 /* Determine the min and max coordinate in plot */ 185 if (rank == 0) xmin = PetscRealPart(xg[0]); 186 if (rank == size - 1) xmax = PetscRealPart(xg[n - 1]); 187 PetscCallMPI(MPI_Bcast(&xmin, 1, MPIU_REAL, 0, comm)); 188 PetscCallMPI(MPI_Bcast(&xmax, 1, MPIU_REAL, size - 1, comm)); 189 190 PetscCall(DMDASelectFields(da, &ndisplayfields, &displayfields)); 191 PetscCall(PetscViewerGetFormat(v, &format)); 192 PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_ports", &useports, NULL)); 193 if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE; 194 if (useports) { 195 PetscCall(PetscViewerDrawGetDraw(v, 0, &draw)); 196 PetscCall(PetscViewerDrawGetDrawAxis(v, 0, &axis)); 197 PetscCall(PetscDrawCheckResizedWindow(draw)); 198 PetscCall(PetscDrawClear(draw)); 199 PetscCall(PetscDrawViewPortsCreate(draw, ndisplayfields, &ports)); 200 } 201 202 /* Loop over each field; drawing each in a different window */ 203 for (k = 0; k < ndisplayfields; k++) { 204 j = displayfields[k]; 205 206 /* determine the min and max value in plot */ 207 PetscCall(VecStrideMin(xin, j, NULL, &min)); 208 PetscCall(VecStrideMax(xin, j, NULL, &max)); 209 if (j < nbounds) { 210 min = PetscMin(min, bounds[2 * j]); 211 max = PetscMax(max, bounds[2 * j + 1]); 212 } 213 if (min == max) { 214 min -= 1.e-5; 215 max += 1.e-5; 216 } 217 218 if (useports) { 219 PetscCall(PetscDrawViewPortsSet(ports, k)); 220 PetscCall(DMDAGetFieldName(da, j, &tlabel)); 221 } else { 222 const char *title; 223 PetscCall(PetscViewerDrawGetHold(v, &hold)); 224 PetscCall(PetscViewerDrawGetDraw(v, k, &draw)); 225 PetscCall(PetscViewerDrawGetDrawAxis(v, k, &axis)); 226 PetscCall(DMDAGetFieldName(da, j, &title)); 227 if (title) PetscCall(PetscDrawSetTitle(draw, title)); 228 PetscCall(PetscDrawCheckResizedWindow(draw)); 229 if (!hold) PetscCall(PetscDrawClear(draw)); 230 } 231 PetscCall(PetscDrawAxisSetLabels(axis, tlabel, xlabel, NULL)); 232 PetscCall(PetscDrawAxisSetLimits(axis, xmin, xmax, min, max)); 233 PetscCall(PetscDrawAxisDraw(axis)); 234 235 /* draw local part of vector */ 236 PetscCall(PetscObjectGetNewTag((PetscObject)xin, &tag)); 237 if (rank < size - 1) { /*send value to right */ 238 PetscCallMPI(MPI_Send((void *)&xg[n - 1], 1, MPIU_REAL, rank + 1, tag, comm)); 239 PetscCallMPI(MPI_Send((void *)&array[j + (n - 1) * dof], 1, MPIU_REAL, rank + 1, tag, comm)); 240 } 241 if (rank) { /* receive value from left */ 242 PetscCallMPI(MPI_Recv(&xgtmp, 1, MPIU_REAL, rank - 1, tag, comm, &status)); 243 PetscCallMPI(MPI_Recv(&tmp, 1, MPIU_REAL, rank - 1, tag, comm, &status)); 244 } 245 PetscDrawCollectiveBegin(draw); 246 if (rank) { 247 PetscCall(PetscDrawLine(draw, xgtmp, tmp, PetscRealPart(xg[0]), PetscRealPart(array[j]), PETSC_DRAW_RED)); 248 if (showmarkers) PetscCall(PetscDrawPoint(draw, xgtmp, tmp, PETSC_DRAW_BLACK)); 249 } 250 for (i = 1; i < n; i++) { 251 PetscCall(PetscDrawLine(draw, PetscRealPart(xg[i - 1]), PetscRealPart(array[j + dof * (i - 1)]), PetscRealPart(xg[i]), PetscRealPart(array[j + dof * i]), PETSC_DRAW_RED)); 252 if (showmarkers) PetscCall(PetscDrawMarker(draw, PetscRealPart(xg[i - 1]), PetscRealPart(array[j + dof * (i - 1)]), PETSC_DRAW_BLACK)); 253 } 254 if (rank == size - 1) { 255 if (showmarkers) PetscCall(PetscDrawMarker(draw, PetscRealPart(xg[n - 1]), PetscRealPart(array[j + dof * (n - 1)]), PETSC_DRAW_BLACK)); 256 } 257 PetscDrawCollectiveEnd(draw); 258 PetscCall(PetscDrawFlush(draw)); 259 PetscCall(PetscDrawPause(draw)); 260 if (!useports) PetscCall(PetscDrawSave(draw)); 261 } 262 if (useports) { 263 PetscCall(PetscViewerDrawGetDraw(v, 0, &draw)); 264 PetscCall(PetscDrawSave(draw)); 265 } 266 267 PetscCall(PetscDrawViewPortsDestroy(ports)); 268 PetscCall(PetscFree(displayfields)); 269 PetscCall(VecRestoreArrayRead(xcoor, &xg)); 270 PetscCall(VecRestoreArrayRead(xin, &array)); 271 PetscFunctionReturn(0); 272 } 273