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:   for (i = 0; i < ctx->nf; i++) VecScatterDestroy(&ctx->scctx[i]);
 20:   PetscFree(ctx->scctx);
 21:   PetscFree(vctx);
 22:   return 0;
 23: }

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

 30:   for (f = 0; f < nf; f++) {
 31:     VecScatterBegin(ctx->scctx[f], (Vec)oX, (Vec)oXfield[f], INSERT_VALUES, SCATTER_FORWARD);
 32:     VecScatterEnd(ctx->scctx[f], (Vec)oX, (Vec)oXfield[f], INSERT_VALUES, SCATTER_FORWARD);
 33:   }
 34:   return 0;
 35: }

 37: /* for FEM, it works for H1 fields only and extracts dofs at cell vertices, discarding any other dof */
 38: PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject odm, PetscViewer viewer)
 39: {
 40:   DM              dm = (DM)odm;
 41:   Vec             xlocal, xfield, *Ufield;
 42:   PetscDS         ds;
 43:   IS              globalNum, isfield;
 44:   PetscBT         vown;
 45:   char          **fieldname = NULL, **fec_type = NULL;
 46:   const PetscInt *gNum;
 47:   PetscInt       *nlocal, *bs, *idxs, *dims;
 48:   PetscInt        f, maxfields, nfields, c, totc, totdofs, Nv, cum, i;
 49:   PetscInt        dim, cStart, cEnd, vStart, vEnd;
 50:   GLVisViewerCtx *ctx;
 51:   PetscSection    s;

 53:   DMGetDimension(dm, &dim);
 54:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
 55:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
 56:   PetscObjectQuery((PetscObject)dm, "_glvis_plex_gnum", (PetscObject *)&globalNum);
 57:   if (!globalNum) {
 58:     DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, &globalNum);
 59:     PetscObjectCompose((PetscObject)dm, "_glvis_plex_gnum", (PetscObject)globalNum);
 60:     PetscObjectDereference((PetscObject)globalNum);
 61:   }
 62:   ISGetIndices(globalNum, &gNum);
 63:   PetscBTCreate(vEnd - vStart, &vown);
 64:   for (c = cStart, totc = 0; c < cEnd; c++) {
 65:     if (gNum[c - cStart] >= 0) {
 66:       PetscInt i, numPoints, *points = NULL;

 68:       totc++;
 69:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &numPoints, &points);
 70:       for (i = 0; i < numPoints * 2; i += 2) {
 71:         if ((points[i] >= vStart) && (points[i] < vEnd)) PetscBTSet(vown, points[i] - vStart);
 72:       }
 73:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &numPoints, &points);
 74:     }
 75:   }
 76:   for (f = 0, Nv = 0; f < vEnd - vStart; f++)
 77:     if (PetscLikely(PetscBTLookup(vown, f))) Nv++;

 79:   DMCreateLocalVector(dm, &xlocal);
 80:   VecGetLocalSize(xlocal, &totdofs);
 81:   DMGetLocalSection(dm, &s);
 82:   PetscSectionGetNumFields(s, &nfields);
 83:   for (f = 0, maxfields = 0; f < nfields; f++) {
 84:     PetscInt bs;

 86:     PetscSectionGetFieldComponents(s, f, &bs);
 87:     maxfields += bs;
 88:   }
 89:   PetscCalloc7(maxfields, &fieldname, maxfields, &nlocal, maxfields, &bs, maxfields, &dims, maxfields, &fec_type, totdofs, &idxs, maxfields, &Ufield);
 90:   PetscNew(&ctx);
 91:   PetscCalloc1(maxfields, &ctx->scctx);
 92:   DMGetDS(dm, &ds);
 93:   if (ds) {
 94:     for (f = 0; f < nfields; f++) {
 95:       const char *fname;
 96:       char        name[256];
 97:       PetscObject disc;
 98:       size_t      len;

100:       PetscSectionGetFieldName(s, f, &fname);
101:       PetscStrlen(fname, &len);
102:       if (len) {
103:         PetscStrcpy(name, fname);
104:       } else {
105:         PetscSNPrintf(name, 256, "Field%" PetscInt_FMT, f);
106:       }
107:       PetscDSGetDiscretization(ds, f, &disc);
108:       if (disc) {
109:         PetscClassId id;
110:         PetscInt     Nc;
111:         char         fec[64];

113:         PetscObjectGetClassId(disc, &id);
114:         if (id == PETSCFE_CLASSID) {
115:           PetscFE            fem = (PetscFE)disc;
116:           PetscDualSpace     sp;
117:           PetscDualSpaceType spname;
118:           PetscInt           order;
119:           PetscBool          islag, continuous, H1 = PETSC_TRUE;

121:           PetscFEGetNumComponents(fem, &Nc);
122:           PetscFEGetDualSpace(fem, &sp);
123:           PetscDualSpaceGetType(sp, &spname);
124:           PetscStrcmp(spname, PETSCDUALSPACELAGRANGE, &islag);
126:           PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
127:           PetscDualSpaceGetOrder(sp, &order);
128:           if (continuous && order > 0) { /* no support for high-order viz, still have to figure out the numbering */
129:             PetscSNPrintf(fec, 64, "FiniteElementCollection: H1_%" PetscInt_FMT "D_P1", dim);
130:           } else {
132:             H1 = PETSC_FALSE;
133:             PetscSNPrintf(fec, 64, "FiniteElementCollection: L2_%" PetscInt_FMT "D_P%" PetscInt_FMT, dim, order);
134:           }
135:           PetscStrallocpy(name, &fieldname[ctx->nf]);
136:           bs[ctx->nf]   = Nc;
137:           dims[ctx->nf] = dim;
138:           if (H1) {
139:             nlocal[ctx->nf] = Nc * Nv;
140:             PetscStrallocpy(fec, &fec_type[ctx->nf]);
141:             VecCreateSeq(PETSC_COMM_SELF, Nv * Nc, &xfield);
142:             for (i = 0, cum = 0; i < vEnd - vStart; i++) {
143:               PetscInt j, off;

145:               if (PetscUnlikely(!PetscBTLookup(vown, i))) continue;
146:               PetscSectionGetFieldOffset(s, i + vStart, f, &off);
147:               for (j = 0; j < Nc; j++) idxs[cum++] = off + j;
148:             }
149:             ISCreateGeneral(PetscObjectComm((PetscObject)xlocal), Nv * Nc, idxs, PETSC_USE_POINTER, &isfield);
150:           } else {
151:             nlocal[ctx->nf] = Nc * totc;
152:             PetscStrallocpy(fec, &fec_type[ctx->nf]);
153:             VecCreateSeq(PETSC_COMM_SELF, Nc * totc, &xfield);
154:             for (i = 0, cum = 0; i < cEnd - cStart; i++) {
155:               PetscInt j, off;

157:               if (PetscUnlikely(gNum[i] < 0)) continue;
158:               PetscSectionGetFieldOffset(s, i + cStart, f, &off);
159:               for (j = 0; j < Nc; j++) idxs[cum++] = off + j;
160:             }
161:             ISCreateGeneral(PetscObjectComm((PetscObject)xlocal), totc * Nc, idxs, PETSC_USE_POINTER, &isfield);
162:           }
163:           VecScatterCreate(xlocal, isfield, xfield, NULL, &ctx->scctx[ctx->nf]);
164:           VecDestroy(&xfield);
165:           ISDestroy(&isfield);
166:           ctx->nf++;
167:         } else if (id == PETSCFV_CLASSID) {
168:           PetscInt c;

170:           PetscFVGetNumComponents((PetscFV)disc, &Nc);
171:           PetscSNPrintf(fec, 64, "FiniteElementCollection: L2_%" PetscInt_FMT "D_P0", dim);
172:           for (c = 0; c < Nc; c++) {
173:             char comp[256];
174:             PetscSNPrintf(comp, 256, "%s-Comp%" PetscInt_FMT, name, c);
175:             PetscStrallocpy(comp, &fieldname[ctx->nf]);
176:             bs[ctx->nf]     = 1; /* Does PetscFV support components with different block size? */
177:             nlocal[ctx->nf] = totc;
178:             dims[ctx->nf]   = dim;
179:             PetscStrallocpy(fec, &fec_type[ctx->nf]);
180:             VecCreateSeq(PETSC_COMM_SELF, totc, &xfield);
181:             for (i = 0, cum = 0; i < cEnd - cStart; i++) {
182:               PetscInt off;

184:               if (PetscUnlikely(gNum[i]) < 0) continue;
185:               PetscSectionGetFieldOffset(s, i + cStart, f, &off);
186:               idxs[cum++] = off + c;
187:             }
188:             ISCreateGeneral(PetscObjectComm((PetscObject)xlocal), totc, idxs, PETSC_USE_POINTER, &isfield);
189:             VecScatterCreate(xlocal, isfield, xfield, NULL, &ctx->scctx[ctx->nf]);
190:             VecDestroy(&xfield);
191:             ISDestroy(&isfield);
192:             ctx->nf++;
193:           }
194:         } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
195:       } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing discretization for field %" PetscInt_FMT, f);
196:     }
197:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Needs a DS attached to the DM");
198:   PetscBTDestroy(&vown);
199:   VecDestroy(&xlocal);
200:   ISRestoreIndices(globalNum, &gNum);

202:   /* create work vectors */
203:   for (f = 0; f < ctx->nf; f++) {
204:     VecCreateMPI(PetscObjectComm((PetscObject)dm), nlocal[f], PETSC_DECIDE, &Ufield[f]);
205:     PetscObjectSetName((PetscObject)Ufield[f], fieldname[f]);
206:     VecSetBlockSize(Ufield[f], bs[f]);
207:     VecSetDM(Ufield[f], dm);
208:   }

210:   /* customize the viewer */
211:   PetscViewerGLVisSetFields(viewer, ctx->nf, (const char **)fec_type, dims, DMPlexSampleGLVisFields_Private, (PetscObject *)Ufield, ctx, DestroyGLVisViewerCtx_Private);
212:   for (f = 0; f < ctx->nf; f++) {
213:     PetscFree(fieldname[f]);
214:     PetscFree(fec_type[f]);
215:     VecDestroy(&Ufield[f]);
216:   }
217:   PetscFree7(fieldname, nlocal, bs, dims, fec_type, idxs, Ufield);
218:   return 0;
219: }

