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