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