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