221: typedef enum {
222:   MFEM_POINT = 0,
223:   MFEM_SEGMENT,
224:   MFEM_TRIANGLE,
225:   MFEM_SQUARE,
226:   MFEM_TETRAHEDRON,
227:   MFEM_CUBE,
228:   MFEM_PRISM,
229:   MFEM_UNDEF
230: } MFEM_cid;

232: MFEM_cid mfem_table_cid[4][7] = {
233:   {MFEM_POINT, MFEM_UNDEF, MFEM_UNDEF,   MFEM_UNDEF,    MFEM_UNDEF,       MFEM_UNDEF, MFEM_UNDEF},
234:   {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_UNDEF,    MFEM_UNDEF,       MFEM_UNDEF, MFEM_UNDEF},
235:   {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_TRIANGLE, MFEM_SQUARE,      MFEM_UNDEF, MFEM_UNDEF},
236:   {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_UNDEF,    MFEM_TETRAHEDRON, MFEM_PRISM, MFEM_CUBE }
237: };

239: MFEM_cid mfem_table_cid_unint[4][9] = {
240:   {MFEM_POINT, MFEM_UNDEF, MFEM_UNDEF,   MFEM_UNDEF,    MFEM_UNDEF,       MFEM_UNDEF, MFEM_PRISM, MFEM_UNDEF, MFEM_UNDEF},
241:   {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_UNDEF,    MFEM_UNDEF,       MFEM_UNDEF, MFEM_PRISM, MFEM_UNDEF, MFEM_UNDEF},
242:   {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_TRIANGLE, MFEM_SQUARE,      MFEM_UNDEF, MFEM_PRISM, MFEM_UNDEF, MFEM_UNDEF},
243:   {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_UNDEF,    MFEM_TETRAHEDRON, MFEM_UNDEF, MFEM_PRISM, MFEM_UNDEF, MFEM_CUBE }
244: };

246: static PetscErrorCode DMPlexGetPointMFEMCellID_Internal(DM dm, DMLabel label, PetscInt minl, PetscInt p, PetscInt *mid, PetscInt *cid)
247: {
248:   DMLabel  dlabel;
249:   PetscInt depth, csize, pdepth, dim;

251:   DMPlexGetDepthLabel(dm, &dlabel);
252:   DMLabelGetValue(dlabel, p, &pdepth);
253:   DMPlexGetConeSize(dm, p, &csize);
254:   DMPlexGetDepth(dm, &depth);
255:   DMGetDimension(dm, &dim);
256:   if (label) {
257:     DMLabelGetValue(label, p, mid);
258:     *mid = *mid - minl + 1; /* MFEM does not like negative markers */
259:   } else *mid = 1;
260:   if (depth >= 0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
264:     *cid = mfem_table_cid_unint[dim][csize];
265:   } else {
268:     *cid = mfem_table_cid[pdepth][csize];
269:   }
270:   return 0;
271: }

