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