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