273: static PetscErrorCode DMPlexGetPointMFEMVertexIDs_Internal(DM dm, PetscInt p, PetscSection csec, PetscInt *nv, PetscInt vids[])
274: {
275:   PetscInt dim, sdim, dof = 0, off = 0, i, q, vStart, vEnd, numPoints, *points = NULL;

277:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
278:   DMGetDimension(dm, &dim);
279:   sdim = dim;
280:   if (csec) {
281:     PetscInt sStart, sEnd;

283:     DMGetCoordinateDim(dm, &sdim);
284:     PetscSectionGetChart(csec, &sStart, &sEnd);
285:     PetscSectionGetOffset(csec, vStart, &off);
286:     off = off / sdim;
287:     if (p >= sStart && p < sEnd) PetscSectionGetDof(csec, p, &dof);
288:   }
289:   if (!dof) {
290:     DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &numPoints, &points);
291:     for (i = 0, q = 0; i < numPoints * 2; i += 2)
292:       if ((points[i] >= vStart) && (points[i] < vEnd)) vids[q++] = points[i] - vStart + off;
293:     DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &numPoints, &points);
294:   } else {
295:     PetscSectionGetOffset(csec, p, &off);
296:     PetscSectionGetDof(csec, p, &dof);
297:     for (q = 0; q < dof / sdim; q++) vids[q] = off / sdim + q;
298:   }
299:   *nv = q;
300:   return 0;
301: }

303: static PetscErrorCode GLVisCreateFE(PetscFE femIn, char name[32], PetscFE *fem, IS *perm)
304: {
305:   DM              K;
306:   PetscSpace      P;
307:   PetscDualSpace  Q;
308:   PetscQuadrature q, fq;
309:   PetscInt        dim, deg, dof;
310:   DMPolytopeType  ptype;
311:   PetscBool       isSimplex, isTensor;
312:   PetscBool       continuity = PETSC_FALSE;
313:   PetscDTNodeType nodeType   = PETSCDTNODES_GAUSSJACOBI;
314:   PetscBool       endpoint   = PETSC_TRUE;
315:   MPI_Comm        comm;

317:   PetscObjectGetComm((PetscObject)femIn, &comm);
318:   PetscFEGetBasisSpace(femIn, &P);
319:   PetscFEGetDualSpace(femIn, &Q);
320:   PetscDualSpaceGetDM(Q, &K);
321:   DMGetDimension(K, &dim);
322:   PetscSpaceGetDegree(P, &deg, NULL);
323:   PetscSpaceGetNumComponents(P, &dof);
324:   DMPlexGetCellType(K, 0, &ptype);
325:   switch (ptype) {
326:   case DM_POLYTOPE_QUADRILATERAL:
327:   case DM_POLYTOPE_HEXAHEDRON:
328:     isSimplex = PETSC_FALSE;
329:     break;
330:   default:
331:     isSimplex = PETSC_TRUE;
332:     break;
333:   }
334:   isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
335:   if (isSimplex) deg = PetscMin(deg, 3); /* Permutation not coded for degree higher than 3 */
336:   /* Create space */
337:   PetscSpaceCreate(comm, &P);
338:   PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
339:   PetscSpacePolynomialSetTensor(P, isTensor);
340:   PetscSpaceSetNumComponents(P, dof);
341:   PetscSpaceSetNumVariables(P, dim);
342:   PetscSpaceSetDegree(P, deg, PETSC_DETERMINE);
343:   PetscSpaceSetUp(P);
344:   /* Create dual space */
345:   PetscDualSpaceCreate(comm, &Q);
346:   PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
347:   PetscDualSpaceLagrangeSetTensor(Q, isTensor);
348:   PetscDualSpaceLagrangeSetContinuity(Q, continuity);
349:   PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0);
350:   PetscDualSpaceSetNumComponents(Q, dof);
351:   PetscDualSpaceSetOrder(Q, deg);
352:   DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K);
353:   PetscDualSpaceSetDM(Q, K);
354:   DMDestroy(&K);
355:   PetscDualSpaceSetUp(Q);
356:   /* Create quadrature */
357:   if (isSimplex) {
358:     PetscDTStroudConicalQuadrature(dim, 1, deg + 1, -1, +1, &q);
359:     PetscDTStroudConicalQuadrature(dim - 1, 1, deg + 1, -1, +1, &fq);
360:   } else {
361:     PetscDTGaussTensorQuadrature(dim, 1, deg + 1, -1, +1, &q);
362:     PetscDTGaussTensorQuadrature(dim - 1, 1, deg + 1, -1, +1, &fq);
363:   }
364:   /* Create finite element */
365:   PetscFECreate(comm, fem);
366:   PetscSNPrintf(name, 32, "L2_T1_%" PetscInt_FMT "D_P%" PetscInt_FMT, dim, deg);
367:   PetscObjectSetName((PetscObject)*fem, name);
368:   PetscFESetType(*fem, PETSCFEBASIC);
369:   PetscFESetNumComponents(*fem, dof);
370:   PetscFESetBasisSpace(*fem, P);
371:   PetscFESetDualSpace(*fem, Q);
372:   PetscFESetQuadrature(*fem, q);
373:   PetscFESetFaceQuadrature(*fem, fq);
374:   PetscFESetUp(*fem);

376:   /* Both MFEM and PETSc are lexicographic, but PLEX stores the swapped cone */
377:   *perm = NULL;
378:   if (isSimplex && dim == 3) {
379:     PetscInt celldofs, *pidx;

381:     PetscDualSpaceGetDimension(Q, &celldofs);
382:     celldofs /= dof;
383:     PetscMalloc1(celldofs, &pidx);
384:     switch (celldofs) {
385:     case 4:
386:       pidx[0] = 2;
387:       pidx[1] = 0;
388:       pidx[2] = 1;
389:       pidx[3] = 3;
390:       break;
391:     case 10:
392:       pidx[0] = 5;
393:       pidx[1] = 3;
394:       pidx[2] = 0;
395:       pidx[3] = 4;
396:       pidx[4] = 1;
397:       pidx[5] = 2;
398:       pidx[6] = 8;
399:       pidx[7] = 6;
400:       pidx[8] = 7;
401:       pidx[9] = 9;
402:       break;
403:     case 20:
404:       pidx[0]  = 9;
405:       pidx[1]  = 7;
406:       pidx[2]  = 4;
407:       pidx[3]  = 0;
408:       pidx[4]  = 8;
409:       pidx[5]  = 5;
410:       pidx[6]  = 1;
411:       pidx[7]  = 6;
412:       pidx[8]  = 2;
413:       pidx[9]  = 3;
414:       pidx[10] = 15;
415:       pidx[11] = 13;
416:       pidx[12] = 10;
417:       pidx[13] = 14;
418:       pidx[14] = 11;
419:       pidx[15] = 12;
420:       pidx[16] = 18;
421:       pidx[17] = 16;
422:       pidx[18] = 17;
423:       pidx[19] = 19;
424:       break;
425:     default:
426:       SETERRQ(comm, PETSC_ERR_SUP, "Unhandled degree,dof pair %" PetscInt_FMT ",%" PetscInt_FMT, deg, celldofs);
427:       break;
428:     }
429:     ISCreateBlock(PETSC_COMM_SELF, dof, celldofs, pidx, PETSC_OWN_POINTER, perm);
430:   }

