Actual source code: plexglvis.c
petsc-3.11.4 2019-09-28
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: DMGetSection(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: if (!continuous) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Discontinuous space visualization currently unsupported");
137: PetscDualSpaceGetOrder(sp,&order);
138: if (continuous && order > 0) {
139: PetscSNPrintf(fec,64,"FiniteElementCollection: H1_%DD_P1",dim);
140: } else {
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;
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: PetscOptionsString("-viewer_glvis_dm_plex_emarker","String for the material id label",NULL,emark,emark,sizeof(emark),&enable_emark);
405: PetscOptionsString("-viewer_glvis_dm_plex_bmarker","String for the boundary id label",NULL,bmark,bmark,sizeof(bmark),&enable_bmark);
406: PetscOptionsEnd();
407: if (enable_bmark) enable_boundary = PETSC_TRUE;
409: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
410: if (enable_ncmesh && size > 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not supported in parallel");
411: DMPlexGetDepth(dm,&depth);
412: if (enable_boundary && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. "
413: "Alternatively, run with -viewer_glvis_dm_plex_enable_boundary 0");
414: if (enable_ncmesh && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. "
415: "Alternatively, run with -viewer_glvis_dm_plex_enable_ncmesh 0");
416: if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
417: if (depth != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported depth %D. You should interpolate the mesh first",depth);
418: cellvertex = PETSC_TRUE;
419: }
421: /* Identify possible cells in the overlap */
422: novl = 0;
423: pown = NULL;
424: if (size > 1) {
425: IS globalNum = NULL;
426: const PetscInt *gNum;
427: PetscBool ovl = PETSC_FALSE;
429: PetscObjectQuery((PetscObject)dm,"_glvis_plex_gnum",(PetscObject*)&globalNum);
430: if (!globalNum) {
431: DMPlexCreateCellNumbering_Internal(dm,PETSC_TRUE,&globalNum);
432: PetscObjectCompose((PetscObject)dm,"_glvis_plex_gnum",(PetscObject)globalNum);
433: PetscObjectDereference((PetscObject)globalNum);
434: }
435: ISGetIndices(globalNum,&gNum);
436: for (p=cStart; p<cEnd; p++) {
437: if (gNum[p-cStart] < 0) {
438: ovl = PETSC_TRUE;
439: novl++;
440: }
441: }
442: if (ovl) {
443: /* it may happen that pown get not destroyed, if the user closes the window while this function is running.
444: TODO: garbage collector? attach pown to dm? */
445: PetscBTCreate(cEnd-cStart,&pown);
446: for (p=cStart; p<cEnd; p++) {
447: if (gNum[p-cStart] < 0) continue;
448: else {
449: PetscBTSet(pown,p-cStart);
450: }
451: }
452: }
453: ISRestoreIndices(globalNum,&gNum);
454: }
456: /* return if this process is disabled */
457: if (!enabled) {
458: PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");
459: PetscViewerASCIIPrintf(viewer,"\ndimension\n");
460: PetscViewerASCIIPrintf(viewer,"%D\n",dim);
461: PetscViewerASCIIPrintf(viewer,"\nelements\n");
462: PetscViewerASCIIPrintf(viewer,"%D\n",0);
463: PetscViewerASCIIPrintf(viewer,"\nboundary\n");
464: PetscViewerASCIIPrintf(viewer,"%D\n",0);
465: PetscViewerASCIIPrintf(viewer,"\nvertices\n");
466: PetscViewerASCIIPrintf(viewer,"%D\n",0);
467: PetscViewerASCIIPrintf(viewer,"%D\n",sdim);
468: PetscBTDestroy(&pown);
469: return(0);
470: }
472: if (enable_mfem) {
473: if (periodic && !hovec) { /* we need to generate a vector of L2 coordinates, as this is how MFEM handles periodic meshes */
474: PetscInt vpc = 0;
475: char fec[64];
476: int vids[8] = {0,1,2,3,4,5,6,7};
477: int hexv[8] = {0,1,3,2,4,5,7,6}, tetv[4] = {0,1,2,3};
478: int quadv[8] = {0,1,3,2}, triv[3] = {0,1,2};
479: int *dof = NULL;
480: PetscScalar *array,*ptr;
482: PetscSNPrintf(fec,sizeof(fec),"FiniteElementCollection: L2_T1_%DD_P1",dim);
483: if (cEndInterior < cEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Support for hybrid meshed not currently implemented");
484: if (cEnd-cStart) {
485: PetscInt fpc;
487: DMPlexGetConeSize(dm,cStart,&fpc);
488: switch(dim) {
489: case 1:
490: vpc = 2;
491: dof = hexv;
492: break;
493: case 2:
494: switch (fpc) {
495: case 3:
496: vpc = 3;
497: dof = triv;
498: break;
499: case 4:
500: vpc = 4;
501: dof = quadv;
502: break;
503: default:
504: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
505: break;
506: }
507: break;
508: case 3:
509: switch (fpc) {
510: case 4: /* TODO: still need to understand L2 ordering for tets */
511: vpc = 4;
512: dof = tetv;
513: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled tethraedral case");
514: break;
515: case 6:
516: if (cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: vertices per cell %D",fpc);
517: vpc = 8;
518: dof = hexv;
519: break;
520: case 8:
521: if (!cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
522: vpc = 8;
523: dof = hexv;
524: break;
525: default:
526: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
527: break;
528: }
529: break;
530: default:
531: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim");
532: break;
533: }
534: DMPlexInvertCell(dim,vpc,vids);
535: }
536: if (!dof) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing dofs");
537: VecCreateSeq(PETSC_COMM_SELF,(cEnd-cStart-novl)*vpc*sdim,&hovec);
538: PetscObjectCompose((PetscObject)dm,"_glvis_mesh_coords",(PetscObject)hovec);
539: PetscObjectDereference((PetscObject)hovec);
540: PetscObjectSetName((PetscObject)hovec,fec);
541: VecGetArray(hovec,&array);
542: ptr = array;
543: for (p=cStart;p<cEnd;p++) {
544: PetscInt csize,v,d;
545: PetscScalar *vals = NULL;
547: if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue;
548: DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);
549: 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);
550: for (v=0;v<vpc;v++) {
551: for (d=0;d<sdim;d++) {
552: ptr[sdim*dof[v]+d] = vals[sdim*vids[v]+d];
553: }
554: }
555: ptr += vpc*sdim;
556: DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);
557: }
558: VecRestoreArray(hovec,&array);
559: }
560: }
561: /* if we have high-order coordinates in 3D, we need to specify the boundary */
562: if (hovec && dim == 3) enable_boundary = PETSC_TRUE;
564: /* header */
565: PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");
567: /* topological dimension */
568: PetscViewerASCIIPrintf(viewer,"\ndimension\n");
569: PetscViewerASCIIPrintf(viewer,"%D\n",dim);
571: /* elements */
572: minl = 1;
573: label = NULL;
574: if (enable_emark) {
575: PetscInt lminl = PETSC_MAX_INT;
577: DMGetLabel(dm,emark,&label);
578: if (label) {
579: IS vals;
580: PetscInt ldef;
582: DMLabelGetDefaultValue(label,&ldef);
583: DMLabelGetValueIS(label,&vals);
584: ISGetMinMax(vals,&lminl,NULL);
585: ISDestroy(&vals);
586: lminl = PetscMin(ldef,lminl);
587: }
588: MPIU_Allreduce(&lminl,&minl,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)dm));
589: if (minl == PETSC_MAX_INT) minl = 1;
590: }
591: PetscViewerASCIIPrintf(viewer,"\nelements\n");
592: PetscViewerASCIIPrintf(viewer,"%D\n",cEnd-cStart-novl);
593: for (p=cStart;p<cEnd;p++) {
594: int vids[8];
595: PetscInt i,nv = 0,cid = -1,mid = 1;
597: if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue;
598: DMPlexGetPointMFEMCellID_Internal(dm,label,minl,p,&mid,&cid);
599: DMPlexGetPointMFEMVertexIDs_Internal(dm,p,(localized && !hovec) ? coordSection : NULL,&nv,vids);
600: DMPlexInvertCell(dim,nv,vids);
601: if (p >= cEndInterior) {
602: DMPlexGlvisInvertHybrid(cid,vids);
603: }
604: PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);
605: for (i=0;i<nv;i++) {
606: PetscViewerASCIIPrintf(viewer," %D",(PetscInt)vids[i]);
607: }
608: PetscViewerASCIIPrintf(viewer,"\n");
609: }
611: /* boundary */
612: PetscViewerASCIIPrintf(viewer,"\nboundary\n");
613: if (!enable_boundary) {
614: PetscViewerASCIIPrintf(viewer,"%D\n",0);
615: } else {
616: DMLabel perLabel;
617: PetscBT bfaces;
618: PetscInt fStart,fEnd,*fcells;
619: PetscInt *faces = NULL,fpc = 0,vpf = 0, vpc = 0;
620: PetscInt *facesH = NULL,fpcH = 0,vpfH = 0, vpcH = 0;
621: PetscInt fv1[] = {0,1},
622: fv2tri[] = {0,1,
623: 1,2,
624: 2,0},
625: fv2quad[] = {0,1,
626: 1,2,
627: 2,3,
628: 3,0},
629: fv2quadH[] = {0,1,
630: 2,3,
631: 0,2,
632: 1,3},
633: fv3tet[] = {0,1,2,
634: 0,3,1,
635: 0,2,3,
636: 2,1,3},
637: fv3wedge[] = {0,1,2,-1,
638: 3,4,5,-1,
639: 0,1,3,4,
640: 1,2,4,5,
641: 2,0,5,3},
642: fv3hex[] = {0,1,2,3,
643: 4,5,6,7,
644: 0,3,5,4,
645: 2,1,7,6,
646: 3,2,6,5,
647: 0,4,7,1};
649: /* determine orientation of boundary mesh */
650: if (cEnd-cStart) {
651: if (cEndInterior < cEnd) {
652: DMPlexGetConeSize(dm,cEndInterior,&fpcH);
653: }
654: if (cEndInterior > cStart) {
655: DMPlexGetConeSize(dm,cStart,&fpc);
656: }
657: switch(dim) {
658: case 1:
659: if (fpc != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case faces per cell %D",fpc);
660: faces = fv1;
661: vpf = 1;
662: vpc = 2;
663: break;
664: case 2:
665: switch (fpc) {
666: case 0:
667: case 3:
668: faces = fv2tri;
669: vpf = 2;
670: vpc = 3;
671: if (fpcH == 4) {
672: facesH = fv2quadH;
673: vpfH = 2;
674: vpcH = 4;
675: } else if (fpcH) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case hybrid faces per cell %D",fpcH);
676: break;
677: case 4:
678: faces = fv2quad;
679: vpf = 2;
680: vpc = 4;
681: if (fpcH) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case hybrid faces per cell %D",fpcH);
682: break;
683: default:
684: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
685: break;
686: }
687: break;
688: case 3:
689: switch (fpc) {
690: case 0:
691: case 4:
692: faces = fv3tet;
693: vpf = 3;
694: vpc = 4;
695: if (fpcH == 5) {
696: facesH = fv3wedge;
697: vpfH = -4;
698: vpcH = 6;
699: } else if (fpcH) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case hybrid faces per cell %D",fpcH);
700: break;
701: case 6:
702: faces = fv3hex;
703: vpf = 4;
704: vpc = 8;
705: if (fpcH) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case hybrid faces per cell %D",fpcH);
706: break;
707: default:
708: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
709: break;
710: }
711: break;
712: default:
713: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim");
714: break;
715: }
716: }
717: DMPlexGetHeightStratum(dm,1,&fStart,&fEnd);
718: PetscBTCreate(fEnd-fStart,&bfaces);
719: DMPlexGetMaxSizes(dm,NULL,&p);
720: PetscMalloc1(p,&fcells);
721: DMGetLabel(dm,"glvis_periodic_cut",&perLabel);
722: if (!perLabel && localized) { /* this periodic cut can be moved up to DMPlex setup */
723: DMCreateLabel(dm,"glvis_periodic_cut");
724: DMGetLabel(dm,"glvis_periodic_cut",&perLabel);
725: DMLabelSetDefaultValue(perLabel,1);
726: for (p=cStart;p<cEnd;p++) {
727: PetscInt dof, uvpc;
729: PetscSectionGetDof(coordSection,p,&dof);
730: if (dof) {
731: PetscInt v,csize,cellClosureSize,*cellClosure = NULL,*vidxs = NULL;
732: PetscScalar *vals = NULL;
733: if (p < cEndInterior) uvpc = vpc;
734: else uvpc = vpcH;
735: if (dof%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Incompatible number of cell dofs %D and space dimension %D",dof,sdim);
736: 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);
737: DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);
738: DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);
739: for (v=0;v<cellClosureSize;v++)
740: if (cellClosure[2*v] >= vStart && cellClosure[2*v] < vEnd) {
741: vidxs = cellClosure + 2*v;
742: break;
743: }
744: if (!vidxs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing vertices");
745: for (v=0;v<uvpc;v++) {
746: PetscInt s;
748: for (s=0;s<sdim;s++) {
749: if (PetscAbsScalar(vals[v*sdim+s]-vals[v*sdim+s+uvpc*sdim])>PETSC_MACHINE_EPSILON) {
750: DMLabelSetValue(perLabel,vidxs[2*v],2);
751: }
752: }
753: }
754: DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);
755: DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);
756: }
757: }
758: if (dim > 1) {
759: PetscInt eEnd,eStart;
761: DMPlexGetDepthStratum(dm,1,&eStart,&eEnd);
762: for (p=eStart;p<eEnd;p++) {
763: const PetscInt *cone;
764: PetscInt coneSize,i;
765: PetscBool ispe = PETSC_TRUE;
767: DMPlexGetCone(dm,p,&cone);
768: DMPlexGetConeSize(dm,p,&coneSize);
769: for (i=0;i<coneSize;i++) {
770: PetscInt v;
772: DMLabelGetValue(perLabel,cone[i],&v);
773: ispe = (PetscBool)(ispe && (v==2));
774: }
775: if (ispe && coneSize) {
776: PetscInt ch, numChildren;
777: const PetscInt *children;
779: DMLabelSetValue(perLabel,p,2);
780: DMPlexGetTreeChildren(dm,p,&numChildren,&children);
781: for (ch = 0; ch < numChildren; ch++) {
782: DMLabelSetValue(perLabel,children[ch],2);
783: }
784: }
785: }
786: if (dim > 2) {
787: for (p=fStart;p<fEnd;p++) {
788: const PetscInt *cone;
789: PetscInt coneSize,i;
790: PetscBool ispe = PETSC_TRUE;
792: DMPlexGetCone(dm,p,&cone);
793: DMPlexGetConeSize(dm,p,&coneSize);
794: for (i=0;i<coneSize;i++) {
795: PetscInt v;
797: DMLabelGetValue(perLabel,cone[i],&v);
798: ispe = (PetscBool)(ispe && (v==2));
799: }
800: if (ispe && coneSize) {
801: PetscInt ch, numChildren;
802: const PetscInt *children;
804: DMLabelSetValue(perLabel,p,2);
805: DMPlexGetTreeChildren(dm,p,&numChildren,&children);
806: for (ch = 0; ch < numChildren; ch++) {
807: DMLabelSetValue(perLabel,children[ch],2);
808: }
809: }
810: }
811: }
812: }
813: }
814: for (p=fStart;p<fEnd;p++) {
815: const PetscInt *support;
816: PetscInt supportSize;
817: PetscBool isbf = PETSC_FALSE;
819: DMPlexGetSupportSize(dm,p,&supportSize);
820: if (pown) {
821: PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE;
822: PetscInt i;
824: DMPlexGetSupport(dm,p,&support);
825: for (i=0;i<supportSize;i++) {
826: if (PetscLikely(PetscBTLookup(pown,support[i]-cStart))) has_owned = PETSC_TRUE;
827: else has_ghost = PETSC_TRUE;
828: }
829: isbf = (PetscBool)((supportSize == 1 && has_owned) || (supportSize > 1 && has_owned && has_ghost));
830: } else {
831: isbf = (PetscBool)(supportSize == 1);
832: }
833: if (!isbf && perLabel) {
834: const PetscInt *cone;
835: PetscInt coneSize,i;
837: DMPlexGetCone(dm,p,&cone);
838: DMPlexGetConeSize(dm,p,&coneSize);
839: isbf = PETSC_TRUE;
840: for (i=0;i<coneSize;i++) {
841: PetscInt v,d;
843: DMLabelGetValue(perLabel,cone[i],&v);
844: DMLabelGetDefaultValue(perLabel,&d);
845: isbf = (PetscBool)(isbf && v != d);
846: }
847: }
848: if (isbf) {
849: PetscBTSet(bfaces,p-fStart);
850: }
851: }
852: /* count boundary faces */
853: for (p=fStart,bf=0;p<fEnd;p++) {
854: if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) {
855: const PetscInt *support;
856: PetscInt supportSize,c;
858: DMPlexGetSupportSize(dm,p,&supportSize);
859: DMPlexGetSupport(dm,p,&support);
860: for (c=0;c<supportSize;c++) {
861: const PetscInt *cone;
862: PetscInt cell,cl,coneSize;
864: cell = support[c];
865: if (pown && PetscUnlikely(!PetscBTLookup(pown,cell-cStart))) continue;
866: DMPlexGetCone(dm,cell,&cone);
867: DMPlexGetConeSize(dm,cell,&coneSize);
868: for (cl=0;cl<coneSize;cl++) {
869: if (cone[cl] == p) {
870: bf += 1;
871: break;
872: }
873: }
874: }
875: }
876: }
877: minl = 1;
878: label = NULL;
879: if (enable_bmark) {
880: PetscInt lminl = PETSC_MAX_INT;
882: DMGetLabel(dm,bmark,&label);
883: if (label) {
884: IS vals;
885: PetscInt ldef;
887: DMLabelGetDefaultValue(label,&ldef);
888: DMLabelGetValueIS(label,&vals);
889: ISGetMinMax(vals,&lminl,NULL);
890: ISDestroy(&vals);
891: lminl = PetscMin(ldef,lminl);
892: }
893: MPIU_Allreduce(&lminl,&minl,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)dm));
894: if (minl == PETSC_MAX_INT) minl = 1;
895: }
896: PetscViewerASCIIPrintf(viewer,"%D\n",bf);
897: for (p=fStart;p<fEnd;p++) {
898: if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) {
899: const PetscInt *support;
900: PetscInt supportSize,c,nc = 0;
902: DMPlexGetSupportSize(dm,p,&supportSize);
903: DMPlexGetSupport(dm,p,&support);
904: if (pown) {
905: for (c=0;c<supportSize;c++) {
906: if (PetscLikely(PetscBTLookup(pown,support[c]-cStart))) {
907: fcells[nc++] = support[c];
908: }
909: }
910: } else for (c=0;c<supportSize;c++) fcells[nc++] = support[c];
911: for (c=0;c<nc;c++) {
912: const PetscInt *cone;
913: int vids[8];
914: PetscInt i,coneSize,cell,cl,nv,cid = -1,mid = -1;
916: cell = fcells[c];
917: DMPlexGetCone(dm,cell,&cone);
918: DMPlexGetConeSize(dm,cell,&coneSize);
919: for (cl=0;cl<coneSize;cl++)
920: if (cone[cl] == p)
921: break;
922: if (cl == coneSize) continue;
924: /* face material id and type */
925: DMPlexGetPointMFEMCellID_Internal(dm,label,minl,p,&mid,&cid);
926: PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);
927: /* vertex ids */
928: DMPlexGetPointMFEMVertexIDs_Internal(dm,cell,(localized && !hovec) ? coordSection : NULL,&nv,vids);
929: if (cell >= cEndInterior) {
930: PetscInt nv = vpfH, inc = vpfH;
931: if (vpfH < 0) { /* Wedge */
932: if (cl == 0 || cl == 1) nv = 3;
933: else nv = 4;
934: inc = -vpfH;
935: }
936: for (i=0;i<nv;i++) {
937: PetscViewerASCIIPrintf(viewer," %d",vids[facesH[cl*inc+i]]);
938: }
939: } else {
940: for (i=0;i<vpf;i++) {
941: PetscViewerASCIIPrintf(viewer," %d",vids[faces[cl*vpf+i]]);
942: }
943: }
944: PetscViewerASCIIPrintf(viewer,"\n");
945: bf -= 1;
946: }
947: }
948: }
949: if (bf) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Remaining boundary faces %D",bf);
950: PetscBTDestroy(&bfaces);
951: PetscFree(fcells);
952: }
954: /* mark owned vertices */
955: vown = NULL;
956: if (pown) {
957: PetscBTCreate(vEnd-vStart,&vown);
958: for (p=cStart;p<cEnd;p++) {
959: PetscInt i,closureSize,*closure = NULL;
961: if (PetscUnlikely(!PetscBTLookup(pown,p-cStart))) continue;
962: DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
963: for (i=0;i<closureSize;i++) {
964: const PetscInt pp = closure[2*i];
966: if (pp >= vStart && pp < vEnd) {
967: PetscBTSet(vown,pp-vStart);
968: }
969: }
970: DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
971: }
972: }
974: /* vertex_parents (Non-conforming meshes) */
975: parentSection = NULL;
976: if (enable_ncmesh) {
977: DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
978: }
979: if (parentSection) {
980: PetscInt vp,gvp;
982: for (vp=0,p=vStart;p<vEnd;p++) {
983: DMLabel dlabel;
984: PetscInt parent,depth;
986: if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue;
987: DMPlexGetDepthLabel(dm,&dlabel);
988: DMLabelGetValue(dlabel,p,&depth);
989: DMPlexGetTreeParent(dm,p,&parent,NULL);
990: if (parent != p) vp++;
991: }
992: MPIU_Allreduce(&vp,&gvp,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
993: if (gvp) {
994: PetscInt maxsupp;
995: PetscBool *skip = NULL;
997: PetscViewerASCIIPrintf(viewer,"\nvertex_parents\n");
998: PetscViewerASCIIPrintf(viewer,"%D\n",vp);
999: DMPlexGetMaxSizes(dm,NULL,&maxsupp);
1000: PetscMalloc1(maxsupp,&skip);
1001: for (p=vStart;p<vEnd;p++) {
1002: DMLabel dlabel;
1003: PetscInt parent;
1005: if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue;
1006: DMPlexGetDepthLabel(dm,&dlabel);
1007: DMPlexGetTreeParent(dm,p,&parent,NULL);
1008: if (parent != p) {
1009: int vids[8] = { -1, -1, -1, -1, -1, -1, -1, -1 }; /* silent overzealous clang static analyzer */
1010: PetscInt i,nv,ssize,n,numChildren,depth = -1;
1011: const PetscInt *children;
1013: DMPlexGetConeSize(dm,parent,&ssize);
1014: switch (ssize) {
1015: case 2: /* edge */
1016: nv = 0;
1017: DMPlexGetPointMFEMVertexIDs_Internal(dm,parent,localized ? coordSection : NULL,&nv,vids);
1018: DMPlexInvertCell(dim,nv,vids);
1019: PetscViewerASCIIPrintf(viewer,"%D",p-vStart);
1020: for (i=0;i<nv;i++) {
1021: PetscViewerASCIIPrintf(viewer," %D",(PetscInt)vids[i]);
1022: }
1023: PetscViewerASCIIPrintf(viewer,"\n");
1024: vp--;
1025: break;
1026: case 4: /* face */
1027: DMPlexGetTreeChildren(dm,parent,&numChildren,&children);
1028: for (n=0;n<numChildren;n++) {
1029: DMLabelGetValue(dlabel,children[n],&depth);
1030: if (!depth) {
1031: const PetscInt *hvsupp,*hesupp,*cone;
1032: PetscInt hvsuppSize,hesuppSize,coneSize;
1033: PetscInt hv = children[n],he = -1,f;
1035: PetscMemzero(skip,maxsupp*sizeof(PetscBool));
1036: DMPlexGetSupportSize(dm,hv,&hvsuppSize);
1037: DMPlexGetSupport(dm,hv,&hvsupp);
1038: for (i=0;i<hvsuppSize;i++) {
1039: PetscInt ep;
1040: DMPlexGetTreeParent(dm,hvsupp[i],&ep,NULL);
1041: if (ep != hvsupp[i]) {
1042: he = hvsupp[i];
1043: } else {
1044: skip[i] = PETSC_TRUE;
1045: }
1046: }
1047: if (he == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vertex %D support size %D: hanging edge not found",hv,hvsuppSize);
1048: DMPlexGetCone(dm,he,&cone);
1049: vids[0] = (int)((cone[0] == hv) ? cone[1] : cone[0]);
1050: DMPlexGetSupportSize(dm,he,&hesuppSize);
1051: DMPlexGetSupport(dm,he,&hesupp);
1052: for (f=0;f<hesuppSize;f++) {
1053: PetscInt j;
1055: DMPlexGetCone(dm,hesupp[f],&cone);
1056: DMPlexGetConeSize(dm,hesupp[f],&coneSize);
1057: for (j=0;j<coneSize;j++) {
1058: PetscInt k;
1059: for (k=0;k<hvsuppSize;k++) {
1060: if (hvsupp[k] == cone[j]) {
1061: skip[k] = PETSC_TRUE;
1062: break;
1063: }
1064: }
1065: }
1066: }
1067: for (i=0;i<hvsuppSize;i++) {
1068: if (!skip[i]) {
1069: DMPlexGetCone(dm,hvsupp[i],&cone);
1070: vids[1] = (int)((cone[0] == hv) ? cone[1] : cone[0]);
1071: }
1072: }
1073: PetscViewerASCIIPrintf(viewer,"%D",hv-vStart);
1074: for (i=0;i<2;i++) {
1075: PetscViewerASCIIPrintf(viewer," %D",(PetscInt)(vids[i]-vStart));
1076: }
1077: PetscViewerASCIIPrintf(viewer,"\n");
1078: vp--;
1079: }
1080: }
1081: break;
1082: default:
1083: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Don't know how to deal with support size %D",ssize);
1084: }
1085: }
1086: }
1087: PetscFree(skip);
1088: }
1089: if (vp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %D hanging vertices",vp);
1090: }
1091: PetscBTDestroy(&pown);
1092: PetscBTDestroy(&vown);
1094: /* vertices */
1095: if (hovec) { /* higher-order meshes */
1096: const char *fec;
1097: PetscInt i,n,s;
1099: PetscViewerASCIIPrintf(viewer,"\nvertices\n");
1100: PetscViewerASCIIPrintf(viewer,"%D\n",vEnd-vStart);
1101: PetscViewerASCIIPrintf(viewer,"nodes\n");
1102: PetscObjectGetName((PetscObject)hovec,&fec);
1103: PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");
1104: PetscViewerASCIIPrintf(viewer,"%s\n",fec);
1105: PetscViewerASCIIPrintf(viewer,"VDim: %D\n",sdim);
1106: PetscViewerASCIIPrintf(viewer,"Ordering: 1\n\n"); /*Ordering::byVDIM*/
1107: VecGetArrayRead(hovec,&array);
1108: VecGetLocalSize(hovec,&n);
1109: if (n%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Size of local coordinate vector %D incompatible with space dimension %D",n,sdim);
1110: for (i=0;i<n/sdim;i++) {
1111: for (s=0;s<sdim;s++) {
1112: PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[i*sdim+s]));
1113: }
1114: PetscViewerASCIIPrintf(viewer,"\n");
1115: }
1116: VecRestoreArrayRead(hovec,&array);
1117: } else {
1118: VecGetLocalSize(coordinates,&nvert);
1119: PetscViewerASCIIPrintf(viewer,"\nvertices\n");
1120: PetscViewerASCIIPrintf(viewer,"%D\n",nvert/sdim);
1121: PetscViewerASCIIPrintf(viewer,"%D\n",sdim);
1122: VecGetArrayRead(coordinates,&array);
1123: for (p=0;p<nvert/sdim;p++) {
1124: PetscInt s;
1125: for (s=0;s<sdim;s++) {
1126: PetscReal v = PetscRealPart(array[p*sdim+s]);
1128: PetscViewerASCIIPrintf(viewer,fmt,PetscIsInfOrNanReal(v) ? 0.0 : v);
1129: }
1130: PetscViewerASCIIPrintf(viewer,"\n");
1131: }
1132: VecRestoreArrayRead(coordinates,&array);
1133: }
1134: return(0);
1135: }
1137: PetscErrorCode DMPlexView_GLVis(DM dm, PetscViewer viewer)
1138: {
1141: DMView_GLVis(dm,viewer,DMPlexView_GLVis_ASCII);
1142: return(0);
1143: }