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

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

 32:   for (f=0;f<nf;f++) {
 33:     VecScatterBegin(ctx->scctx[f],(Vec)oX,(Vec)oXfield[f],INSERT_VALUES,SCATTER_FORWARD);
 34:     VecScatterEnd(ctx->scctx[f],(Vec)oX,(Vec)oXfield[f],INSERT_VALUES,SCATTER_FORWARD);
 35:   }
 36:   return 0;
 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:   DMGetDimension(dm,&dim);
 56:   DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);
 57:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
 58:   PetscObjectQuery((PetscObject)dm,"_glvis_plex_gnum",(PetscObject*)&globalNum);
 59:   if (!globalNum) {
 60:     DMPlexCreateCellNumbering_Internal(dm,PETSC_TRUE,&globalNum);
 61:     PetscObjectCompose((PetscObject)dm,"_glvis_plex_gnum",(PetscObject)globalNum);
 62:     PetscObjectDereference((PetscObject)globalNum);
 63:   }
 64:   ISGetIndices(globalNum,&gNum);
 65:   PetscBTCreate(vEnd-vStart,&vown);
 66:   for (c = cStart, totc = 0; c < cEnd; c++) {
 67:     if (gNum[c-cStart] >= 0) {
 68:       PetscInt i,numPoints,*points = NULL;

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

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

 89:     PetscSectionGetFieldComponents(s,f,&bs);
 90:     maxfields += bs;
 91:   }
 92:   PetscCalloc7(maxfields,&fieldname,maxfields,&nlocal,maxfields,&bs,maxfields,&dims,maxfields,&fec_type,totdofs,&idxs,maxfields,&Ufield);
 93:   PetscNew(&ctx);
 94:   PetscCalloc1(maxfields,&ctx->scctx);
 95:   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:       PetscSectionGetFieldName(s,f,&fname);
104:       PetscStrlen(fname,&len);
105:       if (len) {
106:         PetscStrcpy(name,fname);
107:       } else {
108:         PetscSNPrintf(name,256,"Field%D",f);
109:       }
110:       PetscDSGetDiscretization(ds,f,&disc);
111:       if (disc) {
112:         PetscClassId id;
113:         PetscInt     Nc;
114:         char         fec[64];

116:         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:           PetscFEGetNumComponents(fem,&Nc);
125:           PetscFEGetDualSpace(fem,&sp);
126:           PetscDualSpaceGetType(sp,&spname);
127:           PetscStrcmp(spname,PETSCDUALSPACELAGRANGE,&islag);
129:           PetscDualSpaceLagrangeGetContinuity(sp,&continuous);
130:           PetscDualSpaceGetOrder(sp,&order);
131:           if (continuous && order > 0) { /* no support for high-order viz, still have to figure out the numbering */
132:             PetscSNPrintf(fec,64,"FiniteElementCollection: H1_%DD_P1",dim);
133:           } else {
135:             H1   = PETSC_FALSE;
136:             PetscSNPrintf(fec,64,"FiniteElementCollection: L2_%DD_P%D",dim,order);
137:           }
138:           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:             PetscStrallocpy(fec,&fec_type[ctx->nf]);
144:             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:               PetscSectionGetFieldOffset(s,i+vStart,f,&off);
150:               for (j=0;j<Nc;j++) idxs[cum++] = off + j;
151:             }
152:             ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),Nv*Nc,idxs,PETSC_USE_POINTER,&isfield);
153:           } else {
154:             nlocal[ctx->nf] = Nc * totc;
155:             PetscStrallocpy(fec,&fec_type[ctx->nf]);
156:             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:               PetscSectionGetFieldOffset(s,i+cStart,f,&off);
162:               for (j=0;j<Nc;j++) idxs[cum++] = off + j;
163:             }
164:             ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),totc*Nc,idxs,PETSC_USE_POINTER,&isfield);
165:           }
166:           VecScatterCreate(xlocal,isfield,xfield,NULL,&ctx->scctx[ctx->nf]);
167:           VecDestroy(&xfield);
168:           ISDestroy(&isfield);
169:           ctx->nf++;
170:         } else if (id == PETSCFV_CLASSID) {
171:           PetscInt c;

173:           PetscFVGetNumComponents((PetscFV)disc,&Nc);
174:           PetscSNPrintf(fec,64,"FiniteElementCollection: L2_%DD_P0",dim);
175:           for (c = 0; c < Nc; c++) {
176:             char comp[256];
177:             PetscSNPrintf(comp,256,"%s-Comp%D",name,c);
178:             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:             PetscStrallocpy(fec,&fec_type[ctx->nf]);
183:             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:               PetscSectionGetFieldOffset(s,i+cStart,f,&off);
189:               idxs[cum++] = off + c;
190:             }
191:             ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),totc,idxs,PETSC_USE_POINTER,&isfield);
192:             VecScatterCreate(xlocal,isfield,xfield,NULL,&ctx->scctx[ctx->nf]);
193:             VecDestroy(&xfield);
194:             ISDestroy(&isfield);
195:             ctx->nf++;
196:           }
197:         } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %D",f);
198:       } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing discretization for field %D",f);
199:     }
200:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Needs a DS attached to the DM");
201:   PetscBTDestroy(&vown);
202:   VecDestroy(&xlocal);
203:   ISRestoreIndices(globalNum,&gNum);

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

