xref: /petsc/src/dm/impls/plex/plexglvis.c (revision 2da392cc7c10228af19ad9843ce5155178acb644) 
1 #include <petsc/private/glvisviewerimpl.h>
2 #include <petsc/private/petscimpl.h>
3 #include <petsc/private/dmpleximpl.h>
4 #include <petscbt.h>
5 #include <petscdmplex.h>
6 #include <petscsf.h>
7 #include <petscds.h>
8 
9 typedef struct {
10   PetscInt   nf;
11   VecScatter *scctx;
12 } GLVisViewerCtx;
13 
14 static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
15 {
16   GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
17   PetscInt       i;
18   PetscErrorCode ierr;
19 
20   PetscFunctionBegin;
21   for (i=0;i<ctx->nf;i++) {
22     ierr = VecScatterDestroy(&ctx->scctx[i]);CHKERRQ(ierr);
23   }
24   ierr = PetscFree(ctx->scctx);CHKERRQ(ierr);
25   ierr = PetscFree(vctx);CHKERRQ(ierr);
26   PetscFunctionReturn(0);
27 }
28 
29 static PetscErrorCode DMPlexSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
30 {
31   GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
32   PetscInt       f;
33   PetscErrorCode ierr;
34 
35   PetscFunctionBegin;
36   for (f=0;f<nf;f++) {
37     ierr = VecScatterBegin(ctx->scctx[f],(Vec)oX,(Vec)oXfield[f],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
38     ierr = VecScatterEnd(ctx->scctx[f],(Vec)oX,(Vec)oXfield[f],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
39   }
40   PetscFunctionReturn(0);
41 }
42 
43 /* for FEM, it works for H1 fields only and extracts dofs at cell vertices, discarding any other dof */
44 PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject odm, PetscViewer viewer)
45 {
46   DM             dm = (DM)odm;
47   Vec            xlocal,xfield,*Ufield;
48   PetscDS        ds;
49   IS             globalNum,isfield;
50   PetscBT        vown;
51   char           **fieldname = NULL,**fec_type = NULL;
52   const PetscInt *gNum;
53   PetscInt       *nlocal,*bs,*idxs,*dims;
54   PetscInt       f,maxfields,nfields,c,totc,totdofs,Nv,cum,i;
55   PetscInt       dim,cStart,cEnd,vStart,vEnd;
56   GLVisViewerCtx *ctx;
57   PetscSection   s;
58   PetscErrorCode ierr;
59 
60   PetscFunctionBegin;
61   ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
62   ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr);
63   ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr);
64   ierr = PetscObjectQuery((PetscObject)dm,"_glvis_plex_gnum",(PetscObject*)&globalNum);CHKERRQ(ierr);
65   if (!globalNum) {
66     ierr = DMPlexCreateCellNumbering_Internal(dm,PETSC_TRUE,&globalNum);CHKERRQ(ierr);
67     ierr = PetscObjectCompose((PetscObject)dm,"_glvis_plex_gnum",(PetscObject)globalNum);CHKERRQ(ierr);
68     ierr = PetscObjectDereference((PetscObject)globalNum);CHKERRQ(ierr);
69   }
70   ierr = ISGetIndices(globalNum,&gNum);CHKERRQ(ierr);
71   ierr = PetscBTCreate(vEnd-vStart,&vown);CHKERRQ(ierr);
72   for (c = cStart, totc = 0; c < cEnd; c++) {
73     if (gNum[c-cStart] >= 0) {
74       PetscInt i,numPoints,*points = NULL;
75 
76       totc++;
77       ierr = DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr);
78       for (i=0;i<numPoints*2;i+= 2) {
79         if ((points[i] >= vStart) && (points[i] < vEnd)) {
80           ierr = PetscBTSet(vown,points[i]-vStart);CHKERRQ(ierr);
81         }
82       }
83       ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr);
84     }
85   }
86   for (f=0,Nv=0;f<vEnd-vStart;f++) if (PetscLikely(PetscBTLookup(vown,f))) Nv++;
87 
88   ierr = DMCreateLocalVector(dm,&xlocal);CHKERRQ(ierr);
89   ierr = VecGetLocalSize(xlocal,&totdofs);CHKERRQ(ierr);
90   ierr = DMGetLocalSection(dm,&s);CHKERRQ(ierr);
91   ierr = PetscSectionGetNumFields(s,&nfields);CHKERRQ(ierr);
92   for (f=0,maxfields=0;f<nfields;f++) {
93     PetscInt bs;
94 
95     ierr = PetscSectionGetFieldComponents(s,f,&bs);CHKERRQ(ierr);
96     maxfields += bs;
97   }
98   ierr = PetscCalloc7(maxfields,&fieldname,maxfields,&nlocal,maxfields,&bs,maxfields,&dims,maxfields,&fec_type,totdofs,&idxs,maxfields,&Ufield);CHKERRQ(ierr);
99   ierr = PetscNew(&ctx);CHKERRQ(ierr);
100   ierr = PetscCalloc1(maxfields,&ctx->scctx);CHKERRQ(ierr);
101   ierr = DMGetDS(dm,&ds);CHKERRQ(ierr);
102   if (ds) {
103     for (f=0;f<nfields;f++) {
104       const char* fname;
105       char        name[256];
106       PetscObject disc;
107       size_t      len;
108 
109       ierr = PetscSectionGetFieldName(s,f,&fname);CHKERRQ(ierr);
110       ierr = PetscStrlen(fname,&len);CHKERRQ(ierr);
111       if (len) {
112         ierr = PetscStrcpy(name,fname);CHKERRQ(ierr);
113       } else {
114         ierr = PetscSNPrintf(name,256,"Field%D",f);CHKERRQ(ierr);
115       }
116       ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr);
117       if (disc) {
118         PetscClassId id;
119         PetscInt     Nc;
120         char         fec[64];
121 
122         ierr = PetscObjectGetClassId(disc, &id);CHKERRQ(ierr);
123         if (id == PETSCFE_CLASSID) {
124           PetscFE            fem = (PetscFE)disc;
125           PetscDualSpace     sp;
126           PetscDualSpaceType spname;
127           PetscInt           order;
128           PetscBool          islag,continuous,H1 = PETSC_TRUE;
129 
130           ierr = PetscFEGetNumComponents(fem,&Nc);CHKERRQ(ierr);
131           ierr = PetscFEGetDualSpace(fem,&sp);CHKERRQ(ierr);
132           ierr = PetscDualSpaceGetType(sp,&spname);CHKERRQ(ierr);
133           ierr = PetscStrcmp(spname,PETSCDUALSPACELAGRANGE,&islag);CHKERRQ(ierr);
134           if (!islag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported dual space");
135           ierr = PetscDualSpaceLagrangeGetContinuity(sp,&continuous);CHKERRQ(ierr);
136           ierr = PetscDualSpaceGetOrder(sp,&order);CHKERRQ(ierr);
137           if (continuous && order > 0) { /* no support for high-order viz, still have to figure out the numbering */
138             ierr = PetscSNPrintf(fec,64,"FiniteElementCollection: H1_%DD_P1",dim);CHKERRQ(ierr);
139           } else {
140             if (!continuous && order) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Discontinuous space visualization currently unsupported for order %D",order);
141             H1   = PETSC_FALSE;
142             ierr = PetscSNPrintf(fec,64,"FiniteElementCollection: L2_%DD_P%D",dim,order);CHKERRQ(ierr);
143           }
144           ierr = PetscStrallocpy(name,&fieldname[ctx->nf]);CHKERRQ(ierr);
145           bs[ctx->nf]   = Nc;
146           dims[ctx->nf] = dim;
147           if (H1) {
148             nlocal[ctx->nf] = Nc * Nv;
149             ierr = PetscStrallocpy(fec,&fec_type[ctx->nf]);CHKERRQ(ierr);
150             ierr = VecCreateSeq(PETSC_COMM_SELF,Nv*Nc,&xfield);CHKERRQ(ierr);
151             for (i=0,cum=0;i<vEnd-vStart;i++) {
152               PetscInt j,off;
153 
154               if (PetscUnlikely(!PetscBTLookup(vown,i))) continue;
155               ierr = PetscSectionGetFieldOffset(s,i+vStart,f,&off);CHKERRQ(ierr);
156               for (j=0;j<Nc;j++) idxs[cum++] = off + j;
157             }
158             ierr = ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),Nv*Nc,idxs,PETSC_USE_POINTER,&isfield);CHKERRQ(ierr);
159           } else {
160             nlocal[ctx->nf] = Nc * totc;
161             ierr = PetscStrallocpy(fec,&fec_type[ctx->nf]);CHKERRQ(ierr);
162             ierr = VecCreateSeq(PETSC_COMM_SELF,Nc*totc,&xfield);CHKERRQ(ierr);
163             for (i=0,cum=0;i<cEnd-cStart;i++) {
164               PetscInt j,off;
165 
166               if (PetscUnlikely(gNum[i] < 0)) continue;
167               ierr = PetscSectionGetFieldOffset(s,i+cStart,f,&off);CHKERRQ(ierr);
168               for (j=0;j<Nc;j++) idxs[cum++] = off + j;
169             }
170             ierr = ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),totc*Nc,idxs,PETSC_USE_POINTER,&isfield);CHKERRQ(ierr);
171           }
172           ierr = VecScatterCreate(xlocal,isfield,xfield,NULL,&ctx->scctx[ctx->nf]);CHKERRQ(ierr);
173           ierr = VecDestroy(&xfield);CHKERRQ(ierr);
174           ierr = ISDestroy(&isfield);CHKERRQ(ierr);
175           ctx->nf++;
176         } else if (id == PETSCFV_CLASSID) {
177           PetscInt c;
178 
179           ierr = PetscFVGetNumComponents((PetscFV)disc,&Nc);CHKERRQ(ierr);
180           ierr = PetscSNPrintf(fec,64,"FiniteElementCollection: L2_%DD_P0",dim);CHKERRQ(ierr);
181           for (c = 0; c < Nc; c++) {
182             char comp[256];
183             ierr = PetscSNPrintf(comp,256,"%s-Comp%D",name,c);CHKERRQ(ierr);
184             ierr = PetscStrallocpy(comp,&fieldname[ctx->nf]);CHKERRQ(ierr);
185             bs[ctx->nf] = 1; /* Does PetscFV support components with different block size? */
186             nlocal[ctx->nf] = totc;
187             dims[ctx->nf] = dim;
188             ierr = PetscStrallocpy(fec,&fec_type[ctx->nf]);CHKERRQ(ierr);
189             ierr = VecCreateSeq(PETSC_COMM_SELF,totc,&xfield);CHKERRQ(ierr);
190             for (i=0,cum=0;i<cEnd-cStart;i++) {
191               PetscInt off;
192 
193               if (PetscUnlikely(gNum[i])<0) continue;
194               ierr = PetscSectionGetFieldOffset(s,i+cStart,f,&off);CHKERRQ(ierr);
195               idxs[cum++] = off + c;
196             }
197             ierr = ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),totc,idxs,PETSC_USE_POINTER,&isfield);CHKERRQ(ierr);
198             ierr = VecScatterCreate(xlocal,isfield,xfield,NULL,&ctx->scctx[ctx->nf]);CHKERRQ(ierr);
199             ierr = VecDestroy(&xfield);CHKERRQ(ierr);
200             ierr = ISDestroy(&isfield);CHKERRQ(ierr);
201             ctx->nf++;
202           }
203         } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %D",f);
204       } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing discretization for field %D",f);
205     }
206   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Needs a DS attached to the DM");
207   ierr = PetscBTDestroy(&vown);CHKERRQ(ierr);
208   ierr = VecDestroy(&xlocal);CHKERRQ(ierr);
209   ierr = ISRestoreIndices(globalNum,&gNum);CHKERRQ(ierr);
210 
211   /* create work vectors */
212   for (f=0;f<ctx->nf;f++) {
213     ierr = VecCreateMPI(PetscObjectComm((PetscObject)dm),nlocal[f],PETSC_DECIDE,&Ufield[f]);CHKERRQ(ierr);
214     ierr = PetscObjectSetName((PetscObject)Ufield[f],fieldname[f]);CHKERRQ(ierr);
215     ierr = VecSetBlockSize(Ufield[f],bs[f]);CHKERRQ(ierr);
216     ierr = VecSetDM(Ufield[f],dm);CHKERRQ(ierr);
217   }
218 
219   /* customize the viewer */
220   ierr = PetscViewerGLVisSetFields(viewer,ctx->nf,(const char**)fec_type,dims,DMPlexSampleGLVisFields_Private,(PetscObject*)Ufield,ctx,DestroyGLVisViewerCtx_Private);CHKERRQ(ierr);
221   for (f=0;f<ctx->nf;f++) {
222     ierr = PetscFree(fieldname[f]);CHKERRQ(ierr);
223     ierr = PetscFree(fec_type[f]);CHKERRQ(ierr);
224     ierr = VecDestroy(&Ufield[f]);CHKERRQ(ierr);
225   }
226   ierr = PetscFree7(fieldname,nlocal,bs,dims,fec_type,idxs,Ufield);CHKERRQ(ierr);
227   PetscFunctionReturn(0);
228 }
229 
230 typedef enum {MFEM_POINT=0,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE,MFEM_TETRAHEDRON,MFEM_CUBE,MFEM_PRISM,MFEM_UNDEF} MFEM_cid;
231 
232 MFEM_cid mfem_table_cid[4][7]       = { {MFEM_POINT,MFEM_UNDEF,MFEM_UNDEF  ,MFEM_UNDEF   ,MFEM_UNDEF      ,MFEM_UNDEF,MFEM_UNDEF},
233                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF   ,MFEM_UNDEF      ,MFEM_UNDEF,MFEM_UNDEF},
234                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE     ,MFEM_UNDEF,MFEM_UNDEF},
235                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF   ,MFEM_TETRAHEDRON,MFEM_PRISM,MFEM_CUBE } };
236 
237 MFEM_cid mfem_table_cid_unint[4][9] = { {MFEM_POINT,MFEM_UNDEF,MFEM_UNDEF  ,MFEM_UNDEF   ,MFEM_UNDEF      ,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_UNDEF},
238                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF   ,MFEM_UNDEF      ,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_UNDEF},
239                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE     ,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_UNDEF},
240                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF   ,MFEM_TETRAHEDRON,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_CUBE } };
241 
242 static PetscErrorCode DMPlexGetPointMFEMCellID_Internal(DM dm, DMLabel label, PetscInt minl, PetscInt p, PetscInt *mid, PetscInt *cid)
243 {
244   DMLabel        dlabel;
245   PetscInt       depth,csize,pdepth,dim;
246   PetscErrorCode ierr;
247 
248   PetscFunctionBegin;
249   ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr);
250   ierr = DMLabelGetValue(dlabel,p,&pdepth);CHKERRQ(ierr);
251   ierr = DMPlexGetConeSize(dm,p,&csize);CHKERRQ(ierr);
252   ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr);
253   ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
254   if (label) {
255     ierr = DMLabelGetValue(label,p,mid);CHKERRQ(ierr);
256     *mid = *mid - minl + 1; /* MFEM does not like negative markers */
257   } else *mid = 1;
258   if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
259     if (dim < 0 || dim > 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %D",dim);
260     if (csize > 8) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Found cone size %D for point %D",csize,p);
261     if (depth != 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Found depth %D for point %D. You should interpolate the mesh first",depth,p);
262     *cid = mfem_table_cid_unint[dim][csize];
263   } else {
264     if (csize > 6) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cone size %D for point %D",csize,p);
265     if (pdepth < 0 || pdepth > 3) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Depth %D for point %D",csize,p);
266     *cid = mfem_table_cid[pdepth][csize];
267   }
268   PetscFunctionReturn(0);
269 }
270 
271 static PetscErrorCode DMPlexGetPointMFEMVertexIDs_Internal(DM dm, PetscInt p, PetscSection csec, PetscInt *nv, PetscInt vids[])
272 {
273   PetscInt       dim,sdim,dof = 0,off = 0,i,q,vStart,vEnd,numPoints,*points = NULL;
274   PetscErrorCode ierr;
275 
276   PetscFunctionBegin;
277   ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr);
278   ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
279   sdim = dim;
280   if (csec) {
281     PetscInt sStart,sEnd;
282 
283     ierr = DMGetCoordinateDim(dm,&sdim);CHKERRQ(ierr);
284     ierr = PetscSectionGetChart(csec,&sStart,&sEnd);CHKERRQ(ierr);
285     ierr = PetscSectionGetOffset(csec,vStart,&off);CHKERRQ(ierr);
286     off  = off/sdim;
287     if (p >= sStart && p < sEnd) {
288       ierr = PetscSectionGetDof(csec,p,&dof);CHKERRQ(ierr);
289     }
290   }
291   if (!dof) {
292     ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr);
293     for (i=0,q=0;i<numPoints*2;i+= 2)
294       if ((points[i] >= vStart) && (points[i] < vEnd))
295         vids[q++] = points[i]-vStart+off;
296     ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr);
297   } else {
298     ierr = PetscSectionGetOffset(csec,p,&off);CHKERRQ(ierr);
299     ierr = PetscSectionGetDof(csec,p,&dof);CHKERRQ(ierr);
300     for (q=0;q<dof/sdim;q++) vids[q] = off/sdim + q;
301   }
302   *nv = q;
303   PetscFunctionReturn(0);
304 }
305 
306 static PetscErrorCode GLVisCreateFE(PetscFE femIn,char name[32],PetscFE *fem)
307 {
308   DM              K;
309   PetscSpace      P;
310   PetscDualSpace  Q;
311   PetscQuadrature q,fq;
312   PetscInt        dim,deg,dof;
313   DMPolytopeType  ptype;
314   PetscBool       isSimplex,isTensor;
315   PetscBool       continuity = PETSC_FALSE;
316   PetscDTNodeType nodeType   = PETSCDTNODES_GAUSSJACOBI;
317   PetscBool       endpoint   = PETSC_TRUE;
318   MPI_Comm        comm;
319   PetscErrorCode  ierr;
320 
321   PetscFunctionBegin;
322   comm = PetscObjectComm((PetscObject)femIn);
323   ierr = PetscFEGetBasisSpace(femIn,&P);CHKERRQ(ierr);
324   ierr = PetscFEGetDualSpace(femIn,&Q);CHKERRQ(ierr);
325   ierr = PetscDualSpaceGetDM(Q,&K);CHKERRQ(ierr);
326   ierr = DMGetDimension(K,&dim);CHKERRQ(ierr);
327   ierr = PetscSpaceGetDegree(P,&deg,NULL);CHKERRQ(ierr);
328   ierr = PetscSpaceGetNumComponents(P,&dof);CHKERRQ(ierr);
329   ierr = DMPlexGetCellType(K,0,&ptype);CHKERRQ(ierr);
330   switch (ptype) {
331   case DM_POLYTOPE_QUADRILATERAL:
332   case DM_POLYTOPE_HEXAHEDRON:
333     isSimplex = PETSC_FALSE; break;
334   default:
335     isSimplex = PETSC_TRUE; break;
336   }
337   isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
338   /* Create space */
339   ierr = PetscSpaceCreate(comm,&P);CHKERRQ(ierr);
340   ierr = PetscSpaceSetType(P,PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr);
341   ierr = PetscSpacePolynomialSetTensor(P,isTensor);CHKERRQ(ierr);
342   ierr = PetscSpaceSetNumComponents(P,dof);CHKERRQ(ierr);
343   ierr = PetscSpaceSetNumVariables(P,dim);CHKERRQ(ierr);
344   ierr = PetscSpaceSetDegree(P,deg,PETSC_DETERMINE);CHKERRQ(ierr);
345   ierr = PetscSpaceSetUp(P);CHKERRQ(ierr);
346   /* Create dual space */
347   ierr = PetscDualSpaceCreate(comm,&Q);CHKERRQ(ierr);
348   ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr);
349   ierr = PetscDualSpaceLagrangeSetTensor(Q,isTensor);CHKERRQ(ierr);
350   ierr = PetscDualSpaceLagrangeSetContinuity(Q,continuity);CHKERRQ(ierr);
351   ierr = PetscDualSpaceLagrangeSetNodeType(Q,nodeType,endpoint,0);CHKERRQ(ierr);
352   ierr = PetscDualSpaceSetNumComponents(Q,dof);CHKERRQ(ierr);
353   ierr = PetscDualSpaceSetOrder(Q,deg);CHKERRQ(ierr);
354   ierr = PetscDualSpaceCreateReferenceCell(Q,dim,isSimplex,&K);CHKERRQ(ierr);
355   ierr = PetscDualSpaceSetDM(Q,K);CHKERRQ(ierr);
356   ierr = DMDestroy(&K);CHKERRQ(ierr);
357   ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr);
358   /* Create quadrature */
359   if (isSimplex) {
360     ierr = PetscDTStroudConicalQuadrature(dim,  1,deg+1,-1,+1,&q);CHKERRQ(ierr);
361     ierr = PetscDTStroudConicalQuadrature(dim-1,1,deg+1,-1,+1,&fq);CHKERRQ(ierr);
362   } else {
363     ierr = PetscDTGaussTensorQuadrature(dim,  1,deg+1,-1,+1,&q);CHKERRQ(ierr);
364     ierr = PetscDTGaussTensorQuadrature(dim-1,1,deg+1,-1,+1,&fq);CHKERRQ(ierr);
365   }
366   /* Create finite element */
367   ierr = PetscFECreate(comm,fem);CHKERRQ(ierr);
368   ierr = PetscSNPrintf(name,32,"L2_T1_%DD_P%D",dim,deg);CHKERRQ(ierr);
369   ierr = PetscObjectSetName((PetscObject)*fem,name);CHKERRQ(ierr);
370   ierr = PetscFESetType(*fem,PETSCFEBASIC);CHKERRQ(ierr);
371   ierr = PetscFESetNumComponents(*fem,dof);CHKERRQ(ierr);
372   ierr = PetscFESetBasisSpace(*fem,P);CHKERRQ(ierr);
373   ierr = PetscFESetDualSpace(*fem,Q);CHKERRQ(ierr);
374   ierr = PetscFESetQuadrature(*fem,q);CHKERRQ(ierr);
375   ierr = PetscFESetFaceQuadrature(*fem,fq);CHKERRQ(ierr);
376   ierr = PetscFESetUp(*fem);CHKERRQ(ierr);
377   /* Cleanup */
378   ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr);
379   ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr);
380   ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr);
381   ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr);
382   PetscFunctionReturn(0);
383 }
384 
385 /*
386    ASCII visualization/dump: full support for simplices and tensor product cells. It supports AMR
387    Higher order meshes are also supported
388 */
389 static PetscErrorCode DMPlexView_GLVis_ASCII(DM dm, PetscViewer viewer)
390 {
391   DMLabel              label;
392   PetscSection         coordSection,parentSection;
393   Vec                  coordinates,hovec;
394   const PetscScalar    *array;
395   PetscInt             bf,p,sdim,dim,depth,novl,minl;
396   PetscInt             cStart,cEnd,vStart,vEnd,nvert;
397   PetscMPIInt          size;
398   PetscBool            localized,isascii;
399   PetscBool            enable_mfem,enable_boundary,enable_ncmesh,view_ovl = PETSC_FALSE;
400   PetscBT              pown,vown;
401   PetscErrorCode       ierr;
402   PetscContainer       glvis_container;
403   PetscBool            cellvertex = PETSC_FALSE, periodic, enabled = PETSC_TRUE;
404   PetscBool            enable_emark,enable_bmark;
405   const char           *fmt;
406   char                 emark[64] = "",bmark[64] = "";
407 
408   PetscFunctionBegin;
409   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
410   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
411   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
412   if (!isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERASCII");
413   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)viewer),&size);CHKERRQ(ierr);
414   if (size > 1) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Use single sequential viewers for parallel visualization");
415   ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
416 
417   /* get container: determines if a process visualizes is portion of the data or not */
418   ierr = PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);CHKERRQ(ierr);
419   if (!glvis_container) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing GLVis container");
420   {
421     PetscViewerGLVisInfo glvis_info;
422     ierr    = PetscContainerGetPointer(glvis_container,(void**)&glvis_info);CHKERRQ(ierr);
423     enabled = glvis_info->enabled;
424     fmt     = glvis_info->fmt;
425   }
426 
427   /* Users can attach a coordinate vector to the DM in case they have a higher-order mesh
428      DMPlex does not currently support HO meshes, so there's no API for this */
429   ierr = PetscObjectQuery((PetscObject)dm,"_glvis_mesh_coords",(PetscObject*)&hovec);CHKERRQ(ierr);
430   ierr = PetscObjectReference((PetscObject)hovec);CHKERRQ(ierr);
431   if (!hovec) {
432     DM           cdm;
433     PetscFE      disc;
434     PetscClassId classid;
435 
436     ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr);
437     ierr = DMGetField(cdm,0,NULL,(PetscObject*)&disc);CHKERRQ(ierr);
438     ierr = PetscObjectGetClassId((PetscObject)disc,&classid);CHKERRQ(ierr);
439     if (classid == PETSCFE_CLASSID) {
440       DM      hocdm;
441       PetscFE hodisc;
442       Vec     vec;
443       Mat     mat;
444       char    name[32],fec_type[64];
445 
446       ierr = GLVisCreateFE(disc,name,&hodisc);CHKERRQ(ierr);
447       ierr = DMClone(cdm,&hocdm);CHKERRQ(ierr);
448       ierr = DMSetField(hocdm,0,NULL,(PetscObject)hodisc);CHKERRQ(ierr);
449       ierr = PetscFEDestroy(&hodisc);CHKERRQ(ierr);
450       ierr = DMCreateDS(hocdm);CHKERRQ(ierr);
451 
452       ierr = DMGetCoordinates(dm,&vec);CHKERRQ(ierr);
453       ierr = DMCreateGlobalVector(hocdm,&hovec);CHKERRQ(ierr);
454       ierr = DMCreateInterpolation(cdm,hocdm,&mat,NULL);CHKERRQ(ierr);
455       ierr = MatInterpolate(mat,vec,hovec);CHKERRQ(ierr);
456       ierr = MatDestroy(&mat);CHKERRQ(ierr);
457       ierr = DMDestroy(&hocdm);CHKERRQ(ierr);
458 
459       ierr = PetscSNPrintf(fec_type,sizeof(fec_type),"FiniteElementCollection: %s", name);CHKERRQ(ierr);
460       ierr = PetscObjectSetName((PetscObject)hovec,fec_type);CHKERRQ(ierr);
461     }
462   }
463 
464   ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr);
465   ierr = DMPlexGetGhostCellStratum(dm,&p,NULL);CHKERRQ(ierr);
466   if (p >= 0) cEnd = p;
467   ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr);
468   ierr = DMGetPeriodicity(dm,&periodic,NULL,NULL,NULL);CHKERRQ(ierr);
469   ierr = DMGetCoordinatesLocalized(dm,&localized);CHKERRQ(ierr);
470   if (periodic && !localized && !hovec) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Coordinates need to be localized");
471   ierr = DMGetCoordinateSection(dm,&coordSection);CHKERRQ(ierr);
472   ierr = DMGetCoordinateDim(dm,&sdim);CHKERRQ(ierr);
473   ierr = DMGetCoordinatesLocal(dm,&coordinates);CHKERRQ(ierr);
474   if (!coordinates && !hovec) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
475 
476   /*
477      a couple of sections of the mesh specification are disabled
478        - boundary: the boundary is not needed for proper mesh visualization unless we want to visualize boundary attributes or we have high-order coordinates in 3D (topologically)
479        - vertex_parents: used for non-conforming meshes only when we want to use MFEM as a discretization package
480                          and be able to derefine the mesh (MFEM does not currently have to ability to read ncmeshes in parallel)
481   */
482   enable_boundary = PETSC_FALSE;
483   enable_ncmesh   = PETSC_FALSE;
484   enable_mfem     = PETSC_FALSE;
485   enable_emark    = PETSC_FALSE;
486   enable_bmark    = PETSC_FALSE;
487   /* I'm tired of problems with negative values in the markers, disable them */
488   ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)dm),((PetscObject)dm)->prefix,"GLVis PetscViewer DMPlex Options","PetscViewer");CHKERRQ(ierr);
489   ierr = PetscOptionsBool("-viewer_glvis_dm_plex_enable_boundary","Enable boundary section in mesh representation",NULL,enable_boundary,&enable_boundary,NULL);CHKERRQ(ierr);
490   ierr = PetscOptionsBool("-viewer_glvis_dm_plex_enable_ncmesh","Enable vertex_parents section in mesh representation (allows derefinement)",NULL,enable_ncmesh,&enable_ncmesh,NULL);CHKERRQ(ierr);
491   ierr = PetscOptionsBool("-viewer_glvis_dm_plex_enable_mfem","Dump a mesh that can be used with MFEM's FiniteElementSpaces",NULL,enable_mfem,&enable_mfem,NULL);CHKERRQ(ierr);
492   ierr = PetscOptionsBool("-viewer_glvis_dm_plex_overlap","Include overlap region in local meshes",NULL,view_ovl,&view_ovl,NULL);CHKERRQ(ierr);
493   ierr = PetscOptionsString("-viewer_glvis_dm_plex_emarker","String for the material id label",NULL,emark,emark,sizeof(emark),&enable_emark);CHKERRQ(ierr);
494   ierr = PetscOptionsString("-viewer_glvis_dm_plex_bmarker","String for the boundary id label",NULL,bmark,bmark,sizeof(bmark),&enable_bmark);CHKERRQ(ierr);
495   ierr = PetscOptionsEnd();CHKERRQ(ierr);
496   if (enable_bmark) enable_boundary = PETSC_TRUE;
497 
498   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);CHKERRQ(ierr);
499   if (enable_ncmesh && size > 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not supported in parallel");
500   ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr);
501   if (enable_boundary && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. "
502                                                              "Alternatively, run with -viewer_glvis_dm_plex_enable_boundary 0");
503   if (enable_ncmesh && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. "
504                                                            "Alternatively, run with -viewer_glvis_dm_plex_enable_ncmesh 0");
505   if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
506     if (depth != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported depth %D. You should interpolate the mesh first",depth);
507     cellvertex = PETSC_TRUE;
508   }
509 
510   /* Identify possible cells in the overlap */
511   novl = 0;
512   pown = NULL;
513   if (size > 1) {
514     IS             globalNum = NULL;
515     const PetscInt *gNum;
516     PetscBool      ovl  = PETSC_FALSE;
517 
518     ierr = PetscObjectQuery((PetscObject)dm,"_glvis_plex_gnum",(PetscObject*)&globalNum);CHKERRQ(ierr);
519     if (!globalNum) {
520       if (view_ovl) {
521         ierr = ISCreateStride(PetscObjectComm((PetscObject)dm),cEnd-cStart,0,1,&globalNum);CHKERRQ(ierr);
522       } else {
523         ierr = DMPlexCreateCellNumbering_Internal(dm,PETSC_TRUE,&globalNum);CHKERRQ(ierr);
524       }
525       ierr = PetscObjectCompose((PetscObject)dm,"_glvis_plex_gnum",(PetscObject)globalNum);CHKERRQ(ierr);
526       ierr = PetscObjectDereference((PetscObject)globalNum);CHKERRQ(ierr);
527     }
528     ierr = ISGetIndices(globalNum,&gNum);CHKERRQ(ierr);
529     for (p=cStart; p<cEnd; p++) {
530       if (gNum[p-cStart] < 0) {
531         ovl = PETSC_TRUE;
532         novl++;
533       }
534     }
535     if (ovl) {
536       /* it may happen that pown get not destroyed, if the user closes the window while this function is running.
537          TODO: garbage collector? attach pown to dm?  */
538       ierr = PetscBTCreate(cEnd-cStart,&pown);CHKERRQ(ierr);
539       for (p=cStart; p<cEnd; p++) {
540         if (gNum[p-cStart] < 0) continue;
541         else {
542           ierr = PetscBTSet(pown,p-cStart);CHKERRQ(ierr);
543         }
544       }
545     }
546     ierr = ISRestoreIndices(globalNum,&gNum);CHKERRQ(ierr);
547   }
548 
549   /* return if this process is disabled */
550   if (!enabled) {
551     ierr = PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");CHKERRQ(ierr);
552     ierr = PetscViewerASCIIPrintf(viewer,"\ndimension\n");CHKERRQ(ierr);
553     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",dim);CHKERRQ(ierr);
554     ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr);
555     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr);
556     ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr);
557     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr);
558     ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr);
559     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr);
560     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr);
561     ierr = PetscBTDestroy(&pown);CHKERRQ(ierr);
562     ierr = VecDestroy(&hovec);CHKERRQ(ierr);
563     PetscFunctionReturn(0);
564   }
565 
566   if (enable_mfem) {
567     if (periodic && !hovec) { /* we need to generate a vector of L2 coordinates, as this is how MFEM handles periodic meshes */
568       PetscInt    vpc = 0;
569       char        fec[64];
570       PetscInt    vids[8] = {0,1,2,3,4,5,6,7};
571       PetscInt    hexv[8] = {0,1,3,2,4,5,7,6}, tetv[4] = {0,1,2,3};
572       PetscInt    quadv[8] = {0,1,3,2}, triv[3] = {0,1,2};
573       PetscInt    *dof = NULL;
574       PetscScalar *array,*ptr;
575 
576       ierr = PetscSNPrintf(fec,sizeof(fec),"FiniteElementCollection: L2_T1_%DD_P1",dim);CHKERRQ(ierr);
577       if (cEnd-cStart) {
578         PetscInt fpc;
579 
580         ierr = DMPlexGetConeSize(dm,cStart,&fpc);CHKERRQ(ierr);
581         switch(dim) {
582           case 1:
583             vpc = 2;
584             dof = hexv;
585             break;
586           case 2:
587             switch (fpc) {
588               case 3:
589                 vpc = 3;
590                 dof = triv;
591                 break;
592               case 4:
593                 vpc = 4;
594                 dof = quadv;
595                 break;
596               default:
597                 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
598                 break;
599             }
600             break;
601           case 3:
602             switch (fpc) {
603               case 4: /* TODO: still need to understand L2 ordering for tets */
604                 vpc = 4;
605                 dof = tetv;
606                 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled tethraedral case");
607                 break;
608               case 6:
609                 if (cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: vertices per cell %D",fpc);
610                 vpc = 8;
611                 dof = hexv;
612                 break;
613               case 8:
614                 if (!cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
615                 vpc = 8;
616                 dof = hexv;
617                 break;
618               default:
619                 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
620                 break;
621             }
622             break;
623           default:
624             SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim");
625             break;
626         }
627         ierr = DMPlexReorderCell(dm,cStart,vids);CHKERRQ(ierr);
628       }
629       if (!dof) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing dofs");
630       ierr = VecCreateSeq(PETSC_COMM_SELF,(cEnd-cStart-novl)*vpc*sdim,&hovec);CHKERRQ(ierr);
631       ierr = PetscObjectSetName((PetscObject)hovec,fec);CHKERRQ(ierr);
632       ierr = VecGetArray(hovec,&array);CHKERRQ(ierr);
633       ptr  = array;
634       for (p=cStart;p<cEnd;p++) {
635         PetscInt    csize,v,d;
636         PetscScalar *vals = NULL;
637 
638         if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue;
639         ierr = DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr);
640         if (csize != vpc*sdim && csize != vpc*sdim*2) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported closure size %D (vpc %D, sdim %D)",csize,vpc,sdim);
641         for (v=0;v<vpc;v++) {
642           for (d=0;d<sdim;d++) {
643             ptr[sdim*dof[v]+d] = vals[sdim*vids[v]+d];
644           }
645         }
646         ptr += vpc*sdim;
647         ierr = DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr);
648       }
649       ierr = VecRestoreArray(hovec,&array);CHKERRQ(ierr);
650     }
651   }
652   /* if we have high-order coordinates in 3D, we need to specify the boundary */
653   if (hovec && dim == 3) enable_boundary = PETSC_TRUE;
654 
655   /* header */
656   ierr = PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");CHKERRQ(ierr);
657 
658   /* topological dimension */
659   ierr = PetscViewerASCIIPrintf(viewer,"\ndimension\n");CHKERRQ(ierr);
660   ierr = PetscViewerASCIIPrintf(viewer,"%D\n",dim);CHKERRQ(ierr);
661 
662   /* elements */
663   minl = 1;
664   label = NULL;
665   if (enable_emark) {
666     PetscInt lminl = PETSC_MAX_INT;
667 
668     ierr = DMGetLabel(dm,emark,&label);CHKERRQ(ierr);
669     if (label) {
670       IS       vals;
671       PetscInt ldef;
672 
673       ierr = DMLabelGetDefaultValue(label,&ldef);CHKERRQ(ierr);
674       ierr = DMLabelGetValueIS(label,&vals);CHKERRQ(ierr);
675       ierr = ISGetMinMax(vals,&lminl,NULL);CHKERRQ(ierr);
676       ierr = ISDestroy(&vals);CHKERRQ(ierr);
677       lminl = PetscMin(ldef,lminl);
678     }
679     ierr = MPIU_Allreduce(&lminl,&minl,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
680     if (minl == PETSC_MAX_INT) minl = 1;
681   }
682   ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr);
683   ierr = PetscViewerASCIIPrintf(viewer,"%D\n",cEnd-cStart-novl);CHKERRQ(ierr);
684   for (p=cStart;p<cEnd;p++) {
685     PetscInt       vids[8];
686     PetscInt       i,nv = 0,cid = -1,mid = 1;
687 
688     if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue;
689     ierr = DMPlexGetPointMFEMCellID_Internal(dm,label,minl,p,&mid,&cid);CHKERRQ(ierr);
690     ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,p,(localized && !hovec) ? coordSection : NULL,&nv,vids);CHKERRQ(ierr);
691     ierr = DMPlexReorderCell(dm,p,vids);CHKERRQ(ierr);
692     ierr = PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);CHKERRQ(ierr);
693     for (i=0;i<nv;i++) {
694       ierr = PetscViewerASCIIPrintf(viewer," %D",vids[i]);CHKERRQ(ierr);
695     }
696     ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
697   }
698 
699   /* boundary */
700   ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr);
701   if (!enable_boundary) {
702     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr);
703   } else {
704     DMLabel  perLabel;
705     PetscBT  bfaces;
706     PetscInt fStart,fEnd,*fcells;
707 
708     ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd);CHKERRQ(ierr);
709     ierr = PetscBTCreate(fEnd-fStart,&bfaces);CHKERRQ(ierr);
710     ierr = DMPlexGetMaxSizes(dm,NULL,&p);CHKERRQ(ierr);
711     ierr = PetscMalloc1(p,&fcells);CHKERRQ(ierr);
712     ierr = DMGetLabel(dm,"glvis_periodic_cut",&perLabel);CHKERRQ(ierr);
713     if (!perLabel && localized) { /* this periodic cut can be moved up to DMPlex setup */
714       ierr = DMCreateLabel(dm,"glvis_periodic_cut");CHKERRQ(ierr);
715       ierr = DMGetLabel(dm,"glvis_periodic_cut",&perLabel);CHKERRQ(ierr);
716       ierr = DMLabelSetDefaultValue(perLabel,1);CHKERRQ(ierr);
717       for (p=cStart;p<cEnd;p++) {
718         DMPolytopeType cellType;
719         PetscInt       dof;
720 
721         ierr = DMPlexGetCellType(dm,p,&cellType);CHKERRQ(ierr);
722         ierr = PetscSectionGetDof(coordSection,p,&dof);CHKERRQ(ierr);
723         if (dof) {
724           PetscInt    uvpc, v,csize,cellClosureSize,*cellClosure = NULL,*vidxs = NULL;
725           PetscScalar *vals = NULL;
726 
727           uvpc = DMPolytopeTypeGetNumVertices(cellType);
728           if (dof%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Incompatible number of cell dofs %D and space dimension %D",dof,sdim);
729           ierr = DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr);
730           ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr);
731           for (v=0;v<cellClosureSize;v++)
732             if (cellClosure[2*v] >= vStart && cellClosure[2*v] < vEnd) {
733               vidxs = cellClosure + 2*v;
734               break;
735             }
736           if (!vidxs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing vertices");
737           for (v=0;v<uvpc;v++) {
738             PetscInt s;
739 
740             for (s=0;s<sdim;s++) {
741               if (PetscAbsScalar(vals[v*sdim+s]-vals[v*sdim+s+uvpc*sdim])>PETSC_MACHINE_EPSILON) {
742                 ierr = DMLabelSetValue(perLabel,vidxs[2*v],2);CHKERRQ(ierr);
743               }
744             }
745           }
746           ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr);
747           ierr = DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr);
748         }
749       }
750       if (dim > 1) {
751         PetscInt eEnd,eStart;
752 
753         ierr = DMPlexGetDepthStratum(dm,1,&eStart,&eEnd);CHKERRQ(ierr);
754         for (p=eStart;p<eEnd;p++) {
755           const PetscInt *cone;
756           PetscInt       coneSize,i;
757           PetscBool      ispe = PETSC_TRUE;
758 
759           ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr);
760           ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr);
761           for (i=0;i<coneSize;i++) {
762             PetscInt v;
763 
764             ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr);
765             ispe = (PetscBool)(ispe && (v==2));
766           }
767           if (ispe && coneSize) {
768             PetscInt       ch, numChildren;
769             const PetscInt *children;
770 
771             ierr = DMLabelSetValue(perLabel,p,2);CHKERRQ(ierr);
772             ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr);
773             for (ch = 0; ch < numChildren; ch++) {
774               ierr = DMLabelSetValue(perLabel,children[ch],2);CHKERRQ(ierr);
775             }
776           }
777         }
778         if (dim > 2) {
779           for (p=fStart;p<fEnd;p++) {
780             const PetscInt *cone;
781             PetscInt       coneSize,i;
782             PetscBool      ispe = PETSC_TRUE;
783 
784             ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr);
785             ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr);
786             for (i=0;i<coneSize;i++) {
787               PetscInt v;
788 
789               ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr);
790               ispe = (PetscBool)(ispe && (v==2));
791             }
792             if (ispe && coneSize) {
793               PetscInt       ch, numChildren;
794               const PetscInt *children;
795 
796               ierr = DMLabelSetValue(perLabel,p,2);CHKERRQ(ierr);
797               ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr);
798               for (ch = 0; ch < numChildren; ch++) {
799                 ierr = DMLabelSetValue(perLabel,children[ch],2);CHKERRQ(ierr);
800               }
801             }
802           }
803         }
804       }
805     }
806     for (p=fStart;p<fEnd;p++) {
807       const PetscInt *support;
808       PetscInt       supportSize;
809       PetscBool      isbf = PETSC_FALSE;
810 
811       ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr);
812       if (pown) {
813         PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE;
814         PetscInt  i;
815 
816         ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr);
817         for (i=0;i<supportSize;i++) {
818           if (PetscLikely(PetscBTLookup(pown,support[i]-cStart))) has_owned = PETSC_TRUE;
819           else has_ghost = PETSC_TRUE;
820         }
821         isbf = (PetscBool)((supportSize == 1 && has_owned) || (supportSize > 1 && has_owned && has_ghost));
822       } else {
823         isbf = (PetscBool)(supportSize == 1);
824       }
825       if (!isbf && perLabel) {
826         const PetscInt *cone;
827         PetscInt       coneSize,i;
828 
829         ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr);
830         ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr);
831         isbf = PETSC_TRUE;
832         for (i=0;i<coneSize;i++) {
833           PetscInt v,d;
834 
835           ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr);
836           ierr = DMLabelGetDefaultValue(perLabel,&d);CHKERRQ(ierr);
837           isbf = (PetscBool)(isbf && v != d);
838         }
839       }
840       if (isbf) {
841         ierr = PetscBTSet(bfaces,p-fStart);CHKERRQ(ierr);
842       }
843     }
844     /* count boundary faces */
845     for (p=fStart,bf=0;p<fEnd;p++) {
846       if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) {
847         const PetscInt *support;
848         PetscInt       supportSize,c;
849 
850         ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr);
851         ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr);
852         for (c=0;c<supportSize;c++) {
853           const    PetscInt *cone;
854           PetscInt cell,cl,coneSize;
855 
856           cell = support[c];
857           if (pown && PetscUnlikely(!PetscBTLookup(pown,cell-cStart))) continue;
858           ierr = DMPlexGetCone(dm,cell,&cone);CHKERRQ(ierr);
859           ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr);
860           for (cl=0;cl<coneSize;cl++) {
861             if (cone[cl] == p) {
862               bf += 1;
863               break;
864             }
865           }
866         }
867       }
868     }
869     minl = 1;
870     label = NULL;
871     if (enable_bmark) {
872       PetscInt lminl = PETSC_MAX_INT;
873 
874       ierr = DMGetLabel(dm,bmark,&label);CHKERRQ(ierr);
875       if (label) {
876         IS       vals;
877         PetscInt ldef;
878 
879         ierr = DMLabelGetDefaultValue(label,&ldef);CHKERRQ(ierr);
880         ierr = DMLabelGetValueIS(label,&vals);CHKERRQ(ierr);
881         ierr = ISGetMinMax(vals,&lminl,NULL);CHKERRQ(ierr);
882         ierr = ISDestroy(&vals);CHKERRQ(ierr);
883         lminl = PetscMin(ldef,lminl);
884       }
885       ierr = MPIU_Allreduce(&lminl,&minl,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
886       if (minl == PETSC_MAX_INT) minl = 1;
887     }
888     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",bf);CHKERRQ(ierr);
889     for (p=fStart;p<fEnd;p++) {
890       if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) {
891         const PetscInt *support;
892         PetscInt       supportSize,c,nc = 0;
893 
894         ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr);
895         ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr);
896         if (pown) {
897           for (c=0;c<supportSize;c++) {
898             if (PetscLikely(PetscBTLookup(pown,support[c]-cStart))) {
899               fcells[nc++] = support[c];
900             }
901           }
902         } else for (c=0;c<supportSize;c++) fcells[nc++] = support[c];
903         for (c=0;c<nc;c++) {
904           const DMPolytopeType *faceTypes;
905           DMPolytopeType       cellType;
906           const PetscInt       *faceSizes,*cone;
907           PetscInt             vids[8],*faces,st,i,coneSize,cell,cl,nv,cid = -1,mid = -1;
908 
909           cell = fcells[c];
910           ierr = DMPlexGetCone(dm,cell,&cone);CHKERRQ(ierr);
911           ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr);
912           for (cl=0;cl<coneSize;cl++)
913             if (cone[cl] == p)
914               break;
915           if (cl == coneSize) continue;
916 
917           /* face material id and type */
918           ierr = DMPlexGetPointMFEMCellID_Internal(dm,label,minl,p,&mid,&cid);CHKERRQ(ierr);
919           ierr = PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);CHKERRQ(ierr);
920           /* vertex ids */
921           ierr = DMPlexGetCellType(dm,cell,&cellType);CHKERRQ(ierr);
922           ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,cell,(localized && !hovec) ? coordSection : NULL,&nv,vids);CHKERRQ(ierr);
923           ierr = DMPlexGetRawFaces_Internal(dm,cellType,vids,NULL,&faceTypes,&faceSizes,(const PetscInt**)&faces);CHKERRQ(ierr);
924           st = 0;
925           for (i=0;i<cl;i++) st += faceSizes[i];
926           ierr = DMPlexInvertCell(faceTypes[cl],faces + st);CHKERRQ(ierr);
927           for (i=0;i<faceSizes[cl];i++) {
928             ierr = PetscViewerASCIIPrintf(viewer," %d",faces[st+i]);CHKERRQ(ierr);
929           }
930           ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
931           ierr = DMPlexRestoreRawFaces_Internal(dm,cellType,vids,NULL,&faceTypes,&faceSizes,(const PetscInt**)&faces);CHKERRQ(ierr);
932           bf -= 1;
933         }
934       }
935     }
936     if (bf) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Remaining boundary faces %D",bf);
937     ierr = PetscBTDestroy(&bfaces);CHKERRQ(ierr);
938     ierr = PetscFree(fcells);CHKERRQ(ierr);
939   }
940 
941   /* mark owned vertices */
942   vown = NULL;
943   if (pown) {
944     ierr = PetscBTCreate(vEnd-vStart,&vown);CHKERRQ(ierr);
945     for (p=cStart;p<cEnd;p++) {
946       PetscInt i,closureSize,*closure = NULL;
947 
948       if (PetscUnlikely(!PetscBTLookup(pown,p-cStart))) continue;
949       ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
950       for (i=0;i<closureSize;i++) {
951         const PetscInt pp = closure[2*i];
952 
953         if (pp >= vStart && pp < vEnd) {
954           ierr = PetscBTSet(vown,pp-vStart);CHKERRQ(ierr);
955         }
956       }
957       ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
958     }
959   }
960 
961   /* vertex_parents (Non-conforming meshes) */
962   parentSection  = NULL;
963   if (enable_ncmesh) {
964     ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
965   }
966   if (parentSection) {
967     PetscInt vp,gvp;
968 
969     for (vp=0,p=vStart;p<vEnd;p++) {
970       DMLabel  dlabel;
971       PetscInt parent,depth;
972 
973       if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue;
974       ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr);
975       ierr = DMLabelGetValue(dlabel,p,&depth);CHKERRQ(ierr);
976       ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr);
977       if (parent != p) vp++;
978     }
979     ierr = MPIU_Allreduce(&vp,&gvp,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
980     if (gvp) {
981       PetscInt  maxsupp;
982       PetscBool *skip = NULL;
983 
984       ierr = PetscViewerASCIIPrintf(viewer,"\nvertex_parents\n");CHKERRQ(ierr);
985       ierr = PetscViewerASCIIPrintf(viewer,"%D\n",vp);CHKERRQ(ierr);
986       ierr = DMPlexGetMaxSizes(dm,NULL,&maxsupp);CHKERRQ(ierr);
987       ierr = PetscMalloc1(maxsupp,&skip);CHKERRQ(ierr);
988       for (p=vStart;p<vEnd;p++) {
989         DMLabel  dlabel;
990         PetscInt parent;
991 
992         if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue;
993         ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr);
994         ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr);
995         if (parent != p) {
996           PetscInt       vids[8] = { -1, -1, -1, -1, -1, -1, -1, -1 }; /* silent overzealous clang static analyzer */
997           PetscInt       i,nv,ssize,n,numChildren,depth = -1;
998           const PetscInt *children;
999 
1000           ierr = DMPlexGetConeSize(dm,parent,&ssize);CHKERRQ(ierr);
1001           switch (ssize) {
1002             case 2: /* edge */
1003               nv   = 0;
1004               ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,parent,localized ? coordSection : NULL,&nv,vids);CHKERRQ(ierr);
1005               ierr = PetscViewerASCIIPrintf(viewer,"%D",p-vStart);CHKERRQ(ierr);
1006               for (i=0;i<nv;i++) {
1007                 ierr = PetscViewerASCIIPrintf(viewer," %D",vids[i]);CHKERRQ(ierr);
1008               }
1009               ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
1010               vp--;
1011               break;
1012             case 4: /* face */
1013               ierr = DMPlexGetTreeChildren(dm,parent,&numChildren,&children);CHKERRQ(ierr);
1014               for (n=0;n<numChildren;n++) {
1015                 ierr = DMLabelGetValue(dlabel,children[n],&depth);CHKERRQ(ierr);
1016                 if (!depth) {
1017                   const PetscInt *hvsupp,*hesupp,*cone;
1018                   PetscInt       hvsuppSize,hesuppSize,coneSize;
1019                   PetscInt       hv = children[n],he = -1,f;
1020 
1021                   ierr = PetscArrayzero(skip,maxsupp);CHKERRQ(ierr);
1022                   ierr = DMPlexGetSupportSize(dm,hv,&hvsuppSize);CHKERRQ(ierr);
1023                   ierr = DMPlexGetSupport(dm,hv,&hvsupp);CHKERRQ(ierr);
1024                   for (i=0;i<hvsuppSize;i++) {
1025                     PetscInt ep;
1026                     ierr = DMPlexGetTreeParent(dm,hvsupp[i],&ep,NULL);CHKERRQ(ierr);
1027                     if (ep != hvsupp[i]) {
1028                       he = hvsupp[i];
1029                     } else {
1030                       skip[i] = PETSC_TRUE;
1031                     }
1032                   }
1033                   if (he == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vertex %D support size %D: hanging edge not found",hv,hvsuppSize);
1034                   ierr    = DMPlexGetCone(dm,he,&cone);CHKERRQ(ierr);
1035                   vids[0] = (cone[0] == hv) ? cone[1] : cone[0];
1036                   ierr    = DMPlexGetSupportSize(dm,he,&hesuppSize);CHKERRQ(ierr);
1037                   ierr    = DMPlexGetSupport(dm,he,&hesupp);CHKERRQ(ierr);
1038                   for (f=0;f<hesuppSize;f++) {
1039                     PetscInt j;
1040 
1041                     ierr = DMPlexGetCone(dm,hesupp[f],&cone);CHKERRQ(ierr);
1042                     ierr = DMPlexGetConeSize(dm,hesupp[f],&coneSize);CHKERRQ(ierr);
1043                     for (j=0;j<coneSize;j++) {
1044                       PetscInt k;
1045                       for (k=0;k<hvsuppSize;k++) {
1046                         if (hvsupp[k] == cone[j]) {
1047                           skip[k] = PETSC_TRUE;
1048                           break;
1049                         }
1050                       }
1051                     }
1052                   }
1053                   for (i=0;i<hvsuppSize;i++) {
1054                     if (!skip[i]) {
1055                       ierr = DMPlexGetCone(dm,hvsupp[i],&cone);CHKERRQ(ierr);
1056                       vids[1] = (cone[0] == hv) ? cone[1] : cone[0];
1057                     }
1058                   }
1059                   ierr = PetscViewerASCIIPrintf(viewer,"%D",hv-vStart);CHKERRQ(ierr);
1060                   for (i=0;i<2;i++) {
1061                     ierr = PetscViewerASCIIPrintf(viewer," %D",vids[i]-vStart);CHKERRQ(ierr);
1062                   }
1063                   ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
1064                   vp--;
1065                 }
1066               }
1067               break;
1068             default:
1069               SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Don't know how to deal with support size %D",ssize);
1070           }
1071         }
1072       }
1073       ierr = PetscFree(skip);CHKERRQ(ierr);
1074     }
1075     if (vp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %D hanging vertices",vp);
1076   }
1077   ierr = PetscBTDestroy(&pown);CHKERRQ(ierr);
1078   ierr = PetscBTDestroy(&vown);CHKERRQ(ierr);
1079 
1080   /* vertices */
1081   if (hovec) { /* higher-order meshes */
1082     const char *fec;
1083     PetscInt   i,n,s;
1084 
1085     ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr);
1086     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",vEnd-vStart);CHKERRQ(ierr);
1087     ierr = PetscViewerASCIIPrintf(viewer,"nodes\n");CHKERRQ(ierr);
1088     ierr = PetscObjectGetName((PetscObject)hovec,&fec);CHKERRQ(ierr);
1089     ierr = PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");CHKERRQ(ierr);
1090     ierr = PetscViewerASCIIPrintf(viewer,"%s\n",fec);CHKERRQ(ierr);
1091     ierr = PetscViewerASCIIPrintf(viewer,"VDim: %D\n",sdim);CHKERRQ(ierr);
1092     ierr = PetscViewerASCIIPrintf(viewer,"Ordering: 1\n\n");CHKERRQ(ierr); /*Ordering::byVDIM*/
1093     ierr = VecGetArrayRead(hovec,&array);CHKERRQ(ierr);
1094     ierr = VecGetLocalSize(hovec,&n);CHKERRQ(ierr);
1095     if (n%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Size of local coordinate vector %D incompatible with space dimension %D",n,sdim);
1096     for (i=0;i<n/sdim;i++) {
1097       for (s=0;s<sdim;s++) {
1098         ierr = PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[i*sdim+s]));CHKERRQ(ierr);
1099       }
1100       ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
1101     }
1102     ierr = VecRestoreArrayRead(hovec,&array);CHKERRQ(ierr);
1103   } else {
1104     ierr = VecGetLocalSize(coordinates,&nvert);CHKERRQ(ierr);
1105     ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr);
1106     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",nvert/sdim);CHKERRQ(ierr);
1107     ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr);
1108     ierr = VecGetArrayRead(coordinates,&array);CHKERRQ(ierr);
1109     for (p=0;p<nvert/sdim;p++) {
1110       PetscInt s;
1111       for (s=0;s<sdim;s++) {
1112         PetscReal v = PetscRealPart(array[p*sdim+s]);
1113 
1114         ierr = PetscViewerASCIIPrintf(viewer,fmt,PetscIsInfOrNanReal(v) ? 0.0 : v);CHKERRQ(ierr);
1115       }
1116       ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
1117     }
1118     ierr = VecRestoreArrayRead(coordinates,&array);CHKERRQ(ierr);
1119   }
1120   ierr = VecDestroy(&hovec);CHKERRQ(ierr);
1121   PetscFunctionReturn(0);
1122 }
1123 
1124 PetscErrorCode DMPlexView_GLVis(DM dm, PetscViewer viewer)
1125 {
1126   PetscErrorCode ierr;
1127   PetscFunctionBegin;
1128   ierr = DMView_GLVis(dm,viewer,DMPlexView_GLVis_ASCII);CHKERRQ(ierr);
1129   PetscFunctionReturn(0);
1130 }
1131