1 /* Routines to visualize DMDAs and fields through GLVis */ 2 3 #include <petsc/private/dmdaimpl.h> 4 #include <petsc/private/glvisviewerimpl.h> 5 6 typedef struct { 7 PetscBool ll; 8 } DMDAGhostedGLVisViewerCtx; 9 10 static PetscErrorCode DMDAGhostedDestroyGLVisViewerCtx_Private(void **vctx) { 11 PetscFunctionBegin; 12 PetscCall(PetscFree(*vctx)); 13 PetscFunctionReturn(0); 14 } 15 16 typedef struct { 17 Vec xlocal; 18 } DMDAFieldGLVisViewerCtx; 19 20 static PetscErrorCode DMDAFieldDestroyGLVisViewerCtx_Private(void *vctx) { 21 DMDAFieldGLVisViewerCtx *ctx = (DMDAFieldGLVisViewerCtx *)vctx; 22 23 PetscFunctionBegin; 24 PetscCall(VecDestroy(&ctx->xlocal)); 25 PetscCall(PetscFree(vctx)); 26 PetscFunctionReturn(0); 27 } 28 29 /* 30 dactx->ll is false -> all but the last proc per dimension claim the ghosted node on the right 31 dactx->ll is true -> all but the first proc per dimension claim the ghosted node on the left 32 */ 33 static PetscErrorCode DMDAGetNumElementsGhosted(DM da, PetscInt *nex, PetscInt *ney, PetscInt *nez) { 34 DMDAGhostedGLVisViewerCtx *dactx; 35 PetscInt sx, sy, sz, ien, jen, ken; 36 37 PetscFunctionBegin; 38 /* Appease -Wmaybe-uninitialized */ 39 if (nex) *nex = -1; 40 if (ney) *ney = -1; 41 if (nez) *nez = -1; 42 PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &ien, &jen, &ken)); 43 PetscCall(DMGetApplicationContext(da, &dactx)); 44 if (dactx->ll) { 45 PetscInt dim; 46 47 PetscCall(DMGetDimension(da, &dim)); 48 if (!sx) ien--; 49 if (!sy && dim > 1) jen--; 50 if (!sz && dim > 2) ken--; 51 } else { 52 PetscInt M, N, P; 53 54 PetscCall(DMDAGetInfo(da, NULL, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 55 if (sx + ien == M) ien--; 56 if (sy + jen == N) jen--; 57 if (sz + ken == P) ken--; 58 } 59 if (nex) *nex = ien; 60 if (ney) *ney = jen; 61 if (nez) *nez = ken; 62 PetscFunctionReturn(0); 63 } 64 65 /* inherits number of vertices from DMDAGetNumElementsGhosted */ 66 static PetscErrorCode DMDAGetNumVerticesGhosted(DM da, PetscInt *nvx, PetscInt *nvy, PetscInt *nvz) { 67 PetscInt ien = 0, jen = 0, ken = 0, dim; 68 PetscInt tote; 69 70 PetscFunctionBegin; 71 PetscCall(DMGetDimension(da, &dim)); 72 PetscCall(DMDAGetNumElementsGhosted(da, &ien, &jen, &ken)); 73 tote = ien * (dim > 1 ? jen : 1) * (dim > 2 ? ken : 1); 74 if (tote) { 75 ien = ien + 1; 76 jen = dim > 1 ? jen + 1 : 1; 77 ken = dim > 2 ? ken + 1 : 1; 78 } 79 if (nvx) *nvx = ien; 80 if (nvy) *nvy = jen; 81 if (nvz) *nvz = ken; 82 PetscFunctionReturn(0); 83 } 84 85 static PetscErrorCode DMDASampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXf[], void *vctx) { 86 DM da; 87 DMDAFieldGLVisViewerCtx *ctx = (DMDAFieldGLVisViewerCtx *)vctx; 88 DMDAGhostedGLVisViewerCtx *dactx; 89 const PetscScalar *array; 90 PetscScalar **arrayf; 91 PetscInt i, f, ii, ien, jen, ken, ie, je, ke, bs, *bss; 92 PetscInt sx, sy, sz, gsx, gsy, gsz, ist, jst, kst, gm, gn, gp; 93 94 PetscFunctionBegin; 95 PetscCall(VecGetDM(ctx->xlocal, &da)); 96 PetscCheck(da, PetscObjectComm(oX), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMDA"); 97 PetscCall(DMGetApplicationContext(da, &dactx)); 98 PetscCall(VecGetBlockSize(ctx->xlocal, &bs)); 99 PetscCall(DMGlobalToLocalBegin(da, (Vec)oX, INSERT_VALUES, ctx->xlocal)); 100 PetscCall(DMGlobalToLocalEnd(da, (Vec)oX, INSERT_VALUES, ctx->xlocal)); 101 PetscCall(DMDAGetNumVerticesGhosted(da, &ien, &jen, &ken)); 102 PetscCall(DMDAGetGhostCorners(da, &gsx, &gsy, &gsz, &gm, &gn, &gp)); 103 PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, NULL, NULL, NULL)); 104 if (dactx->ll) { 105 kst = jst = ist = 0; 106 } else { 107 kst = gsz != sz ? 1 : 0; 108 jst = gsy != sy ? 1 : 0; 109 ist = gsx != sx ? 1 : 0; 110 } 111 PetscCall(PetscMalloc2(nf, &arrayf, nf, &bss)); 112 PetscCall(VecGetArrayRead(ctx->xlocal, &array)); 113 for (f = 0; f < nf; f++) { 114 PetscCall(VecGetBlockSize((Vec)oXf[f], &bss[f])); 115 PetscCall(VecGetArray((Vec)oXf[f], &arrayf[f])); 116 } 117 for (ke = kst, ii = 0; ke < kst + ken; ke++) { 118 for (je = jst; je < jst + jen; je++) { 119 for (ie = ist; ie < ist + ien; ie++) { 120 PetscInt cf, b; 121 i = ke * gm * gn + je * gm + ie; 122 for (f = 0, cf = 0; f < nf; f++) 123 for (b = 0; b < bss[f]; b++) arrayf[f][bss[f] * ii + b] = array[i * bs + cf++]; 124 ii++; 125 } 126 } 127 } 128 for (f = 0; f < nf; f++) PetscCall(VecRestoreArray((Vec)oXf[f], &arrayf[f])); 129 PetscCall(VecRestoreArrayRead(ctx->xlocal, &array)); 130 PetscCall(PetscFree2(arrayf, bss)); 131 PetscFunctionReturn(0); 132 } 133 134 PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_DMDA(PetscObject oda, PetscViewer viewer) { 135 DM da = (DM)oda, daview; 136 137 PetscFunctionBegin; 138 PetscCall(PetscObjectQuery(oda, "GLVisGraphicsDMDAGhosted", (PetscObject *)&daview)); 139 if (!daview) { 140 DMDAGhostedGLVisViewerCtx *dactx; 141 DM dacoord = NULL; 142 Vec xcoor, xcoorl; 143 PetscBool hashocoord = PETSC_FALSE; 144 const PetscInt *lx, *ly, *lz; 145 PetscInt dim, M, N, P, m, n, p, dof, s, i; 146 147 PetscCall(PetscNew(&dactx)); 148 dactx->ll = PETSC_TRUE; /* default to match elements layout obtained by DMDAGetElements */ 149 PetscOptionsBegin(PetscObjectComm(oda), oda->prefix, "GLVis PetscViewer DMDA Options", "PetscViewer"); 150 PetscCall(PetscOptionsBool("-viewer_glvis_dm_da_ll", "Left-looking subdomain view", NULL, dactx->ll, &dactx->ll, NULL)); 151 PetscOptionsEnd(); 152 /* Create a properly ghosted DMDA to visualize the mesh and the fields associated with */ 153 PetscCall(DMDAGetInfo(da, &dim, &M, &N, &P, &m, &n, &p, &dof, &s, NULL, NULL, NULL, NULL)); 154 PetscCall(DMDAGetOwnershipRanges(da, &lx, &ly, &lz)); 155 PetscCall(PetscObjectQuery((PetscObject)da, "_glvis_mesh_coords", (PetscObject *)&xcoor)); 156 if (!xcoor) { 157 PetscCall(DMGetCoordinates(da, &xcoor)); 158 } else { 159 hashocoord = PETSC_TRUE; 160 } 161 PetscCall(PetscInfo(da, "Creating auxiliary DMDA for managing GLVis graphics\n")); 162 switch (dim) { 163 case 1: 164 PetscCall(DMDACreate1d(PetscObjectComm((PetscObject)da), DM_BOUNDARY_NONE, M, dof, 1, lx, &daview)); 165 if (!hashocoord) PetscCall(DMDACreate1d(PetscObjectComm((PetscObject)da), DM_BOUNDARY_NONE, M, 1, 1, lx, &dacoord)); 166 break; 167 case 2: 168 PetscCall(DMDACreate2d(PetscObjectComm((PetscObject)da), DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, M, N, m, n, dof, 1, lx, ly, &daview)); 169 if (!hashocoord) PetscCall(DMDACreate2d(PetscObjectComm((PetscObject)da), DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, M, N, m, n, 2, 1, lx, ly, &dacoord)); 170 break; 171 case 3: 172 PetscCall(DMDACreate3d(PetscObjectComm((PetscObject)da), DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, M, N, P, m, n, p, dof, 1, lx, ly, lz, &daview)); 173 if (!hashocoord) PetscCall(DMDACreate3d(PetscObjectComm((PetscObject)da), DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, M, N, P, m, n, p, 3, 1, lx, ly, lz, &dacoord)); 174 break; 175 default: SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim); 176 } 177 PetscCall(DMSetApplicationContext(daview, dactx)); 178 PetscCall(DMSetApplicationContextDestroy(daview, DMDAGhostedDestroyGLVisViewerCtx_Private)); 179 PetscCall(DMSetUp(daview)); 180 if (!xcoor) { 181 PetscCall(DMDASetUniformCoordinates(daview, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 182 PetscCall(DMGetCoordinates(daview, &xcoor)); 183 } 184 if (dacoord) { 185 PetscCall(DMSetUp(dacoord)); 186 PetscCall(DMCreateLocalVector(dacoord, &xcoorl)); 187 PetscCall(DMGlobalToLocalBegin(dacoord, xcoor, INSERT_VALUES, xcoorl)); 188 PetscCall(DMGlobalToLocalEnd(dacoord, xcoor, INSERT_VALUES, xcoorl)); 189 PetscCall(DMDestroy(&dacoord)); 190 } else { 191 PetscInt ien, jen, ken, nc, nl, cdof, deg; 192 char fecmesh[64]; 193 const char *name; 194 PetscBool flg; 195 196 PetscCall(DMDAGetNumElementsGhosted(daview, &ien, &jen, &ken)); 197 nc = ien * (jen > 0 ? jen : 1) * (ken > 0 ? ken : 1); 198 199 PetscCall(VecGetLocalSize(xcoor, &nl)); 200 PetscCheck(!nc || (nl % nc) == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Incompatible local coordinate size %" PetscInt_FMT " and number of cells %" PetscInt_FMT, nl, nc); 201 PetscCall(VecDuplicate(xcoor, &xcoorl)); 202 PetscCall(VecCopy(xcoor, xcoorl)); 203 PetscCall(VecSetDM(xcoorl, NULL)); 204 PetscCall(PetscObjectGetName((PetscObject)xcoor, &name)); 205 PetscCall(PetscStrbeginswith(name, "FiniteElementCollection:", &flg)); 206 if (!flg) { 207 deg = 0; 208 if (nc && nl) { 209 cdof = nl / (nc * dim); 210 deg = 1; 211 while (1) { 212 PetscInt degd = 1; 213 for (i = 0; i < dim; i++) degd *= (deg + 1); 214 PetscCheck(degd <= cdof, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell dofs %" PetscInt_FMT, cdof); 215 if (degd == cdof) break; 216 deg++; 217 } 218 } 219 PetscCall(PetscSNPrintf(fecmesh, sizeof(fecmesh), "FiniteElementCollection: L2_T1_%" PetscInt_FMT "D_P%" PetscInt_FMT, dim, deg)); 220 PetscCall(PetscObjectSetName((PetscObject)xcoorl, fecmesh)); 221 } else PetscCall(PetscObjectSetName((PetscObject)xcoorl, name)); 222 } 223 224 /* xcoorl is composed with the ghosted DMDA, the ghosted coordinate DMDA (if present) is only available through this vector */ 225 PetscCall(PetscObjectCompose((PetscObject)daview, "GLVisGraphicsCoordsGhosted", (PetscObject)xcoorl)); 226 PetscCall(PetscObjectDereference((PetscObject)xcoorl)); 227 228 /* daview is composed with the original DMDA */ 229 PetscCall(PetscObjectCompose(oda, "GLVisGraphicsDMDAGhosted", (PetscObject)daview)); 230 PetscCall(PetscObjectDereference((PetscObject)daview)); 231 } 232 233 /* customize the viewer if present */ 234 if (viewer) { 235 DMDAFieldGLVisViewerCtx *ctx; 236 DMDAGhostedGLVisViewerCtx *dactx; 237 char fec[64]; 238 Vec xlocal, *Ufield; 239 const char **dafieldname; 240 char **fec_type, **fieldname; 241 PetscInt *nlocal, *bss, *dims; 242 PetscInt dim, M, N, P, dof, s, i, nf; 243 PetscBool bsset; 244 245 PetscCall(DMDAGetInfo(daview, &dim, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 246 PetscCall(DMGetApplicationContext(daview, &dactx)); 247 PetscCall(DMCreateLocalVector(daview, &xlocal)); 248 PetscCall(DMDAGetFieldNames(da, (const char *const **)&dafieldname)); 249 PetscCall(DMDAGetNumVerticesGhosted(daview, &M, &N, &P)); 250 PetscCall(PetscSNPrintf(fec, sizeof(fec), "FiniteElementCollection: H1_%" PetscInt_FMT "D_P1", dim)); 251 PetscCall(PetscMalloc6(dof, &fec_type, dof, &nlocal, dof, &bss, dof, &dims, dof, &fieldname, dof, &Ufield)); 252 for (i = 0; i < dof; i++) bss[i] = 1; 253 nf = dof; 254 255 PetscOptionsBegin(PetscObjectComm(oda), oda->prefix, "GLVis PetscViewer DMDA Field options", "PetscViewer"); 256 PetscCall(PetscOptionsIntArray("-viewer_glvis_dm_da_bs", "Block sizes for subfields; enable vector representation", NULL, bss, &nf, &bsset)); 257 PetscOptionsEnd(); 258 if (bsset) { 259 PetscInt t; 260 for (i = 0, t = 0; i < nf; i++) t += bss[i]; 261 PetscCheck(t == dof, PetscObjectComm(oda), PETSC_ERR_USER, "Sum of block sizes %" PetscInt_FMT " should equal %" PetscInt_FMT, t, dof); 262 } else nf = dof; 263 264 for (i = 0, s = 0; i < nf; i++) { 265 PetscCall(PetscStrallocpy(fec, &fec_type[i])); 266 if (bss[i] == 1) { 267 PetscCall(PetscStrallocpy(dafieldname[s], &fieldname[i])); 268 } else { 269 PetscInt b; 270 size_t tlen = 9; /* "Vector-" + end */ 271 for (b = 0; b < bss[i]; b++) { 272 size_t len; 273 PetscCall(PetscStrlen(dafieldname[s + b], &len)); 274 tlen += len + 1; /* field + "-" */ 275 } 276 PetscCall(PetscMalloc1(tlen, &fieldname[i])); 277 PetscCall(PetscStrcpy(fieldname[i], "Vector-")); 278 for (b = 0; b < bss[i] - 1; b++) { 279 PetscCall(PetscStrcat(fieldname[i], dafieldname[s + b])); 280 PetscCall(PetscStrcat(fieldname[i], "-")); 281 } 282 PetscCall(PetscStrcat(fieldname[i], dafieldname[s + b])); 283 } 284 dims[i] = dim; 285 nlocal[i] = M * N * P * bss[i]; 286 s += bss[i]; 287 } 288 289 /* the viewer context takes ownership of xlocal and destroys it in DMDAFieldDestroyGLVisViewerCtx_Private */ 290 PetscCall(PetscNew(&ctx)); 291 ctx->xlocal = xlocal; 292 293 /* create work vectors */ 294 for (i = 0; i < nf; i++) { 295 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)da), nlocal[i], PETSC_DECIDE, &Ufield[i])); 296 PetscCall(PetscObjectSetName((PetscObject)Ufield[i], fieldname[i])); 297 PetscCall(VecSetBlockSize(Ufield[i], bss[i])); 298 PetscCall(VecSetDM(Ufield[i], da)); 299 } 300 301 PetscCall(PetscViewerGLVisSetFields(viewer, nf, (const char **)fec_type, dims, DMDASampleGLVisFields_Private, (PetscObject *)Ufield, ctx, DMDAFieldDestroyGLVisViewerCtx_Private)); 302 for (i = 0; i < nf; i++) { 303 PetscCall(PetscFree(fec_type[i])); 304 PetscCall(PetscFree(fieldname[i])); 305 PetscCall(VecDestroy(&Ufield[i])); 306 } 307 PetscCall(PetscFree6(fec_type, nlocal, bss, dims, fieldname, Ufield)); 308 } 309 PetscFunctionReturn(0); 310 } 311 312 static PetscErrorCode DMDAView_GLVis_ASCII(DM dm, PetscViewer viewer) { 313 DM da, cda; 314 Vec xcoorl; 315 PetscMPIInt size; 316 const PetscScalar *array; 317 PetscContainer glvis_container; 318 PetscInt dim, sdim, i, vid[8], mid, cid, cdof; 319 PetscInt sx, sy, sz, ie, je, ke, ien, jen, ken, nel; 320 PetscInt gsx, gsy, gsz, gm, gn, gp, kst, jst, ist; 321 PetscBool enabled = PETSC_TRUE, isascii; 322 const char *fmt; 323 324 PetscFunctionBegin; 325 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 326 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 327 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 328 PetscCheck(isascii, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Viewer must be of type VIEWERASCII"); 329 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size)); 330 PetscCheck(size <= 1, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Use single sequential viewers for parallel visualization"); 331 PetscCall(DMGetDimension(dm, &dim)); 332 333 /* get container: determines if a process visualizes is portion of the data or not */ 334 PetscCall(PetscObjectQuery((PetscObject)viewer, "_glvis_info_container", (PetscObject *)&glvis_container)); 335 PetscCheck(glvis_container, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing GLVis container"); 336 { 337 PetscViewerGLVisInfo glvis_info; 338 PetscCall(PetscContainerGetPointer(glvis_container, (void **)&glvis_info)); 339 enabled = glvis_info->enabled; 340 fmt = glvis_info->fmt; 341 } 342 /* this can happen if we are calling DMView outside of VecView_GLVis */ 343 PetscCall(PetscObjectQuery((PetscObject)dm, "GLVisGraphicsDMDAGhosted", (PetscObject *)&da)); 344 if (!da) PetscCall(DMSetUpGLVisViewer_DMDA((PetscObject)dm, NULL)); 345 PetscCall(PetscObjectQuery((PetscObject)dm, "GLVisGraphicsDMDAGhosted", (PetscObject *)&da)); 346 PetscCheck(da, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing GLVis ghosted DMDA"); 347 PetscCall(DMGetCoordinateDim(da, &sdim)); 348 349 PetscCall(PetscViewerASCIIPrintf(viewer, "MFEM mesh v1.0\n")); 350 PetscCall(PetscViewerASCIIPrintf(viewer, "\ndimension\n")); 351 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", dim)); 352 353 if (!enabled) { 354 PetscCall(PetscViewerASCIIPrintf(viewer, "\nelements\n")); 355 PetscCall(PetscViewerASCIIPrintf(viewer, "0\n")); 356 PetscCall(PetscViewerASCIIPrintf(viewer, "\nboundary\n")); 357 PetscCall(PetscViewerASCIIPrintf(viewer, "0\n")); 358 PetscCall(PetscViewerASCIIPrintf(viewer, "\nvertices\n")); 359 PetscCall(PetscViewerASCIIPrintf(viewer, "0\n")); 360 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", sdim)); 361 PetscFunctionReturn(0); 362 } 363 364 PetscCall(DMDAGetNumElementsGhosted(da, &ien, &jen, &ken)); 365 nel = ien; 366 if (dim > 1) nel *= jen; 367 if (dim > 2) nel *= ken; 368 PetscCall(PetscViewerASCIIPrintf(viewer, "\nelements\n")); 369 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", nel)); 370 switch (dim) { 371 case 1: 372 for (ie = 0; ie < ien; ie++) { 373 vid[0] = ie; 374 vid[1] = ie + 1; 375 mid = 1; /* material id */ 376 cid = 1; /* segment */ 377 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", mid, cid, vid[0], vid[1])); 378 } 379 break; 380 case 2: 381 for (je = 0; je < jen; je++) { 382 for (ie = 0; ie < ien; ie++) { 383 vid[0] = je * (ien + 1) + ie; 384 vid[1] = je * (ien + 1) + ie + 1; 385 vid[2] = (je + 1) * (ien + 1) + ie + 1; 386 vid[3] = (je + 1) * (ien + 1) + ie; 387 mid = 1; /* material id */ 388 cid = 3; /* quad */ 389 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", mid, cid, vid[0], vid[1], vid[2], vid[3])); 390 } 391 } 392 break; 393 case 3: 394 for (ke = 0; ke < ken; ke++) { 395 for (je = 0; je < jen; je++) { 396 for (ie = 0; ie < ien; ie++) { 397 vid[0] = ke * (jen + 1) * (ien + 1) + je * (ien + 1) + ie; 398 vid[1] = ke * (jen + 1) * (ien + 1) + je * (ien + 1) + ie + 1; 399 vid[2] = ke * (jen + 1) * (ien + 1) + (je + 1) * (ien + 1) + ie + 1; 400 vid[3] = ke * (jen + 1) * (ien + 1) + (je + 1) * (ien + 1) + ie; 401 vid[4] = (ke + 1) * (jen + 1) * (ien + 1) + je * (ien + 1) + ie; 402 vid[5] = (ke + 1) * (jen + 1) * (ien + 1) + je * (ien + 1) + ie + 1; 403 vid[6] = (ke + 1) * (jen + 1) * (ien + 1) + (je + 1) * (ien + 1) + ie + 1; 404 vid[7] = (ke + 1) * (jen + 1) * (ien + 1) + (je + 1) * (ien + 1) + ie; 405 mid = 1; /* material id */ 406 cid = 5; /* hex */ 407 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", mid, cid, vid[0], vid[1], vid[2], vid[3], vid[4], vid[5], vid[6], vid[7])); 408 } 409 } 410 } 411 break; 412 default: SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim); 413 } 414 PetscCall(PetscViewerASCIIPrintf(viewer, "\nboundary\n")); 415 PetscCall(PetscViewerASCIIPrintf(viewer, "0\n")); 416 417 /* vertex coordinates */ 418 PetscCall(PetscObjectQuery((PetscObject)da, "GLVisGraphicsCoordsGhosted", (PetscObject *)&xcoorl)); 419 PetscCheck(xcoorl, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing GLVis ghosted coords"); 420 PetscCall(DMDAGetNumVerticesGhosted(da, &ien, &jen, &ken)); 421 PetscCall(PetscViewerASCIIPrintf(viewer, "\nvertices\n")); 422 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", ien * jen * ken)); 423 if (nel) { 424 PetscCall(VecGetDM(xcoorl, &cda)); 425 PetscCall(VecGetArrayRead(xcoorl, &array)); 426 if (!cda) { /* HO viz */ 427 const char *fecname; 428 PetscInt nc, nl; 429 430 PetscCall(PetscObjectGetName((PetscObject)xcoorl, &fecname)); 431 PetscCall(PetscViewerASCIIPrintf(viewer, "nodes\n")); 432 PetscCall(PetscViewerASCIIPrintf(viewer, "FiniteElementSpace\n")); 433 PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", fecname)); 434 PetscCall(PetscViewerASCIIPrintf(viewer, "VDim: %" PetscInt_FMT "\n", sdim)); 435 PetscCall(PetscViewerASCIIPrintf(viewer, "Ordering: 1\n\n")); /*Ordering::byVDIM*/ 436 /* L2 coordinates */ 437 PetscCall(DMDAGetNumElementsGhosted(da, &ien, &jen, &ken)); 438 PetscCall(VecGetLocalSize(xcoorl, &nl)); 439 nc = ien * (jen > 0 ? jen : 1) * (ken > 0 ? ken : 1); 440 cdof = nc ? nl / nc : 0; 441 if (!ien) ien++; 442 if (!jen) jen++; 443 if (!ken) ken++; 444 ist = jst = kst = 0; 445 gm = ien; 446 gn = jen; 447 gp = ken; 448 } else { 449 DMDAGhostedGLVisViewerCtx *dactx; 450 451 PetscCall(DMGetApplicationContext(da, &dactx)); 452 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", sdim)); 453 cdof = sdim; 454 PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, NULL, NULL, NULL)); 455 PetscCall(DMDAGetGhostCorners(da, &gsx, &gsy, &gsz, &gm, &gn, &gp)); 456 if (dactx->ll) { 457 kst = jst = ist = 0; 458 } else { 459 kst = gsz != sz ? 1 : 0; 460 jst = gsy != sy ? 1 : 0; 461 ist = gsx != sx ? 1 : 0; 462 } 463 } 464 for (ke = kst; ke < kst + ken; ke++) { 465 for (je = jst; je < jst + jen; je++) { 466 for (ie = ist; ie < ist + ien; ie++) { 467 PetscInt c; 468 469 i = ke * gm * gn + je * gm + ie; 470 for (c = 0; c < cdof / sdim; c++) { 471 PetscInt d; 472 for (d = 0; d < sdim; d++) PetscCall(PetscViewerASCIIPrintf(viewer, fmt, PetscRealPart(array[cdof * i + c * sdim + d]))); 473 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 474 } 475 } 476 } 477 } 478 PetscCall(VecRestoreArrayRead(xcoorl, &array)); 479 } 480 PetscFunctionReturn(0); 481 } 482 483 PetscErrorCode DMView_DA_GLVis(DM dm, PetscViewer viewer) { 484 PetscFunctionBegin; 485 PetscCall(DMView_GLVis(dm, viewer, DMDAView_GLVis_ASCII)); 486 PetscFunctionReturn(0); 487 } 488