432:   /* Cleanup */
433:   PetscSpaceDestroy(&P);
434:   PetscDualSpaceDestroy(&Q);
435:   PetscQuadratureDestroy(&q);
436:   PetscQuadratureDestroy(&fq);
437:   return 0;
438: }

440: /*
441:    ASCII visualization/dump: full support for simplices and tensor product cells. It supports AMR
442:    Higher order meshes are also supported
443: */
444: static PetscErrorCode DMPlexView_GLVis_ASCII(DM dm, PetscViewer viewer)
445: {
446:   DMLabel            label;
447:   PetscSection       coordSection, coordSectionCell, parentSection, hoSection = NULL;
448:   Vec                coordinates, coordinatesCell, hovec;
449:   const PetscScalar *array;
450:   PetscInt           bf, p, sdim, dim, depth, novl, minl;
451:   PetscInt           cStart, cEnd, vStart, vEnd, nvert;
452:   PetscMPIInt        size;
453:   PetscBool          localized, isascii;
454:   PetscBool          enable_mfem, enable_boundary, enable_ncmesh, view_ovl = PETSC_FALSE;
455:   PetscBT            pown, vown;
456:   PetscContainer     glvis_container;
457:   PetscBool          cellvertex = PETSC_FALSE, enabled = PETSC_TRUE;
458:   PetscBool          enable_emark, enable_bmark;
459:   const char        *fmt;
460:   char               emark[64] = "", bmark[64] = "";

464:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
466:   MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size);
468:   DMGetDimension(dm, &dim);
469:   DMPlexGetDepth(dm, &depth);

471:   /* get container: determines if a process visualizes is portion of the data or not */
472:   PetscObjectQuery((PetscObject)viewer, "_glvis_info_container", (PetscObject *)&glvis_container);
474:   {
475:     PetscViewerGLVisInfo glvis_info;
476:     PetscContainerGetPointer(glvis_container, (void **)&glvis_info);
477:     enabled = glvis_info->enabled;
478:     fmt     = glvis_info->fmt;
479:   }

481:   /* Users can attach a coordinate vector to the DM in case they have a higher-order mesh */
482:   PetscObjectQuery((PetscObject)dm, "_glvis_mesh_coords", (PetscObject *)&hovec);
483:   PetscObjectReference((PetscObject)hovec);
484:   if (!hovec) {
485:     DM           cdm;
486:     PetscFE      disc;
487:     PetscClassId classid;

489:     DMGetCoordinateDM(dm, &cdm);
490:     DMGetField(cdm, 0, NULL, (PetscObject *)&disc);
491:     PetscObjectGetClassId((PetscObject)disc, &classid);
492:     if (classid == PETSCFE_CLASSID) {
493:       DM      hocdm;
494:       PetscFE hodisc;
495:       Vec     vec;
496:       Mat     mat;
497:       char    name[32], fec_type[64];
498:       IS      perm = NULL;

500:       GLVisCreateFE(disc, name, &hodisc, &perm);
501:       DMClone(cdm, &hocdm);
502:       DMSetField(hocdm, 0, NULL, (PetscObject)hodisc);
503:       PetscFEDestroy(&hodisc);
504:       DMCreateDS(hocdm);

506:       DMGetCoordinates(dm, &vec);
507:       DMCreateGlobalVector(hocdm, &hovec);
508:       DMCreateInterpolation(cdm, hocdm, &mat, NULL);
509:       MatInterpolate(mat, vec, hovec);
510:       MatDestroy(&mat);
511:       DMGetLocalSection(hocdm, &hoSection);
512:       PetscSectionSetClosurePermutation(hoSection, (PetscObject)hocdm, depth, perm);
513:       ISDestroy(&perm);
514:       DMDestroy(&hocdm);
515:       PetscSNPrintf(fec_type, sizeof(fec_type), "FiniteElementCollection: %s", name);
516:       PetscObjectSetName((PetscObject)hovec, fec_type);
517:     }
518:   }

520:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
521:   DMPlexGetGhostCellStratum(dm, &p, NULL);
522:   if (p >= 0) cEnd = p;
523:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
524:   DMGetCoordinatesLocalized(dm, &localized);
525:   DMGetCoordinateSection(dm, &coordSection);
526:   DMGetCoordinateDim(dm, &sdim);
527:   DMGetCoordinatesLocal(dm, &coordinates);

530:   /*
531:      a couple of sections of the mesh specification are disabled
532:        - 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)
533:        - vertex_parents: used for non-conforming meshes only when we want to use MFEM as a discretization package
534:                          and be able to derefine the mesh (MFEM does not currently have to ability to read ncmeshes in parallel)
535:   */
536:   enable_boundary = PETSC_FALSE;
537:   enable_ncmesh   = PETSC_FALSE;
538:   enable_mfem     = PETSC_FALSE;
539:   enable_emark    = PETSC_FALSE;
540:   enable_bmark    = PETSC_FALSE;
541:   /* I'm tired of problems with negative values in the markers, disable them */
542:   PetscOptionsBegin(PetscObjectComm((PetscObject)dm), ((PetscObject)dm)->prefix, "GLVis PetscViewer DMPlex Options", "PetscViewer");
543:   PetscOptionsBool("-viewer_glvis_dm_plex_enable_boundary", "Enable boundary section in mesh representation", NULL, enable_boundary, &enable_boundary, NULL);
544:   PetscOptionsBool("-viewer_glvis_dm_plex_enable_ncmesh", "Enable vertex_parents section in mesh representation (allows derefinement)", NULL, enable_ncmesh, &enable_ncmesh, NULL);
545:   PetscOptionsBool("-viewer_glvis_dm_plex_enable_mfem", "Dump a mesh that can be used with MFEM's FiniteElementSpaces", NULL, enable_mfem, &enable_mfem, NULL);
546:   PetscOptionsBool("-viewer_glvis_dm_plex_overlap", "Include overlap region in local meshes", NULL, view_ovl, &view_ovl, NULL);
547:   PetscOptionsString("-viewer_glvis_dm_plex_emarker", "String for the material id label", NULL, emark, emark, sizeof(emark), &enable_emark);
548:   PetscOptionsString("-viewer_glvis_dm_plex_bmarker", "String for the boundary id label", NULL, bmark, bmark, sizeof(bmark), &enable_bmark);
549:   PetscOptionsEnd();
550:   if (enable_bmark) enable_boundary = PETSC_TRUE;

552:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
555:              "Mesh must be interpolated. "
556:              "Alternatively, run with -viewer_glvis_dm_plex_enable_boundary 0");
558:              "Mesh must be interpolated. "
559:              "Alternatively, run with -viewer_glvis_dm_plex_enable_ncmesh 0");
560:   if (depth >= 0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
562:     cellvertex = PETSC_TRUE;
563:   }

