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