xref: /petsc/src/dm/impls/da/grglvis.c (revision 12801b3912e16a0709a58df21c3d316f86673188) 
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,*Ufield;
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 = PetscMalloc6(dof,&fec_type,dof,&nlocal,dof,&bss,dof,&dims,dof,&fieldname,dof,&Ufield);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     /* create work vectors */
302     for (i=0;i<nf;i++) {
303       ierr = VecCreateMPI(PetscObjectComm((PetscObject)da),nlocal[i],PETSC_DECIDE,&Ufield[i]);CHKERRQ(ierr);
304       ierr = PetscObjectSetName((PetscObject)Ufield[i],fieldname[i]);CHKERRQ(ierr);
305       ierr = VecSetBlockSize(Ufield[i],bss[i]);CHKERRQ(ierr);
306       ierr = VecSetDM(Ufield[i],da);CHKERRQ(ierr);
307     }
308 
309     ierr = PetscViewerGLVisSetFields(viewer,nf,(const char**)fec_type,dims,DMDASampleGLVisFields_Private,(PetscObject*)Ufield,ctx,DMDAFieldDestroyGLVisViewerCtx_Private);CHKERRQ(ierr);
310     for (i=0;i<nf;i++) {
311       ierr = PetscFree(fec_type[i]);CHKERRQ(ierr);
312       ierr = PetscFree(fieldname[i]);CHKERRQ(ierr);
313       ierr = VecDestroy(&Ufield[i]);CHKERRQ(ierr);
314     }
315     ierr = PetscFree6(fec_type,nlocal,bss,dims,fieldname,Ufield);CHKERRQ(ierr);
316   }
317   PetscFunctionReturn(0);
318 }
319 
320 static PetscErrorCode DMDAView_GLVis_ASCII(DM dm, PetscViewer viewer)
321 {
322   DM                da,cda;
323   Vec               xcoorl;
324   PetscMPIInt       commsize;
325   const PetscScalar *array;
326   PetscContainer    glvis_container;
327   PetscInt          dim,sdim,i,vid[8],mid,cid,cdof;
328   PetscInt          sx,sy,sz,ie,je,ke,ien,jen,ken;
329   PetscInt          gsx,gsy,gsz,gm,gn,gp,kst,jst,ist;
330   PetscBool         enabled = PETSC_TRUE, isascii;
331   PetscErrorCode    ierr;
332   const char        *fmt;
333 
334   PetscFunctionBegin;
335   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
336   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
337   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
338   if (!isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERASCII");
339   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)viewer),&commsize);CHKERRQ(ierr);
340   if (commsize > 1) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Use single sequential viewers for parallel visualization");
341   ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
342 
343   /* get container: determines if a process visualizes is portion of the data or not */
344   ierr = PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);CHKERRQ(ierr);
345   if (!glvis_container) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing GLVis container");
346   {
347     PetscViewerGLVisInfo glvis_info;
348     ierr = PetscContainerGetPointer(glvis_container,(void**)&glvis_info);CHKERRQ(ierr);
349     enabled = glvis_info->enabled;
350     fmt = glvis_info->fmt;
351   }
352   /* this can happen if we are calling DMView outside of VecView_GLVis */
353   ierr = PetscObjectQuery((PetscObject)dm,"GLVisGraphicsDMDAGhosted",(PetscObject*)&da);CHKERRQ(ierr);
354   if (!da) {ierr = DMSetUpGLVisViewer_DMDA((PetscObject)dm,NULL);CHKERRQ(ierr);}
355   ierr = PetscObjectQuery((PetscObject)dm,"GLVisGraphicsDMDAGhosted",(PetscObject*)&da);CHKERRQ(ierr);
356   if (!da) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing GLVis ghosted DMDA");
357   ierr = DMGetCoordinateDim(da,&sdim);CHKERRQ(ierr);
358 
359   ierr = PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");CHKERRQ(ierr);
360   ierr = PetscViewerASCIIPrintf(viewer,"\ndimension\n");CHKERRQ(ierr);
361   ierr = PetscViewerASCIIPrintf(viewer,"%D\n",dim);CHKERRQ(ierr);
362 
363   if (!enabled) {
364     ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr);
365     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr);
366     ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr);
367     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr);
368     ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr);
369     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr);
370     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr);
371     PetscFunctionReturn(0);
372   }
373 
374   ierr = DMDAGetNumElementsGhosted(da,&ien,&jen,&ken);CHKERRQ(ierr);
375   i    = ien;
376   if (dim > 1) i *= jen;
377   if (dim > 2) i *= ken;
378   ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr);
379   ierr = PetscViewerASCIIPrintf(viewer,"%D\n",i);CHKERRQ(ierr);
380   switch (dim) {
381   case 1:
382     for (ie = 0; ie < ien; ie++) {
383       vid[0] = ie;
384       vid[1] = ie+1;
385       mid    = 1; /* material id */
386       cid    = 1; /* segment */
387       ierr   = PetscViewerASCIIPrintf(viewer,"%D %D %D %D\n",mid,cid,vid[0],vid[1]);CHKERRQ(ierr);
388     }
389     break;
390   case 2:
391     for (je = 0; je < jen; je++) {
392       for (ie = 0; ie < ien; ie++) {
393         vid[0] =     je*(ien+1) + ie;
394         vid[1] =     je*(ien+1) + ie+1;
395         vid[2] = (je+1)*(ien+1) + ie+1;
396         vid[3] = (je+1)*(ien+1) + ie;
397         mid    = 1; /* material id */
398         cid    = 3; /* quad */
399         ierr   = PetscViewerASCIIPrintf(viewer,"%D %D %D %D %D %D\n",mid,cid,vid[0],vid[1],vid[2],vid[3]);CHKERRQ(ierr);
400       }
401     }
402     break;
403   case 3:
404     for (ke = 0; ke < ken; ke++) {
405       for (je = 0; je < jen; je++) {
406         for (ie = 0; ie < ien; ie++) {
407           vid[0] =     ke*(jen+1)*(ien+1) +     je*(ien+1) + ie;
408           vid[1] =     ke*(jen+1)*(ien+1) +     je*(ien+1) + ie+1;
409           vid[2] =     ke*(jen+1)*(ien+1) + (je+1)*(ien+1) + ie+1;
410           vid[3] =     ke*(jen+1)*(ien+1) + (je+1)*(ien+1) + ie;
411           vid[4] = (ke+1)*(jen+1)*(ien+1) +     je*(ien+1) + ie;
412           vid[5] = (ke+1)*(jen+1)*(ien+1) +     je*(ien+1) + ie+1;
413           vid[6] = (ke+1)*(jen+1)*(ien+1) + (je+1)*(ien+1) + ie+1;
414           vid[7] = (ke+1)*(jen+1)*(ien+1) + (je+1)*(ien+1) + ie;
415           mid    = 1; /* material id */
416           cid    = 5; /* hex */
417           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);
418         }
419       }
420     }
421     break;
422   default:
423     SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Unsupported dimension %D",dim);
424     break;
425   }
426   ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr);
427   ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr);
428 
429   /* vertex coordinates */
430   ierr = PetscObjectQuery((PetscObject)da,"GLVisGraphicsCoordsGhosted",(PetscObject*)&xcoorl);CHKERRQ(ierr);
431   if (!xcoorl) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing GLVis ghosted coords");
432   ierr = DMDAGetNumVerticesGhosted(da,&ien,&jen,&ken);CHKERRQ(ierr);
433   ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr);
434   ierr = PetscViewerASCIIPrintf(viewer,"%D\n",ien*jen*ken);CHKERRQ(ierr);
435   ierr = VecGetDM(xcoorl,&cda);CHKERRQ(ierr);
436   ierr = VecGetArrayRead(xcoorl,&array);CHKERRQ(ierr);
437   if (!cda) { /* HO viz */
438     const char *fecname;
439     PetscInt   nc,nl;
440 
441     ierr = PetscObjectGetName((PetscObject)xcoorl,&fecname);CHKERRQ(ierr);
442     ierr = PetscViewerASCIIPrintf(viewer,"nodes\n");CHKERRQ(ierr);
443     ierr = PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");CHKERRQ(ierr);
444     ierr = PetscViewerASCIIPrintf(viewer,"FiniteElementCollection: %s\n",fecname);CHKERRQ(ierr);
445     ierr = PetscViewerASCIIPrintf(viewer,"VDim: %D\n",sdim);CHKERRQ(ierr);
446     ierr = PetscViewerASCIIPrintf(viewer,"Ordering: 1\n\n");CHKERRQ(ierr); /*Ordering::byVDIM*/
447     /* L2 coordinates */
448     ierr = DMDAGetNumElementsGhosted(da,&ien,&jen,&ken);CHKERRQ(ierr);
449     ierr = VecGetLocalSize(xcoorl,&nl);CHKERRQ(ierr);
450     nc   = ien*(jen>0 ? jen : 1)*(ken>0 ? ken : 1);
451     cdof = nl/nc;
452     if (!ien) ien++;
453     if (!jen) jen++;
454     if (!ken) ken++;
455     ist = jst = kst = 0;
456     gm = ien;
457     gn = jen;
458     gp = ken;
459   } else {
460     DMDAGhostedGLVisViewerCtx *dactx;
461 
462     ierr = DMGetApplicationContext(da,(void**)&dactx);CHKERRQ(ierr);
463     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr);
464     cdof = sdim;
465     ierr = DMDAGetCorners(da,&sx,&sy,&sz,NULL,NULL,NULL);CHKERRQ(ierr);
466     ierr = DMDAGetGhostCorners(da,&gsx,&gsy,&gsz,&gm,&gn,&gp);CHKERRQ(ierr);
467     if (dactx->ll) {
468       kst = jst = ist = 0;
469     } else {
470       kst  = gsz != sz ? 1 : 0;
471       jst  = gsy != sy ? 1 : 0;
472       ist  = gsx != sx ? 1 : 0;
473     }
474   }
475   for (ke = kst; ke < kst + ken; ke++) {
476     for (je = jst; je < jst + jen; je++) {
477       for (ie = ist; ie < ist + ien; ie++) {
478         PetscInt c;
479 
480         i = ke * gm * gn + je * gm + ie;
481         for (c=0;c<cdof/sdim;c++) {
482           PetscInt d;
483           for (d=0;d<sdim;d++) {
484             ierr = PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[cdof*i+c*sdim+d]));CHKERRQ(ierr);
485           }
486           ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
487         }
488       }
489     }
490   }
491   ierr = VecRestoreArrayRead(xcoorl,&array);CHKERRQ(ierr);
492   PetscFunctionReturn(0);
493 }
494 
495 /* 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 */
496 PETSC_INTERN PetscErrorCode DMView_DA_GLVis(DM dm, PetscViewer viewer)
497 {
498   PetscErrorCode ierr;
499   PetscBool      isglvis,isascii;
500 
501   PetscFunctionBegin;
502   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
503   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,1);
504   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);CHKERRQ(ierr);
505   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
506   if (!isglvis && !isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERGLVIS or VIEWERASCII");
507   if (isglvis) {
508     PetscViewer          view;
509     PetscViewerGLVisType type;
510 
511     ierr = PetscViewerGLVisGetType_Private(viewer,&type);CHKERRQ(ierr);
512     ierr = PetscViewerGLVisGetDMWindow_Private(viewer,&view);CHKERRQ(ierr);
513     if (view) { /* in the socket case, it may happen that the connection failed */
514       if (type == PETSC_VIEWER_GLVIS_SOCKET) {
515         PetscMPIInt size,rank;
516         ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);CHKERRQ(ierr);
517         ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
518         ierr = PetscViewerASCIIPrintf(view,"parallel %D %D\nmesh\n",size,rank);CHKERRQ(ierr);
519       }
520       ierr = DMDAView_GLVis_ASCII(dm,view);CHKERRQ(ierr);
521       ierr = PetscViewerFlush(view);CHKERRQ(ierr);
522       if (type == PETSC_VIEWER_GLVIS_SOCKET) {
523         PetscInt    dim;
524         const char* name;
525 
526         ierr = PetscObjectGetName((PetscObject)dm,&name);CHKERRQ(ierr);
527         ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
528         ierr = PetscViewerGLVisInitWindow_Private(view,PETSC_TRUE,dim,name);CHKERRQ(ierr);
529         ierr = PetscBarrier((PetscObject)dm);CHKERRQ(ierr);
530       }
531     }
532   } else {
533     ierr = DMDAView_GLVis_ASCII(dm,viewer);CHKERRQ(ierr);
534   }
535   PetscFunctionReturn(0);
536 }
537