565:   /* Identify possible cells in the overlap */
566:   novl = 0;
567:   pown = NULL;
568:   if (size > 1) {
569:     IS              globalNum = NULL;
570:     const PetscInt *gNum;
571:     PetscBool       ovl = PETSC_FALSE;

573:     PetscObjectQuery((PetscObject)dm, "_glvis_plex_gnum", (PetscObject *)&globalNum);
574:     if (!globalNum) {
575:       if (view_ovl) {
576:         ISCreateStride(PetscObjectComm((PetscObject)dm), cEnd - cStart, 0, 1, &globalNum);
577:       } else {
578:         DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, &globalNum);
579:       }
580:       PetscObjectCompose((PetscObject)dm, "_glvis_plex_gnum", (PetscObject)globalNum);
581:       PetscObjectDereference((PetscObject)globalNum);
582:     }
583:     ISGetIndices(globalNum, &gNum);
584:     for (p = cStart; p < cEnd; p++) {
585:       if (gNum[p - cStart] < 0) {
586:         ovl = PETSC_TRUE;
587:         novl++;
588:       }
589:     }
590:     if (ovl) {
591:       /* it may happen that pown get not destroyed, if the user closes the window while this function is running.
592:          TODO: garbage collector? attach pown to dm?  */
593:       PetscBTCreate(cEnd - cStart, &pown);
594:       for (p = cStart; p < cEnd; p++) {
595:         if (gNum[p - cStart] < 0) continue;
596:         else PetscBTSet(pown, p - cStart);
597:       }
598:     }
599:     ISRestoreIndices(globalNum, &gNum);
600:   }

602:   /* vertex_parents (Non-conforming meshes) */
603:   parentSection = NULL;
604:   if (enable_ncmesh) {
605:     DMPlexGetTree(dm, &parentSection, NULL, NULL, NULL, NULL);
606:     enable_ncmesh = (PetscBool)(enable_ncmesh && parentSection);
607:   }
608:   /* return if this process is disabled */
609:   if (!enabled) {
610:     PetscViewerASCIIPrintf(viewer, "MFEM mesh %s\n", enable_ncmesh ? "v1.1" : "v1.0");
611:     PetscViewerASCIIPrintf(viewer, "\ndimension\n");
612:     PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", dim);
613:     PetscViewerASCIIPrintf(viewer, "\nelements\n");
614:     PetscViewerASCIIPrintf(viewer, "0\n");
615:     PetscViewerASCIIPrintf(viewer, "\nboundary\n");
616:     PetscViewerASCIIPrintf(viewer, "0\n");
617:     PetscViewerASCIIPrintf(viewer, "\nvertices\n");
618:     PetscViewerASCIIPrintf(viewer, "0\n");
619:     PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", sdim);
620:     PetscBTDestroy(&pown);
621:     VecDestroy(&hovec);
622:     return 0;
623:   }

625:   if (enable_mfem) {
626:     if (localized && !hovec) { /* we need to generate a vector of L2 coordinates, as this is how MFEM handles periodic meshes */
627:       PetscInt     vpc = 0;
628:       char         fec[64];
629:       PetscInt     vids[8] = {0, 1, 2, 3, 4, 5, 6, 7};
630:       PetscInt     hexv[8] = {0, 1, 3, 2, 4, 5, 7, 6}, tetv[4] = {0, 1, 2, 3};
631:       PetscInt     quadv[8] = {0, 1, 3, 2}, triv[3] = {0, 1, 2};
632:       PetscInt    *dof = NULL;
633:       PetscScalar *array, *ptr;

635:       PetscSNPrintf(fec, sizeof(fec), "FiniteElementCollection: L2_T1_%" PetscInt_FMT "D_P1", dim);
636:       if (cEnd - cStart) {
637:         PetscInt fpc;

639:         DMPlexGetConeSize(dm, cStart, &fpc);
640:         switch (dim) {
641:         case 1:
642:           vpc = 2;
643:           dof = hexv;
644:           break;
645:         case 2:
646:           switch (fpc) {
647:           case 3:
648:             vpc = 3;
649:             dof = triv;
650:             break;
651:           case 4:
652:             vpc = 4;
653:             dof = quadv;
654:             break;
655:           default:
656:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unhandled case: faces per cell %" PetscInt_FMT, fpc);
657:           }
658:           break;
659:         case 3:
660:           switch (fpc) {
661:           case 4: /* TODO: still need to understand L2 ordering for tets */
662:             vpc = 4;
663:             dof = tetv;
664:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unhandled tethraedral case");
665:           case 6:
667:             vpc = 8;
668:             dof = hexv;
669:             break;
670:           case 8:
672:             vpc = 8;
673:             dof = hexv;
674:             break;
675:           default:
676:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unhandled case: faces per cell %" PetscInt_FMT, fpc);
677:           }
678:           break;
679:         default:
680:           SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unhandled dim");
681:         }
682:         DMPlexReorderCell(dm, cStart, vids);
683:       }
685:       VecCreateSeq(PETSC_COMM_SELF, (cEnd - cStart - novl) * vpc * sdim, &hovec);
686:       PetscObjectSetName((PetscObject)hovec, fec);
687:       VecGetArray(hovec, &array);
688:       ptr = array;
689:       for (p = cStart; p < cEnd; p++) {
690:         PetscInt     csize, v, d;
691:         PetscScalar *vals = NULL;

693:         if (PetscUnlikely(pown && !PetscBTLookup(pown, p - cStart))) continue;
694:         DMPlexVecGetClosure(dm, coordSection, coordinates, p, &csize, &vals);
696:         for (v = 0; v < vpc; v++) {
697:           for (d = 0; d < sdim; d++) ptr[sdim * dof[v] + d] = vals[sdim * vids[v] + d];
698:         }
699:         ptr += vpc * sdim;
700:         DMPlexVecRestoreClosure(dm, coordSection, coordinates, p, &csize, &vals);
701:       }
702:       VecRestoreArray(hovec, &array);
703:     }
704:   }
705:   /* if we have high-order coordinates in 3D, we need to specify the boundary */
706:   if (hovec && dim == 3) enable_boundary = PETSC_TRUE;

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

711:   /* topological dimension */
712:   PetscViewerASCIIPrintf(viewer, "\ndimension\n");
713:   PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", dim);

715:   /* elements */
716:   minl  = 1;
717:   label = NULL;
718:   if (enable_emark) {
719:     PetscInt lminl = PETSC_MAX_INT;

721:     DMGetLabel(dm, emark, &label);
722:     if (label) {
723:       IS       vals;
724:       PetscInt ldef;

726:       DMLabelGetDefaultValue(label, &ldef);
727:       DMLabelGetValueIS(label, &vals);
728:       ISGetMinMax(vals, &lminl, NULL);
729:       ISDestroy(&vals);
730:       lminl = PetscMin(ldef, lminl);
731:     }
732:     MPIU_Allreduce(&lminl, &minl, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
733:     if (minl == PETSC_MAX_INT) minl = 1;
734:   }
735:   PetscViewerASCIIPrintf(viewer, "\nelements\n");
736:   PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", cEnd - cStart - novl);
737:   for (p = cStart; p < cEnd; p++) {
738:     PetscInt vids[8];
739:     PetscInt i, nv = 0, cid = -1, mid = 1;

741:     if (PetscUnlikely(pown && !PetscBTLookup(pown, p - cStart))) continue;
742:     DMPlexGetPointMFEMCellID_Internal(dm, label, minl, p, &mid, &cid);
743:     DMPlexGetPointMFEMVertexIDs_Internal(dm, p, (localized && !hovec) ? coordSection : NULL, &nv, vids);
744:     DMPlexReorderCell(dm, p, vids);
745:     PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT, mid, cid);
746:     for (i = 0; i < nv; i++) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, vids[i]);
747:     PetscViewerASCIIPrintf(viewer, "\n");
748:   }