213:   /* customize the viewer */
214:   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:     PetscFree(fieldname[f]);
217:     PetscFree(fec_type[f]);
218:     VecDestroy(&Ufield[f]);
219:   }
220:   PetscFree7(fieldname,nlocal,bs,dims,fec_type,idxs,Ufield);
221:   return 0;
222: }

224: typedef enum {MFEM_POINT=0,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE,MFEM_TETRAHEDRON,MFEM_CUBE,MFEM_PRISM,MFEM_UNDEF} MFEM_cid;

226: MFEM_cid mfem_table_cid[4][7]       = { {MFEM_POINT,MFEM_UNDEF,MFEM_UNDEF  ,MFEM_UNDEF   ,MFEM_UNDEF      ,MFEM_UNDEF,MFEM_UNDEF},
227:                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF   ,MFEM_UNDEF      ,MFEM_UNDEF,MFEM_UNDEF},
228:                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE     ,MFEM_UNDEF,MFEM_UNDEF},
229:                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF   ,MFEM_TETRAHEDRON,MFEM_PRISM,MFEM_CUBE } };

231: MFEM_cid mfem_table_cid_unint[4][9] = { {MFEM_POINT,MFEM_UNDEF,MFEM_UNDEF  ,MFEM_UNDEF   ,MFEM_UNDEF      ,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_UNDEF},
232:                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF   ,MFEM_UNDEF      ,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_UNDEF},
233:                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE     ,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_UNDEF},
234:                                         {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF   ,MFEM_TETRAHEDRON,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_CUBE } };

236: static PetscErrorCode DMPlexGetPointMFEMCellID_Internal(DM dm, DMLabel label, PetscInt minl, PetscInt p, PetscInt *mid, PetscInt *cid)
237: {
238:   DMLabel        dlabel;
239:   PetscInt       depth,csize,pdepth,dim;

241:   DMPlexGetDepthLabel(dm,&dlabel);
242:   DMLabelGetValue(dlabel,p,&pdepth);
243:   DMPlexGetConeSize(dm,p,&csize);
244:   DMPlexGetDepth(dm,&depth);
245:   DMGetDimension(dm,&dim);
246:   if (label) {
247:     DMLabelGetValue(label,p,mid);
248:     *mid = *mid - minl + 1; /* MFEM does not like negative markers */
249:   } else *mid = 1;
250:   if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
254:     *cid = mfem_table_cid_unint[dim][csize];
255:   } else {
258:     *cid = mfem_table_cid[pdepth][csize];
259:   }
260:   return 0;
261: }

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

