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