750:   /* boundary */
751:   PetscViewerASCIIPrintf(viewer, "\nboundary\n");
752:   if (!enable_boundary) {
753:     PetscViewerASCIIPrintf(viewer, "0\n");
754:   } else {
755:     DMLabel  perLabel;
756:     PetscBT  bfaces;
757:     PetscInt fStart, fEnd, *fcells;

759:     DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
760:     PetscBTCreate(fEnd - fStart, &bfaces);
761:     DMPlexGetMaxSizes(dm, NULL, &p);
762:     PetscMalloc1(p, &fcells);
763:     DMGetLabel(dm, "glvis_periodic_cut", &perLabel);
764:     if (!perLabel && localized) { /* this periodic cut can be moved up to DMPlex setup */
765:       DMCreateLabel(dm, "glvis_periodic_cut");
766:       DMGetLabel(dm, "glvis_periodic_cut", &perLabel);
767:       DMLabelSetDefaultValue(perLabel, 1);
768:       DMGetCellCoordinateSection(dm, &coordSectionCell);
769:       DMGetCellCoordinatesLocal(dm, &coordinatesCell);
770:       for (p = cStart; p < cEnd; p++) {
771:         DMPolytopeType cellType;
772:         PetscInt       dof;

774:         DMPlexGetCellType(dm, p, &cellType);
775:         PetscSectionGetDof(coordSectionCell, p, &dof);
776:         if (dof) {
777:           PetscInt     uvpc, v, csize, csizeCell, cellClosureSize, *cellClosure = NULL, *vidxs = NULL;
778:           PetscScalar *vals = NULL, *valsCell = NULL;

780:           uvpc = DMPolytopeTypeGetNumVertices(cellType);
782:           DMPlexVecGetClosure(dm, coordSection, coordinates, p, &csize, &vals);
783:           DMPlexVecGetClosure(dm, coordSectionCell, coordinatesCell, p, &csizeCell, &valsCell);
785:           DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &cellClosureSize, &cellClosure);
786:           for (v = 0; v < cellClosureSize; v++)
787:             if (cellClosure[2 * v] >= vStart && cellClosure[2 * v] < vEnd) {
788:               vidxs = cellClosure + 2 * v;
789:               break;
790:             }
792:           for (v = 0; v < uvpc; v++) {
793:             PetscInt s;

795:             for (s = 0; s < sdim; s++) {
796:               if (PetscAbsScalar(vals[v * sdim + s] - valsCell[v * sdim + s]) > PETSC_MACHINE_EPSILON) DMLabelSetValue(perLabel, vidxs[2 * v], 2);
797:             }
798:           }
799:           DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &cellClosureSize, &cellClosure);
800:           DMPlexVecRestoreClosure(dm, coordSection, coordinates, p, &csize, &vals);
801:           DMPlexVecRestoreClosure(dm, coordSectionCell, coordinatesCell, p, &csizeCell, &valsCell);
802:         }
803:       }
804:       if (dim > 1) {
805:         PetscInt eEnd, eStart;

807:         DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
808:         for (p = eStart; p < eEnd; p++) {
809:           const PetscInt *cone;
810:           PetscInt        coneSize, i;
811:           PetscBool       ispe = PETSC_TRUE;

813:           DMPlexGetCone(dm, p, &cone);
814:           DMPlexGetConeSize(dm, p, &coneSize);
815:           for (i = 0; i < coneSize; i++) {
816:             PetscInt v;

818:             DMLabelGetValue(perLabel, cone[i], &v);
819:             ispe = (PetscBool)(ispe && (v == 2));
820:           }
821:           if (ispe && coneSize) {
822:             PetscInt        ch, numChildren;
823:             const PetscInt *children;

825:             DMLabelSetValue(perLabel, p, 2);
826:             DMPlexGetTreeChildren(dm, p, &numChildren, &children);
827:             for (ch = 0; ch < numChildren; ch++) DMLabelSetValue(perLabel, children[ch], 2);
828:           }
829:         }
830:         if (dim > 2) {
831:           for (p = fStart; p < fEnd; p++) {
832:             const PetscInt *cone;
833:             PetscInt        coneSize, i;
834:             PetscBool       ispe = PETSC_TRUE;

836:             DMPlexGetCone(dm, p, &cone);
837:             DMPlexGetConeSize(dm, p, &coneSize);
838:             for (i = 0; i < coneSize; i++) {
839:               PetscInt v;

841:               DMLabelGetValue(perLabel, cone[i], &v);
842:               ispe = (PetscBool)(ispe && (v == 2));
843:             }
844:             if (ispe && coneSize) {
845:               PetscInt        ch, numChildren;
846:               const PetscInt *children;

848:               DMLabelSetValue(perLabel, p, 2);
849:               DMPlexGetTreeChildren(dm, p, &numChildren, &children);
850:               for (ch = 0; ch < numChildren; ch++) DMLabelSetValue(perLabel, children[ch], 2);
851:             }
852:           }
853:         }
854:       }
855:     }
856:     for (p = fStart; p < fEnd; p++) {
857:       const PetscInt *support;
858:       PetscInt        supportSize;
859:       PetscBool       isbf = PETSC_FALSE;

861:       DMPlexGetSupportSize(dm, p, &supportSize);
862:       if (pown) {
863:         PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE;
864:         PetscInt  i;

866:         DMPlexGetSupport(dm, p, &support);
867:         for (i = 0; i < supportSize; i++) {
868:           if (PetscLikely(PetscBTLookup(pown, support[i] - cStart))) has_owned = PETSC_TRUE;
869:           else has_ghost = PETSC_TRUE;
870:         }
871:         isbf = (PetscBool)((supportSize == 1 && has_owned) || (supportSize > 1 && has_owned && has_ghost));
872:       } else {
873:         isbf = (PetscBool)(supportSize == 1);
874:       }
875:       if (!isbf && perLabel) {
876:         const PetscInt *cone;
877:         PetscInt        coneSize, i;

879:         DMPlexGetCone(dm, p, &cone);
880:         DMPlexGetConeSize(dm, p, &coneSize);
881:         isbf = PETSC_TRUE;
882:         for (i = 0; i < coneSize; i++) {
883:           PetscInt v, d;

885:           DMLabelGetValue(perLabel, cone[i], &v);
886:           DMLabelGetDefaultValue(perLabel, &d);
887:           isbf = (PetscBool)(isbf && v != d);
888:         }
889:       }
890:       if (isbf) PetscBTSet(bfaces, p - fStart);
891:     }
892:     /* count boundary faces */
893:     for (p = fStart, bf = 0; p < fEnd; p++) {
894:       if (PetscUnlikely(PetscBTLookup(bfaces, p - fStart))) {
895:         const PetscInt *support;
896:         PetscInt        supportSize, c;

898:         DMPlexGetSupportSize(dm, p, &supportSize);
899:         DMPlexGetSupport(dm, p, &support);
900:         for (c = 0; c < supportSize; c++) {
901:           const PetscInt *cone;
902:           PetscInt        cell, cl, coneSize;

904:           cell = support[c];
905:           if (pown && PetscUnlikely(!PetscBTLookup(pown, cell - cStart))) continue;
906:           DMPlexGetCone(dm, cell, &cone);
907:           DMPlexGetConeSize(dm, cell, &coneSize);
908:           for (cl = 0; cl < coneSize; cl++) {
909:             if (cone[cl] == p) {
910:               bf += 1;
911:               break;
912:             }
913:           }
914:         }
915:       }
916:     }
917:     minl  = 1;
918:     label = NULL;
919:     if (enable_bmark) {
920:       PetscInt lminl = PETSC_MAX_INT;

922:       DMGetLabel(dm, bmark, &label);
923:       if (label) {
924:         IS       vals;
925:         PetscInt ldef;

927:         DMLabelGetDefaultValue(label, &ldef);
928:         DMLabelGetValueIS(label, &vals);
929:         ISGetMinMax(vals, &lminl, NULL);
930:         ISDestroy(&vals);
931:         lminl = PetscMin(ldef, lminl);
932:       }
933:       MPIU_Allreduce(&lminl, &minl, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
934:       if (minl == PETSC_MAX_INT) minl = 1;
935:     }
936:     PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", bf);
937:     for (p = fStart; p < fEnd; p++) {
938:       if (PetscUnlikely(PetscBTLookup(bfaces, p - fStart))) {
939:         const PetscInt *support;
940:         PetscInt        supportSize, c, nc = 0;

942:         DMPlexGetSupportSize(dm, p, &supportSize);
943:         DMPlexGetSupport(dm, p, &support);
944:         if (pown) {
945:           for (c = 0; c < supportSize; c++) {
946:             if (PetscLikely(PetscBTLookup(pown, support[c] - cStart))) fcells[nc++] = support[c];
947:           }
948:         } else
949:           for (c = 0; c < supportSize; c++) fcells[nc++] = support[c];
950:         for (c = 0; c < nc; c++) {
951:           const DMPolytopeType *faceTypes;
952:           DMPolytopeType        cellType;
953:           const PetscInt       *faceSizes, *cone;
954:           PetscInt              vids[8], *faces, st, i, coneSize, cell, cl, nv, cid = -1, mid = -1;

956:           cell = fcells[c];
957:           DMPlexGetCone(dm, cell, &cone);
958:           DMPlexGetConeSize(dm, cell, &coneSize);
959:           for (cl = 0; cl < coneSize; cl++)
960:             if (cone[cl] == p) break;
961:           if (cl == coneSize) continue;

963:           /* face material id and type */
964:           DMPlexGetPointMFEMCellID_Internal(dm, label, minl, p, &mid, &cid);
965:           PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT, mid, cid);
966:           /* vertex ids */
967:           DMPlexGetCellType(dm, cell, &cellType);
968:           DMPlexGetPointMFEMVertexIDs_Internal(dm, cell, (localized && !hovec) ? coordSection : NULL, &nv, vids);
969:           DMPlexGetRawFaces_Internal(dm, cellType, vids, NULL, &faceTypes, &faceSizes, (const PetscInt **)&faces);
970:           st = 0;
971:           for (i = 0; i < cl; i++) st += faceSizes[i];
972:           DMPlexInvertCell(faceTypes[cl], faces + st);
973:           for (i = 0; i < faceSizes[cl]; i++) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, faces[st + i]);
974:           PetscViewerASCIIPrintf(viewer, "\n");
975:           DMPlexRestoreRawFaces_Internal(dm, cellType, vids, NULL, &faceTypes, &faceSizes, (const PetscInt **)&faces);
976:           bf -= 1;
977:         }
978:       }
979:     }
981:     PetscBTDestroy(&bfaces);
982:     PetscFree(fcells);
983:   }

