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