Actual source code: plexglvis.c

  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>

  9: typedef struct {
 10:   PetscInt    nf;
 11:   VecScatter *scctx;
 12: } GLVisViewerCtx;

 14: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
 15: {
 16:   GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
 17:   PetscInt        i;

 19:   PetscFunctionBegin;
 20:   for (i = 0; i < ctx->nf; i++) PetscCall(VecScatterDestroy(&ctx->scctx[i]));
 21:   PetscCall(PetscFree(ctx->scctx));
 22:   PetscCall(PetscFree(vctx));
 23:   PetscFunctionReturn(PETSC_SUCCESS);
 24: }

 26: static PetscErrorCode DMPlexSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
 27: {
 28:   GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
 29:   PetscInt        f;

 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: }

 39: /* for FEM, it works for H1 fields only and extracts dofs at cell vertices, discarding any other dof */
 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;

 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_Internal(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;

 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++;

 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;

 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;

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];

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;

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;

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;

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;

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;

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));

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:   }

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: }

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;

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: };

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: };

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;

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: }

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;

281:   PetscFunctionBegin;
282:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
283:   PetscCall(DMGetDimension(dm, &dim));
284:   sdim = dim;
285:   if (csec) {
286:     PetscInt sStart, sEnd;

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: }

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;

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));

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;

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:   }

438:   /* Cleanup */
439:   PetscCall(PetscSpaceDestroy(&P));
440:   PetscCall(PetscDualSpaceDestroy(&Q));
441:   PetscCall(PetscQuadratureDestroy(&q));
442:   PetscCall(PetscQuadratureDestroy(&fq));
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: /*
447:    ASCII visualization/dump: full support for simplices and tensor product cells. It supports AMR
448:    Higher order meshes are also supported
449: */
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] = "";

468:   PetscFunctionBegin;
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));

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, (void **)&glvis_info));
484:     enabled = glvis_info->enabled;
485:     fmt     = glvis_info->fmt;
486:   }

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;

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;

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));

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:   }

527:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
528:   PetscCall(DMPlexGetGhostCellStratum(dm, &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");

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;

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:   }

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;

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_Internal(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:   }

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:   }

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;

642:       PetscCall(PetscSNPrintf(fec, sizeof(fec), "FiniteElementCollection: L2_T1_%" PetscInt_FMT "D_P1", dim));
643:       if (cEnd - cStart) {
644:         PetscInt fpc;

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;

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;

715:   /* header */
716:   PetscCall(PetscViewerASCIIPrintf(viewer, "MFEM mesh %s\n", enable_ncmesh ? "v1.1" : "v1.0"));

718:   /* topological dimension */
719:   PetscCall(PetscViewerASCIIPrintf(viewer, "\ndimension\n"));
720:   PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", dim));

722:   /* elements */
723:   minl  = 1;
724:   label = NULL;
725:   if (enable_emark) {
726:     PetscInt lminl = PETSC_MAX_INT;

728:     PetscCall(DMGetLabel(dm, emark, &label));
729:     if (label) {
730:       IS       vals;
731:       PetscInt ldef;

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:     PetscCall(MPIU_Allreduce(&lminl, &minl, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm)));
740:     if (minl == PETSC_MAX_INT) 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;

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:   }

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;

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;

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;

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;

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;

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;

820:           PetscCall(DMPlexGetCone(dm, p, &cone));
821:           PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
822:           for (i = 0; i < coneSize; i++) {
823:             PetscInt v;

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;

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;

843:             PetscCall(DMPlexGetCone(dm, p, &cone));
844:             PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
845:             for (i = 0; i < coneSize; i++) {
846:               PetscInt v;

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;

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;

868:       PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
869:       if (pown) {
870:         PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE;
871:         PetscInt  i;

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;

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;

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;

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;

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_MAX_INT;

929:       PetscCall(DMGetLabel(dm, bmark, &label));
930:       if (label) {
931:         IS       vals;
932:         PetscInt ldef;

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:       PetscCall(MPIU_Allreduce(&lminl, &minl, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm)));
941:       if (minl == PETSC_MAX_INT) 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;

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;

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;

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:   }

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;

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];

1004:         if (pp >= vStart && pp < vEnd) PetscCall(PetscBTSet(vown, pp - vStart));
1005:       }
1006:       PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1007:     }
1008:   }

1010:   if (parentSection) {
1011:     PetscInt vp, gvp;

1013:     for (vp = 0, p = vStart; p < vEnd; p++) {
1014:       DMLabel  dlabel;
1015:       PetscInt parent, depth;

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:     PetscCall(MPIU_Allreduce(&vp, &gvp, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
1024:     if (gvp) {
1025:       PetscInt   maxsupp;
1026:       PetscBool *skip = NULL;

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;

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;

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;

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;

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));

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;

1134:       PetscCall(VecGetDM(hovec, &cdm));
1135:       for (p = cStart; p < cEnd; p++) {
1136:         PetscScalar *vals = NULL;
1137:         PetscInt     csize;

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]);

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: }

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: }