985:   /* mark owned vertices */
986:   vown = NULL;
987:   if (pown) {
988:     PetscBTCreate(vEnd - vStart, &vown);
989:     for (p = cStart; p < cEnd; p++) {
990:       PetscInt i, closureSize, *closure = NULL;

992:       if (PetscUnlikely(!PetscBTLookup(pown, p - cStart))) continue;
993:       DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
994:       for (i = 0; i < closureSize; i++) {
995:         const PetscInt pp = closure[2 * i];

997:         if (pp >= vStart && pp < vEnd) PetscBTSet(vown, pp - vStart);
998:       }
999:       DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
1000:     }
1001:   }

1003:   if (parentSection) {
1004:     PetscInt vp, gvp;

1006:     for (vp = 0, p = vStart; p < vEnd; p++) {
1007:       DMLabel  dlabel;
1008:       PetscInt parent, depth;

1010:       if (PetscUnlikely(vown && !PetscBTLookup(vown, p - vStart))) continue;
1011:       DMPlexGetDepthLabel(dm, &dlabel);
1012:       DMLabelGetValue(dlabel, p, &depth);
1013:       DMPlexGetTreeParent(dm, p, &parent, NULL);
1014:       if (parent != p) vp++;
1015:     }
1016:     MPIU_Allreduce(&vp, &gvp, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm));
1017:     if (gvp) {
1018:       PetscInt   maxsupp;
1019:       PetscBool *skip = NULL;

1021:       PetscViewerASCIIPrintf(viewer, "\nvertex_parents\n");
1022:       PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", vp);
1023:       DMPlexGetMaxSizes(dm, NULL, &maxsupp);
1024:       PetscMalloc1(maxsupp, &skip);
1025:       for (p = vStart; p < vEnd; p++) {
1026:         DMLabel  dlabel;
1027:         PetscInt parent;

1029:         if (PetscUnlikely(vown && !PetscBTLookup(vown, p - vStart))) continue;
1030:         DMPlexGetDepthLabel(dm, &dlabel);
1031:         DMPlexGetTreeParent(dm, p, &parent, NULL);
1032:         if (parent != p) {
1033:           PetscInt        vids[8] = {-1, -1, -1, -1, -1, -1, -1, -1}; /* silent overzealous clang static analyzer */
1034:           PetscInt        i, nv, ssize, n, numChildren, depth = -1;
1035:           const PetscInt *children;

1037:           DMPlexGetConeSize(dm, parent, &ssize);
1038:           switch (ssize) {
1039:           case 2: /* edge */
1040:             nv = 0;
1041:             DMPlexGetPointMFEMVertexIDs_Internal(dm, parent, localized ? coordSection : NULL, &nv, vids);
1042:             PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, p - vStart);
1043:             for (i = 0; i < nv; i++) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, vids[i]);
1044:             PetscViewerASCIIPrintf(viewer, "\n");
1045:             vp--;
1046:             break;
1047:           case 4: /* face */
1048:             DMPlexGetTreeChildren(dm, parent, &numChildren, &children);
1049:             for (n = 0; n < numChildren; n++) {
1050:               DMLabelGetValue(dlabel, children[n], &depth);
1051:               if (!depth) {
1052:                 const PetscInt *hvsupp, *hesupp, *cone;
1053:                 PetscInt        hvsuppSize, hesuppSize, coneSize;
1054:                 PetscInt        hv = children[n], he = -1, f;

1056:                 PetscArrayzero(skip, maxsupp);
1057:                 DMPlexGetSupportSize(dm, hv, &hvsuppSize);
1058:                 DMPlexGetSupport(dm, hv, &hvsupp);
1059:                 for (i = 0; i < hvsuppSize; i++) {
1060:                   PetscInt ep;
1061:                   DMPlexGetTreeParent(dm, hvsupp[i], &ep, NULL);
1062:                   if (ep != hvsupp[i]) {
1063:                     he = hvsupp[i];
1064:                   } else {
1065:                     skip[i] = PETSC_TRUE;
1066:                   }
1067:                 }
1069:                 DMPlexGetCone(dm, he, &cone);
1070:                 vids[0] = (cone[0] == hv) ? cone[1] : cone[0];
1071:                 DMPlexGetSupportSize(dm, he, &hesuppSize);
1072:                 DMPlexGetSupport(dm, he, &hesupp);
1073:                 for (f = 0; f < hesuppSize; f++) {
1074:                   PetscInt j;

1076:                   DMPlexGetCone(dm, hesupp[f], &cone);
1077:                   DMPlexGetConeSize(dm, hesupp[f], &coneSize);
1078:                   for (j = 0; j < coneSize; j++) {
1079:                     PetscInt k;
1080:                     for (k = 0; k < hvsuppSize; k++) {
1081:                       if (hvsupp[k] == cone[j]) {
1082:                         skip[k] = PETSC_TRUE;
1083:                         break;
1084:                       }
1085:                     }
1086:                   }
1087:                 }
1088:                 for (i = 0; i < hvsuppSize; i++) {
1089:                   if (!skip[i]) {
1090:                     DMPlexGetCone(dm, hvsupp[i], &cone);
1091:                     vids[1] = (cone[0] == hv) ? cone[1] : cone[0];
1092:                   }
1093:                 }
1094:                 PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, hv - vStart);
1095:                 for (i = 0; i < 2; i++) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, vids[i] - vStart);
1096:                 PetscViewerASCIIPrintf(viewer, "\n");
1097:                 vp--;
1098:               }
1099:             }
1100:             break;
1101:           default:
1102:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Don't know how to deal with support size %" PetscInt_FMT, ssize);
1103:           }
1104:         }
1105:       }
1106:       PetscFree(skip);
1107:     }
1109:   }
1110:   PetscBTDestroy(&vown);

