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