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