1112:   /* vertices */
1113:   if (hovec) { /* higher-order meshes */
1114:     const char *fec;
1115:     PetscInt    i, n, s;
1116:     PetscViewerASCIIPrintf(viewer, "\nvertices\n");
1117:     PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", vEnd - vStart);
1118:     PetscViewerASCIIPrintf(viewer, "nodes\n");
1119:     PetscObjectGetName((PetscObject)hovec, &fec);
1120:     PetscViewerASCIIPrintf(viewer, "FiniteElementSpace\n");
1121:     PetscViewerASCIIPrintf(viewer, "%s\n", fec);
1122:     PetscViewerASCIIPrintf(viewer, "VDim: %" PetscInt_FMT "\n", sdim);
1123:     PetscViewerASCIIPrintf(viewer, "Ordering: 1\n\n"); /*Ordering::byVDIM*/
1124:     if (hoSection) {
1125:       DM cdm;

1127:       VecGetDM(hovec, &cdm);
1128:       for (p = cStart; p < cEnd; p++) {
1129:         PetscScalar *vals = NULL;
1130:         PetscInt     csize;

1132:         if (PetscUnlikely(pown && !PetscBTLookup(pown, p - cStart))) continue;
1133:         DMPlexVecGetClosure(cdm, hoSection, hovec, p, &csize, &vals);
1135:         for (i = 0; i < csize / sdim; i++) {
1136:           for (s = 0; s < sdim; s++) PetscViewerASCIIPrintf(viewer, fmt, (double)PetscRealPart(vals[i * sdim + s]));
1137:           PetscViewerASCIIPrintf(viewer, "\n");
1138:         }
1139:         DMPlexVecRestoreClosure(cdm, hoSection, hovec, p, &csize, &vals);
1140:       }
1141:     } else {
1142:       VecGetArrayRead(hovec, &array);
1143:       VecGetLocalSize(hovec, &n);
1145:       for (i = 0; i < n / sdim; i++) {
1146:         for (s = 0; s < sdim; s++) PetscViewerASCIIPrintf(viewer, fmt, (double)PetscRealPart(array[i * sdim + s]));
1147:         PetscViewerASCIIPrintf(viewer, "\n");
1148:       }
1149:       VecRestoreArrayRead(hovec, &array);
1150:     }
1151:   } else {
1152:     VecGetLocalSize(coordinates, &nvert);
1153:     PetscViewerASCIIPrintf(viewer, "\nvertices\n");
1154:     PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", nvert / sdim);
1155:     PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", sdim);
1156:     VecGetArrayRead(coordinates, &array);
1157:     for (p = 0; p < nvert / sdim; p++) {
1158:       PetscInt s;
1159:       for (s = 0; s < sdim; s++) {
1160:         PetscReal v = PetscRealPart(array[p * sdim + s]);

1162:         PetscViewerASCIIPrintf(viewer, fmt, PetscIsInfOrNanReal(v) ? 0.0 : (double)v);
1163:       }
1164:       PetscViewerASCIIPrintf(viewer, "\n");
1165:     }
1166:     VecRestoreArrayRead(coordinates, &array);
1167:   }
1168:   PetscBTDestroy(&pown);
1169:   VecDestroy(&hovec);
1170:   return 0;
1171: }

1173: PetscErrorCode DMPlexView_GLVis(DM dm, PetscViewer viewer)
1174: {
1175:   DMView_GLVis(dm, viewer, DMPlexView_GLVis_ASCII);
1176:   return 0;
1177: }