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