267:   DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);
268:   DMGetDimension(dm,&dim);
269:   sdim = dim;
270:   if (csec) {
271:     PetscInt sStart,sEnd;

273:     DMGetCoordinateDim(dm,&sdim);
274:     PetscSectionGetChart(csec,&sStart,&sEnd);
275:     PetscSectionGetOffset(csec,vStart,&off);
276:     off  = off/sdim;
277:     if (p >= sStart && p < sEnd) {
278:       PetscSectionGetDof(csec,p,&dof);
279:     }
280:   }
281:   if (!dof) {
282:     DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);
283:     for (i=0,q=0;i<numPoints*2;i+= 2)
284:       if ((points[i] >= vStart) && (points[i] < vEnd))
285:         vids[q++] = points[i]-vStart+off;
286:     DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);
287:   } else {
288:     PetscSectionGetOffset(csec,p,&off);
289:     PetscSectionGetDof(csec,p,&dof);
290:     for (q=0;q<dof/sdim;q++) vids[q] = off/sdim + q;
291:   }
292:   *nv = q;
293:   return 0;
294: }

296: static PetscErrorCode GLVisCreateFE(PetscFE femIn,char name[32],PetscFE *fem,IS *perm)
297: {
298:   DM              K;
299:   PetscSpace      P;
300:   PetscDualSpace  Q;
301:   PetscQuadrature q,fq;
302:   PetscInt        dim,deg,dof;
303:   DMPolytopeType  ptype;
304:   PetscBool       isSimplex,isTensor;
305:   PetscBool       continuity = PETSC_FALSE;
306:   PetscDTNodeType nodeType   = PETSCDTNODES_GAUSSJACOBI;
307:   PetscBool       endpoint   = PETSC_TRUE;
308:   MPI_Comm        comm;

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

367:   /* Both MFEM and PETSc are lexicographic, but PLEX stores the swapped cone */
368:   *perm = NULL;
369:   if (isSimplex && dim == 3) {
370:     PetscInt celldofs,*pidx;

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

423:   /* Cleanup */
424:   PetscSpaceDestroy(&P);
425:   PetscDualSpaceDestroy(&Q);
426:   PetscQuadratureDestroy(&q);
427:   PetscQuadratureDestroy(&fq);
428:   return 0;
429: }

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

456:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
458:   MPI_Comm_size(PetscObjectComm((PetscObject)viewer),&size);
460:   DMGetDimension(dm,&dim);
461:   DMPlexGetDepth(dm,&depth);

463:   /* get container: determines if a process visualizes is portion of the data or not */
464:   PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);
466:   {
467:     PetscViewerGLVisInfo glvis_info;
468:     PetscContainerGetPointer(glvis_container,(void**)&glvis_info);
469:     enabled = glvis_info->enabled;
470:     fmt     = glvis_info->fmt;
471:   }

473:   /* Users can attach a coordinate vector to the DM in case they have a higher-order mesh */
474:   PetscObjectQuery((PetscObject)dm,"_glvis_mesh_coords",(PetscObject*)&hovec);
475:   PetscObjectReference((PetscObject)hovec);
476:   if (!hovec) {
477:     DM           cdm;
478:     PetscFE      disc;
479:     PetscClassId classid;

481:     DMGetCoordinateDM(dm,&cdm);
482:     DMGetField(cdm,0,NULL,(PetscObject*)&disc);
483:     PetscObjectGetClassId((PetscObject)disc,&classid);
484:     if (classid == PETSCFE_CLASSID) {
485:       DM      hocdm;
486:       PetscFE hodisc;
487:       Vec     vec;
488:       Mat     mat;
489:       char    name[32],fec_type[64];
490:       IS      perm = NULL;

492:       GLVisCreateFE(disc,name,&hodisc,&perm);
493:       DMClone(cdm,&hocdm);
494:       DMSetField(hocdm,0,NULL,(PetscObject)hodisc);
495:       PetscFEDestroy(&hodisc);
496:       DMCreateDS(hocdm);

498:       DMGetCoordinates(dm,&vec);
499:       DMCreateGlobalVector(hocdm,&hovec);
500:       DMCreateInterpolation(cdm,hocdm,&mat,NULL);
501:       MatInterpolate(mat,vec,hovec);
502:       MatDestroy(&mat);
503:       DMGetLocalSection(hocdm,&hoSection);
504:       PetscSectionSetClosurePermutation(hoSection, (PetscObject)hocdm, depth, perm);
505:       ISDestroy(&perm);
506:       DMDestroy(&hocdm);
507:       PetscSNPrintf(fec_type,sizeof(fec_type),"FiniteElementCollection: %s", name);
508:       PetscObjectSetName((PetscObject)hovec,fec_type);
509:     }
510:   }

512:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
513:   DMPlexGetGhostCellStratum(dm,&p,NULL);
514:   if (p >= 0) cEnd = p;
515:   DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);
516:   DMGetPeriodicity(dm,&periodic,NULL,NULL,NULL);
517:   DMGetCoordinatesLocalized(dm,&localized);
519:   DMGetCoordinateSection(dm,&coordSection);
520:   DMGetCoordinateDim(dm,&sdim);
521:   DMGetCoordinatesLocal(dm,&coordinates);

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

546:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
549:                                                              "Alternatively, run with -viewer_glvis_dm_plex_enable_boundary 0");
551:                                                            "Alternatively, run with -viewer_glvis_dm_plex_enable_ncmesh 0");
552:   if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
554:     cellvertex = PETSC_TRUE;
555:   }

557:   /* Identify possible cells in the overlap */
558:   novl = 0;
559:   pown = NULL;
560:   if (size > 1) {
561:     IS             globalNum = NULL;
562:     const PetscInt *gNum;
563:     PetscBool      ovl  = PETSC_FALSE;

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

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

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

629:       PetscSNPrintf(fec,sizeof(fec),"FiniteElementCollection: L2_T1_%DD_P1",dim);
630:       if (cEnd-cStart) {
631:         PetscInt fpc;

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

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

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

707:   /* topological dimension */
708:   PetscViewerASCIIPrintf(viewer,"\ndimension\n");
709:   PetscViewerASCIIPrintf(viewer,"%D\n",dim);

711:   /* elements */
712:   minl = 1;
713:   label = NULL;
714:   if (enable_emark) {
715:     PetscInt lminl = PETSC_MAX_INT;

717:     DMGetLabel(dm,emark,&label);
718:     if (label) {
719:       IS       vals;
720:       PetscInt ldef;

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

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

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

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

770:         DMPlexGetCellType(dm,p,&cellType);
771:         PetscSectionGetDof(coordSection,p,&dof);
772:         if (dof) {
773:           PetscInt    uvpc, v,csize,cellClosureSize,*cellClosure = NULL,*vidxs = NULL;
774:           PetscScalar *vals = NULL;

776:           uvpc = DMPolytopeTypeGetNumVertices(cellType);
778:           DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);
779:           DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);
780:           for (v=0;v<cellClosureSize;v++)
781:             if (cellClosure[2*v] >= vStart && cellClosure[2*v] < vEnd) {
782:               vidxs = cellClosure + 2*v;
783:               break;
784:             }
786:           for (v=0;v<uvpc;v++) {
787:             PetscInt s;

789:             for (s=0;s<sdim;s++) {
790:               if (PetscAbsScalar(vals[v*sdim+s]-vals[v*sdim+s+uvpc*sdim])>PETSC_MACHINE_EPSILON) {
791:                 DMLabelSetValue(perLabel,vidxs[2*v],2);
792:               }
793:             }
794:           }
795:           DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);
796:           DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);
797:         }
798:       }
799:       if (dim > 1) {
800:         PetscInt eEnd,eStart;

802:         DMPlexGetDepthStratum(dm,1,&eStart,&eEnd);
803:         for (p=eStart;p<eEnd;p++) {
804:           const PetscInt *cone;
805:           PetscInt       coneSize,i;
806:           PetscBool      ispe = PETSC_TRUE;

808:           DMPlexGetCone(dm,p,&cone);
809:           DMPlexGetConeSize(dm,p,&coneSize);
810:           for (i=0;i<coneSize;i++) {
811:             PetscInt v;

813:             DMLabelGetValue(perLabel,cone[i],&v);
814:             ispe = (PetscBool)(ispe && (v==2));
815:           }
816:           if (ispe && coneSize) {
817:             PetscInt       ch, numChildren;
818:             const PetscInt *children;

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

833:             DMPlexGetCone(dm,p,&cone);
834:             DMPlexGetConeSize(dm,p,&coneSize);
835:             for (i=0;i<coneSize;i++) {
836:               PetscInt v;

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

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

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

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

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

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

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

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

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

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

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

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

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

990:   /* mark owned vertices */
991:   vown = NULL;
992:   if (pown) {
993:     PetscBTCreate(vEnd-vStart,&vown);
994:     for (p=cStart;p<cEnd;p++) {
995:       PetscInt i,closureSize,*closure = NULL;

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

1002:         if (pp >= vStart && pp < vEnd) {
1003:           PetscBTSet(vown,pp-vStart);
1004:         }
1005:       }
1006:       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:       DMPlexGetDepthLabel(dm,&dlabel);
1019:       DMLabelGetValue(dlabel,p,&depth);
1020:       DMPlexGetTreeParent(dm,p,&parent,NULL);
1021:       if (parent != p) vp++;
1022:     }
1023:     MPIU_Allreduce(&vp,&gvp,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
1024:     if (gvp) {
1025:       PetscInt  maxsupp;
1026:       PetscBool *skip = NULL;

1028:       PetscViewerASCIIPrintf(viewer,"\nvertex_parents\n");
1029:       PetscViewerASCIIPrintf(viewer,"%D\n",vp);
1030:       DMPlexGetMaxSizes(dm,NULL,&maxsupp);
1031:       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:         DMPlexGetDepthLabel(dm,&dlabel);
1038:         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:           DMPlexGetConeSize(dm,parent,&ssize);
1045:           switch (ssize) {
1046:             case 2: /* edge */
1047:               nv   = 0;
1048:               DMPlexGetPointMFEMVertexIDs_Internal(dm,parent,localized ? coordSection : NULL,&nv,vids);
1049:               PetscViewerASCIIPrintf(viewer,"%D",p-vStart);
1050:               for (i=0;i<nv;i++) {
1051:                 PetscViewerASCIIPrintf(viewer," %D",vids[i]);
1052:               }
1053:               PetscViewerASCIIPrintf(viewer,"\n");
1054:               vp--;
1055:               break;
1056:             case 4: /* face */
1057:               DMPlexGetTreeChildren(dm,parent,&numChildren,&children);
1058:               for (n=0;n<numChildren;n++) {
1059:                 DMLabelGetValue(dlabel,children[n],&depth);
1060:                 if (!depth) {
1061:                   const PetscInt *hvsupp,*hesupp,*cone;
1062:                   PetscInt       hvsuppSize,hesuppSize,coneSize;
1063:                   PetscInt       hv = children[n],he = -1,f;

1065:                   PetscArrayzero(skip,maxsupp);
1066:                   DMPlexGetSupportSize(dm,hv,&hvsuppSize);
1067:                   DMPlexGetSupport(dm,hv,&hvsupp);
1068:                   for (i=0;i<hvsuppSize;i++) {
1069:                     PetscInt ep;
1070:                     DMPlexGetTreeParent(dm,hvsupp[i],&ep,NULL);
1071:                     if (ep != hvsupp[i]) {
1072:                       he = hvsupp[i];
1073:                     } else {
1074:                       skip[i] = PETSC_TRUE;
1075:                     }
1076:                   }
1078:                   DMPlexGetCone(dm,he,&cone);
1079:                   vids[0] = (cone[0] == hv) ? cone[1] : cone[0];
1080:                   DMPlexGetSupportSize(dm,he,&hesuppSize);
1081:                   DMPlexGetSupport(dm,he,&hesupp);
1082:                   for (f=0;f<hesuppSize;f++) {
1083:                     PetscInt j;

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

1123:   /* vertices */
1124:   if (hovec) { /* higher-order meshes */
1125:     const char *fec;
1126:     PetscInt   i,n,s;
1127:     PetscViewerASCIIPrintf(viewer,"\nvertices\n");
1128:     PetscViewerASCIIPrintf(viewer,"%D\n",vEnd-vStart);
1129:     PetscViewerASCIIPrintf(viewer,"nodes\n");
1130:     PetscObjectGetName((PetscObject)hovec,&fec);
1131:     PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");
1132:     PetscViewerASCIIPrintf(viewer,"%s\n",fec);
1133:     PetscViewerASCIIPrintf(viewer,"VDim: %D\n",sdim);
1134:     PetscViewerASCIIPrintf(viewer,"Ordering: 1\n\n"); /*Ordering::byVDIM*/
1135:     if (hoSection) {
1136:       DM cdm;

1138:       VecGetDM(hovec,&cdm);
1139:       for (p=cStart;p<cEnd;p++) {
1140:         PetscScalar *vals = NULL;
1141:         PetscInt    csize;

1143:         if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue;
1144:         DMPlexVecGetClosure(cdm,hoSection,hovec,p,&csize,&vals);
1146:         for (i=0;i<csize/sdim;i++) {
1147:           for (s=0;s<sdim;s++) {
1148:             PetscViewerASCIIPrintf(viewer,fmt,(double) PetscRealPart(vals[i*sdim+s]));
1149:           }
1150:           PetscViewerASCIIPrintf(viewer,"\n");
1151:         }
1152:         DMPlexVecRestoreClosure(cdm,hoSection,hovec,p,&csize,&vals);
1153:       }
1154:     } else {
1155:       VecGetArrayRead(hovec,&array);
1156:       VecGetLocalSize(hovec,&n);
1158:       for (i=0;i<n/sdim;i++) {
1159:         for (s=0;s<sdim;s++) {
1160:           PetscViewerASCIIPrintf(viewer,fmt,(double) PetscRealPart(array[i*sdim+s]));
1161:         }
1162:         PetscViewerASCIIPrintf(viewer,"\n");
1163:       }
1164:       VecRestoreArrayRead(hovec,&array);
1165:     }
1166:   } else {
1167:     VecGetLocalSize(coordinates,&nvert);
1168:     PetscViewerASCIIPrintf(viewer,"\nvertices\n");
1169:     PetscViewerASCIIPrintf(viewer,"%D\n",nvert/sdim);
1170:     PetscViewerASCIIPrintf(viewer,"%D\n",sdim);
1171:     VecGetArrayRead(coordinates,&array);
1172:     for (p=0;p<nvert/sdim;p++) {
1173:       PetscInt s;
1174:       for (s=0;s<sdim;s++) {
1175:         PetscReal v = PetscRealPart(array[p*sdim+s]);

1177:         PetscViewerASCIIPrintf(viewer,fmt,PetscIsInfOrNanReal(v) ? 0.0 : (double) v);
1178:       }
1179:       PetscViewerASCIIPrintf(viewer,"\n");
1180:     }
1181:     VecRestoreArrayRead(coordinates,&array);
1182:   }
1183:   PetscBTDestroy(&pown);
1184:   VecDestroy(&hovec);
1185:   return 0;
1186: }

1188: PetscErrorCode DMPlexView_GLVis(DM dm, PetscViewer viewer)
1189: {
1190:   DMView_GLVis(dm,viewer,DMPlexView_GLVis_ASCII);
1191:   return 0;
1192: }