Actual source code: plexglvis.c
petsc-3.10.5 2019-03-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,cEndInterior,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: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
65: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
66: DMPlexGetCellNumbering(dm,&globalNum);
67: ISGetIndices(globalNum,&gNum);
68: PetscBTCreate(vEnd-vStart,&vown);
69: for (c = cStart, totc = 0; c < cEnd; c++) {
70: if (gNum[c-cStart] >= 0) {
71: PetscInt i,numPoints,*points = NULL;
73: totc++;
74: DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&numPoints,&points);
75: for (i=0;i<numPoints*2;i+= 2) {
76: if ((points[i] >= vStart) && (points[i] < vEnd)) {
77: PetscBTSet(vown,points[i]-vStart);
78: }
79: }
80: DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&numPoints,&points);
81: }
82: }
83: for (f=0,Nv=0;f<vEnd-vStart;f++) if (PetscLikely(PetscBTLookup(vown,f))) Nv++;
85: DMCreateLocalVector(dm,&xlocal);
86: VecGetLocalSize(xlocal,&totdofs);
87: DMGetSection(dm,&s);
88: PetscSectionGetNumFields(s,&nfields);
89: for (f=0,maxfields=0;f<nfields;f++) {
90: PetscInt bs;
92: PetscSectionGetFieldComponents(s,f,&bs);
93: maxfields += bs;
94: }
95: PetscCalloc7(maxfields,&fieldname,maxfields,&nlocal,maxfields,&bs,maxfields,&dims,maxfields,&fec_type,totdofs,&idxs,maxfields,&Ufield);
96: PetscNew(&ctx);
97: PetscCalloc1(maxfields,&ctx->scctx);
98: DMGetDS(dm,&ds);
99: if (ds) {
100: for (f=0;f<nfields;f++) {
101: const char* fname;
102: char name[256];
103: PetscObject disc;
104: size_t len;
106: PetscSectionGetFieldName(s,f,&fname);
107: PetscStrlen(fname,&len);
108: if (len) {
109: PetscStrcpy(name,fname);
110: } else {
111: PetscSNPrintf(name,256,"Field%D",f);
112: }
113: PetscDSGetDiscretization(ds,f,&disc);
114: if (disc) {
115: PetscClassId id;
116: PetscInt Nc;
117: char fec[64];
119: PetscObjectGetClassId(disc, &id);
120: if (id == PETSCFE_CLASSID) {
121: PetscFE fem = (PetscFE)disc;
122: PetscDualSpace sp;
123: PetscDualSpaceType spname;
124: PetscInt order;
125: PetscBool islag,continuous,H1 = PETSC_TRUE;
127: PetscFEGetNumComponents(fem,&Nc);
128: PetscFEGetDualSpace(fem,&sp);
129: PetscDualSpaceGetType(sp,&spname);
130: PetscStrcmp(spname,PETSCDUALSPACELAGRANGE,&islag);
131: if (!islag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported dual space");
132: PetscDualSpaceLagrangeGetContinuity(sp,&continuous);
133: if (!continuous) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Discontinuous space visualization currently unsupported");
134: PetscDualSpaceGetOrder(sp,&order);
135: if (continuous && order > 0) {
136: PetscSNPrintf(fec,64,"FiniteElementCollection: H1_%DD_P1",dim);
137: } else {
138: H1 = PETSC_FALSE;
139: PetscSNPrintf(fec,64,"FiniteElementCollection: L2_%DD_P%D",dim,order);
140: }
141: PetscStrallocpy(name,&fieldname[ctx->nf]);
142: bs[ctx->nf] = Nc;
143: dims[ctx->nf] = dim;
144: if (H1) {
145: nlocal[ctx->nf] = Nc * Nv;
146: PetscStrallocpy(fec,&fec_type[ctx->nf]);
147: VecCreateSeq(PETSC_COMM_SELF,Nv*Nc,&xfield);
148: for (i=0,cum=0;i<vEnd-vStart;i++) {
149: PetscInt j,off;
151: if (PetscUnlikely(!PetscBTLookup(vown,i))) continue;
152: PetscSectionGetFieldOffset(s,i+vStart,f,&off);
153: for (j=0;j<Nc;j++) idxs[cum++] = off + j;
154: }
155: ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),Nv*Nc,idxs,PETSC_USE_POINTER,&isfield);
156: } else {
157: nlocal[ctx->nf] = Nc * totc;
158: PetscStrallocpy(fec,&fec_type[ctx->nf]);
159: VecCreateSeq(PETSC_COMM_SELF,Nc*totc,&xfield);
160: for (i=0,cum=0;i<cEnd-cStart;i++) {
161: PetscInt j,off;
163: if (PetscUnlikely(gNum[i] < 0)) continue;
164: PetscSectionGetFieldOffset(s,i+cStart,f,&off);
165: for (j=0;j<Nc;j++) idxs[cum++] = off + j;
166: }
167: ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),totc*Nc,idxs,PETSC_USE_POINTER,&isfield);
168: }
169: VecScatterCreate(xlocal,isfield,xfield,NULL,&ctx->scctx[ctx->nf]);
170: VecDestroy(&xfield);
171: ISDestroy(&isfield);
172: ctx->nf++;
173: } else if (id == PETSCFV_CLASSID) {
174: PetscInt c;
176: PetscFVGetNumComponents((PetscFV)disc,&Nc);
177: PetscSNPrintf(fec,64,"FiniteElementCollection: L2_%DD_P0",dim);
178: for (c = 0; c < Nc; c++) {
179: char comp[256];
180: PetscSNPrintf(comp,256,"%s-Comp%D",name,c);
181: PetscStrallocpy(comp,&fieldname[ctx->nf]);
182: bs[ctx->nf] = 1; /* Does PetscFV support components with different block size? */
183: nlocal[ctx->nf] = totc;
184: dims[ctx->nf] = dim;
185: PetscStrallocpy(fec,&fec_type[ctx->nf]);
186: VecCreateSeq(PETSC_COMM_SELF,totc,&xfield);
187: for (i=0,cum=0;i<cEnd-cStart;i++) {
188: PetscInt off;
190: if (PetscUnlikely(gNum[i])<0) continue;
191: PetscSectionGetFieldOffset(s,i+cStart,f,&off);
192: idxs[cum++] = off + c;
193: }
194: ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),totc,idxs,PETSC_USE_POINTER,&isfield);
195: VecScatterCreate(xlocal,isfield,xfield,NULL,&ctx->scctx[ctx->nf]);
196: VecDestroy(&xfield);
197: ISDestroy(&isfield);
198: ctx->nf++;
199: }
200: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %D",f);
201: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing discretization for field %D",f);
202: }
203: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Needs a DS attached to the DM");
204: PetscBTDestroy(&vown);
205: VecDestroy(&xlocal);
206: ISRestoreIndices(globalNum,&gNum);
208: /* create work vectors */
209: for (f=0;f<ctx->nf;f++) {
210: VecCreateMPI(PetscObjectComm((PetscObject)dm),nlocal[f],PETSC_DECIDE,&Ufield[f]);
211: PetscObjectSetName((PetscObject)Ufield[f],fieldname[f]);
212: VecSetBlockSize(Ufield[f],bs[f]);
213: VecSetDM(Ufield[f],dm);
214: }
216: /* customize the viewer */
217: PetscViewerGLVisSetFields(viewer,ctx->nf,(const char**)fec_type,dims,DMPlexSampleGLVisFields_Private,(PetscObject*)Ufield,ctx,DestroyGLVisViewerCtx_Private);
218: for (f=0;f<ctx->nf;f++) {
219: PetscFree(fieldname[f]);
220: PetscFree(fec_type[f]);
221: VecDestroy(&Ufield[f]);
222: }
223: PetscFree7(fieldname,nlocal,bs,dims,fec_type,idxs,Ufield);
224: return(0);
225: }
227: typedef enum {MFEM_POINT=0,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE,MFEM_TETRAHEDRON,MFEM_CUBE,MFEM_UNDEF} MFEM_cid;
229: MFEM_cid mfem_table_cid[4][7] = { {MFEM_POINT,MFEM_UNDEF,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_UNDEF},
230: {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_UNDEF},
231: {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE ,MFEM_UNDEF,MFEM_UNDEF},
232: {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_TETRAHEDRON,MFEM_UNDEF,MFEM_CUBE } };
234: MFEM_cid mfem_table_cid_unint[4][9] = { {MFEM_POINT,MFEM_UNDEF,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF},
235: {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF},
236: {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE ,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF},
237: {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_TETRAHEDRON,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF,MFEM_CUBE } };
239: static PetscErrorCode DMPlexGetPointMFEMCellID_Internal(DM dm, DMLabel label, PetscInt p, PetscInt *mid, PetscInt *cid)
240: {
241: DMLabel dlabel;
242: PetscInt depth,csize,pdepth,dim;
246: DMPlexGetDepthLabel(dm,&dlabel);
247: DMLabelGetValue(dlabel,p,&pdepth);
248: DMPlexGetConeSize(dm,p,&csize);
249: DMPlexGetDepth(dm,&depth);
250: DMGetDimension(dm,&dim);
251: if (label) {
252: DMLabelGetValue(label,p,mid);
253: } else *mid = 1;
254: if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
255: #if defined PETSC_USE_DEBUG
256: if (dim < 0 || dim > 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %D",dim);
257: if (csize > 8) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Found cone size %D for point %D",csize,p);
258: if (depth != 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Found depth %D for point %D. You should interpolate the mesh first",depth,p);
259: #endif
260: *cid = mfem_table_cid_unint[dim][csize];
261: } else {
262: #if defined PETSC_USE_DEBUG
263: if (csize > 6) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cone size %D for point %D",csize,p);
264: if (pdepth < 0 || pdepth > 3) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Depth %D for point %D",csize,p);
265: #endif
266: *cid = mfem_table_cid[pdepth][csize];
267: }
268: return(0);
269: }
271: static PetscErrorCode DMPlexGetPointMFEMVertexIDs_Internal(DM dm, PetscInt p, PetscSection csec, PetscInt *nv, int vids[])
272: {
273: PetscInt dim,sdim,dof = 0,off = 0,i,q,vStart,vEnd,numPoints,*points = NULL;
277: DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);
278: DMGetDimension(dm,&dim);
279: sdim = dim;
280: if (csec) {
281: PetscInt sStart,sEnd;
283: DMGetCoordinateDim(dm,&sdim);
284: PetscSectionGetChart(csec,&sStart,&sEnd);
285: PetscSectionGetOffset(csec,vStart,&off);
286: off = off/sdim;
287: if (p >= sStart && p < sEnd) {
288: PetscSectionGetDof(csec,p,&dof);
289: }
290: }
291: if (!dof) {
292: DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);
293: for (i=0,q=0;i<numPoints*2;i+= 2)
294: if ((points[i] >= vStart) && (points[i] < vEnd))
295: vids[q++] = (int)(points[i]-vStart+off);
296: DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);
297: } else {
298: PetscSectionGetOffset(csec,p,&off);
299: PetscSectionGetDof(csec,p,&dof);
300: for (q=0;q<dof/sdim;q++) vids[q] = (int)(off/sdim + q);
301: }
302: *nv = q;
303: return(0);
304: }
306: /*
307: ASCII visualization/dump: full support for simplices and tensor product cells. It supports AMR
308: Higher order meshes are also supported
309: */
310: static PetscErrorCode DMPlexView_GLVis_ASCII(DM dm, PetscViewer viewer)
311: {
312: DMLabel label;
313: PetscSection coordSection,parentSection;
314: Vec coordinates,hovec;
315: const PetscScalar *array;
316: PetscInt bf,p,sdim,dim,depth,novl;
317: PetscInt cStart,cEnd,cEndInterior,vStart,vEnd,nvert;
318: PetscMPIInt size;
319: PetscBool localized,isascii;
320: PetscBool enable_mfem,enable_boundary,enable_ncmesh;
321: PetscBT pown,vown;
322: PetscErrorCode ierr;
323: PetscContainer glvis_container;
324: PetscBool cellvertex = PETSC_FALSE, periodic, enabled = PETSC_TRUE;
325: const char *fmt;
326: char emark[64] = "",bmark[64] = "";
331: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
332: if (!isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERASCII");
333: MPI_Comm_size(PetscObjectComm((PetscObject)viewer),&size);
334: if (size > 1) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Use single sequential viewers for parallel visualization");
335: DMGetDimension(dm,&dim);
337: /* get container: determines if a process visualizes is portion of the data or not */
338: PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);
339: if (!glvis_container) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing GLVis container");
340: {
341: PetscViewerGLVisInfo glvis_info;
342: PetscContainerGetPointer(glvis_container,(void**)&glvis_info);
343: enabled = glvis_info->enabled;
344: fmt = glvis_info->fmt;
345: }
347: /* Users can attach a coordinate vector to the DM in case they have a higher-order mesh
348: DMPlex does not currently support HO meshes, so there's no API for this */
349: PetscObjectQuery((PetscObject)dm,"_glvis_mesh_coords",(PetscObject*)&hovec);
351: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
352: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
353: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
354: DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);
355: DMGetPeriodicity(dm,&periodic,NULL,NULL,NULL);
356: DMGetCoordinatesLocalized(dm,&localized);
357: if (periodic && !localized && !hovec) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Coordinates need to be localized");
358: DMGetCoordinateSection(dm,&coordSection);
359: DMGetCoordinateDim(dm,&sdim);
360: DMGetCoordinatesLocal(dm,&coordinates);
361: if (!coordinates && !hovec) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
363: /*
364: a couple of sections of the mesh specification are disabled
365: - boundary: unless we want to visualize boundary attributes or we have a periodic mesh
366: the boundary is not needed for proper mesh visualization
367: - vertex_parents: used for non-conforming meshes only when we want to use MFEM as a discretization package
368: and be able to derefine the mesh (MFEM does not currently have to ability to read ncmeshes in parallel)
369: */
370: enable_boundary = periodic;
371: enable_ncmesh = PETSC_FALSE;
372: enable_mfem = PETSC_FALSE;
373: /* I'm tired of problems with negative values in the markers, disable them */
374: PetscOptionsBegin(PetscObjectComm((PetscObject)dm),((PetscObject)dm)->prefix,"GLVis PetscViewer DMPlex Options","PetscViewer");
375: PetscOptionsBool("-viewer_glvis_dm_plex_enable_boundary","Enable boundary section in mesh representation",NULL,enable_boundary,&enable_boundary,NULL);
376: PetscOptionsBool("-viewer_glvis_dm_plex_enable_ncmesh","Enable vertex_parents section in mesh representation (allows derefinement)",NULL,enable_ncmesh,&enable_ncmesh,NULL);
377: PetscOptionsBool("-viewer_glvis_dm_plex_enable_mfem","Dump a mesh that can be used with MFEM's FiniteElementSpaces",NULL,enable_mfem,&enable_mfem,NULL);
378: PetscOptionsString("-viewer_glvis_dm_plex_emarker","String for the material id label",NULL,emark,emark,sizeof(emark),NULL);
379: PetscOptionsString("-viewer_glvis_dm_plex_bmarker","String for the boundary id label",NULL,bmark,bmark,sizeof(bmark),NULL);
380: PetscOptionsEnd();
381: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
382: if (enable_ncmesh && size > 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not supported in parallel");
383: DMPlexGetDepth(dm,&depth);
384: if (enable_boundary && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. "
385: "Alternatively, run with -viewer_glvis_dm_plex_enable_boundary 0");
386: if (enable_ncmesh && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. "
387: "Alternatively, run with -viewer_glvis_dm_plex_enable_ncmesh 0");
388: if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
389: if (depth != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported depth %D. You should interpolate the mesh first",depth);
390: cellvertex = PETSC_TRUE;
391: }
393: /* Identify possible cells in the overlap */
394: novl = 0;
395: pown = NULL;
396: if (size > 1) {
397: IS globalNum = NULL;
398: const PetscInt *gNum;
399: PetscBool ovl = PETSC_FALSE;
401: DMPlexGetCellNumbering(dm,&globalNum);
402: ISGetIndices(globalNum,&gNum);
403: for (p=cStart; p<cEnd; p++) {
404: if (gNum[p-cStart] < 0) {
405: ovl = PETSC_TRUE;
406: novl++;
407: }
408: }
409: if (ovl) {
410: /* it may happen that pown get not destroyed, if the user closes the window while this function is running.
411: TODO: garbage collector? attach pown to dm? */
412: PetscBTCreate(cEnd-cStart,&pown);
413: for (p=cStart; p<cEnd; p++) {
414: if (gNum[p-cStart] < 0) continue;
415: else {
416: PetscBTSet(pown,p-cStart);
417: }
418: }
419: }
420: ISRestoreIndices(globalNum,&gNum);
421: }
423: /* return if this process is disabled */
424: if (!enabled) {
425: PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");
426: PetscViewerASCIIPrintf(viewer,"\ndimension\n");
427: PetscViewerASCIIPrintf(viewer,"%D\n",dim);
428: PetscViewerASCIIPrintf(viewer,"\nelements\n");
429: PetscViewerASCIIPrintf(viewer,"%D\n",0);
430: PetscViewerASCIIPrintf(viewer,"\nboundary\n");
431: PetscViewerASCIIPrintf(viewer,"%D\n",0);
432: PetscViewerASCIIPrintf(viewer,"\nvertices\n");
433: PetscViewerASCIIPrintf(viewer,"%D\n",0);
434: PetscViewerASCIIPrintf(viewer,"%D\n",sdim);
435: PetscBTDestroy(&pown);
436: return(0);
437: }
439: if (enable_mfem) {
440: if (periodic && !hovec) { /* we need to generate a vector of L2 coordinates, as this is how MFEM handles periodic meshes */
441: PetscInt vpc = 0;
442: char fec[64];
443: int vids[8] = {0,1,2,3,4,5,6,7};
444: int hexv[8] = {0,1,3,2,4,5,7,6}, tetv[4] = {0,1,2,3};
445: int quadv[8] = {0,1,3,2}, triv[3] = {0,1,2};
446: int *dof = NULL;
447: PetscScalar *array,*ptr;
449: PetscSNPrintf(fec,sizeof(fec),"FiniteElementCollection: L2_T1_%DD_P1",dim);
450: if (cEnd-cStart) {
451: PetscInt fpc;
453: DMPlexGetConeSize(dm,cStart,&fpc);
454: switch(dim) {
455: case 1:
456: vpc = 2;
457: dof = hexv;
458: break;
459: case 2:
460: switch (fpc) {
461: case 3:
462: vpc = 3;
463: dof = triv;
464: break;
465: case 4:
466: vpc = 4;
467: dof = quadv;
468: break;
469: default:
470: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
471: break;
472: }
473: break;
474: case 3:
475: switch (fpc) {
476: case 4: /* TODO: still need to understand L2 ordering for tets */
477: vpc = 4;
478: dof = tetv;
479: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled tethraedral case");
480: break;
481: case 6:
482: if (cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: vertices per cell %D",fpc);
483: vpc = 8;
484: dof = hexv;
485: break;
486: case 8:
487: if (!cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
488: vpc = 8;
489: dof = hexv;
490: break;
491: default:
492: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
493: break;
494: }
495: break;
496: default:
497: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim");
498: break;
499: }
500: DMPlexInvertCell(dim,vpc,vids);
501: }
502: if (!dof) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing dofs");
503: VecCreateSeq(PETSC_COMM_SELF,(cEnd-cStart-novl)*vpc*sdim,&hovec);
504: PetscObjectCompose((PetscObject)dm,"_glvis_mesh_coords",(PetscObject)hovec);
505: PetscObjectDereference((PetscObject)hovec);
506: PetscObjectSetName((PetscObject)hovec,fec);
507: VecGetArray(hovec,&array);
508: ptr = array;
509: for (p=cStart;p<cEnd;p++) {
510: PetscInt csize,v,d;
511: PetscScalar *vals = NULL;
513: if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue;
514: DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);
515: 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);
516: for (v=0;v<vpc;v++) {
517: for (d=0;d<sdim;d++) {
518: ptr[sdim*dof[v]+d] = vals[sdim*vids[v]+d];
519: }
520: }
521: ptr += vpc*sdim;
522: DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);
523: }
524: VecRestoreArray(hovec,&array);
525: }
526: }
529: /* header */
530: PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");
532: /* topological dimension */
533: PetscViewerASCIIPrintf(viewer,"\ndimension\n");
534: PetscViewerASCIIPrintf(viewer,"%D\n",dim);
536: /* elements */
537: DMGetLabel(dm,emark,&label);
538: PetscViewerASCIIPrintf(viewer,"\nelements\n");
539: PetscViewerASCIIPrintf(viewer,"%D\n",cEnd-cStart-novl);
540: for (p=cStart;p<cEnd;p++) {
541: int vids[8];
542: PetscInt i,nv = 0,cid = -1,mid = 1;
544: if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue;
545: DMPlexGetPointMFEMCellID_Internal(dm,label,p,&mid,&cid);
546: DMPlexGetPointMFEMVertexIDs_Internal(dm,p,(localized && !hovec) ? coordSection : NULL,&nv,vids);
547: DMPlexInvertCell(dim,nv,vids);
548: PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);
549: for (i=0;i<nv;i++) {
550: PetscViewerASCIIPrintf(viewer," %D",(PetscInt)vids[i]);
551: }
552: PetscViewerASCIIPrintf(viewer,"\n");
553: }
555: /* boundary */
556: PetscViewerASCIIPrintf(viewer,"\nboundary\n");
557: if (!enable_boundary) {
558: PetscViewerASCIIPrintf(viewer,"%D\n",0);
559: } else {
560: DMLabel perLabel;
561: PetscBT bfaces;
562: PetscInt fStart,fEnd,fEndInterior,*fcells;
563: PetscInt *faces = NULL,fpc = 0,vpf = 0, vpc = 0;
564: PetscInt fv1[] = {0,1},
565: fv2tri[] = {0,1,
566: 1,2,
567: 2,0},
568: fv2quad[] = {0,1,
569: 1,2,
570: 2,3,
571: 3,0},
572: fv3tet[] = {0,1,2,
573: 0,3,1,
574: 0,2,3,
575: 2,1,3},
576: fv3hex[] = {0,1,2,3,
577: 4,5,6,7,
578: 0,3,5,4,
579: 2,1,7,6,
580: 3,2,6,5,
581: 0,4,7,1};
583: /* determine orientation of boundary mesh */
584: if (cEnd-cStart) {
585: DMPlexGetConeSize(dm,cStart,&fpc);
586: switch(dim) {
587: case 1:
588: if (fpc != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case faces per cell %D",fpc);
589: faces = fv1;
590: vpf = 1;
591: vpc = 2;
592: break;
593: case 2:
594: switch (fpc) {
595: case 3:
596: faces = fv2tri;
597: vpf = 2;
598: vpc = 3;
599: break;
600: case 4:
601: faces = fv2quad;
602: vpf = 2;
603: vpc = 4;
604: break;
605: default:
606: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
607: break;
608: }
609: break;
610: case 3:
611: switch (fpc) {
612: case 4:
613: faces = fv3tet;
614: vpf = 3;
615: vpc = 4;
616: break;
617: case 6:
618: faces = fv3hex;
619: vpf = 4;
620: vpc = 8;
621: break;
622: default:
623: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
624: break;
625: }
626: break;
627: default:
628: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim");
629: break;
630: }
631: }
632: DMPlexGetHybridBounds(dm,NULL,&fEndInterior,NULL,NULL);
633: DMPlexGetHeightStratum(dm,1,&fStart,&fEnd);
634: fEnd = fEndInterior < 0 ? fEnd : fEndInterior;
635: PetscBTCreate(fEnd-fStart,&bfaces);
636: DMPlexGetMaxSizes(dm,NULL,&p);
637: PetscMalloc1(p,&fcells);
638: DMGetLabel(dm,"glvis_periodic_cut",&perLabel);
639: if (!perLabel && localized) { /* this periodic cut can be moved up to DMPlex setup */
640: DMCreateLabel(dm,"glvis_periodic_cut");
641: DMGetLabel(dm,"glvis_periodic_cut",&perLabel);
642: DMLabelSetDefaultValue(perLabel,1);
643: for (p=cStart;p<cEnd;p++) {
644: PetscInt dof;
645: PetscSectionGetDof(coordSection,p,&dof);
646: if (dof) {
647: PetscInt v,csize,cellClosureSize,*cellClosure = NULL,*vidxs = NULL;
648: PetscScalar *vals = NULL;
650: if (dof%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Incompatible number of cell dofs %D and space dimension %D",dof,sdim);
651: if (dof/sdim != vpc) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Incompatible number of cell dofs %D, vertices %D and space dim %D",dof/sdim,vpc,sdim);
652: DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);
653: DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);
654: for (v=0;v<cellClosureSize;v++)
655: if (cellClosure[2*v] >= vStart && cellClosure[2*v] < vEnd) {
656: vidxs = cellClosure + 2*v;
657: break;
658: }
659: if (!vidxs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing vertices");
660: for (v=0;v<vpc;v++) {
661: PetscInt s;
663: for (s=0;s<sdim;s++) {
664: if (PetscAbsScalar(vals[v*sdim+s]-vals[v*sdim+s+vpc*sdim])>PETSC_MACHINE_EPSILON) {
665: DMLabelSetValue(perLabel,vidxs[2*v],2);
666: }
667: }
668: }
669: DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);
670: DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);
671: }
672: }
673: if (dim > 1) {
674: PetscInt eEnd,eStart,eEndInterior;
676: DMPlexGetHybridBounds(dm,NULL,NULL,&eEndInterior,NULL);
677: DMPlexGetDepthStratum(dm,1,&eStart,&eEnd);
678: eEnd = eEndInterior < 0 ? eEnd : eEndInterior;
679: for (p=eStart;p<eEnd;p++) {
680: const PetscInt *cone;
681: PetscInt coneSize,i;
682: PetscBool ispe = PETSC_TRUE;
684: DMPlexGetCone(dm,p,&cone);
685: DMPlexGetConeSize(dm,p,&coneSize);
686: for (i=0;i<coneSize;i++) {
687: PetscInt v;
689: DMLabelGetValue(perLabel,cone[i],&v);
690: ispe = (PetscBool)(ispe && (v==2));
691: }
692: if (ispe && coneSize) {
693: DMLabelSetValue(perLabel,p,2);
694: }
695: }
696: if (dim > 2) {
697: for (p=fStart;p<fEnd;p++) {
698: const PetscInt *cone;
699: PetscInt coneSize,i;
700: PetscBool ispe = PETSC_TRUE;
702: DMPlexGetCone(dm,p,&cone);
703: DMPlexGetConeSize(dm,p,&coneSize);
704: for (i=0;i<coneSize;i++) {
705: PetscInt v;
707: DMLabelGetValue(perLabel,cone[i],&v);
708: ispe = (PetscBool)(ispe && (v==2));
709: }
710: if (ispe && coneSize) {
711: DMLabelSetValue(perLabel,p,2);
712: }
713: }
714: }
715: }
716: }
717: for (p=fStart;p<fEnd;p++) {
718: const PetscInt *support;
719: PetscInt supportSize;
720: PetscBool isbf = PETSC_FALSE;
722: DMPlexGetSupportSize(dm,p,&supportSize);
723: if (pown) {
724: PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE;
725: PetscInt i;
727: DMPlexGetSupport(dm,p,&support);
728: for (i=0;i<supportSize;i++) {
729: if (PetscLikely(PetscBTLookup(pown,support[i]-cStart))) has_owned = PETSC_TRUE;
730: else has_ghost = PETSC_TRUE;
731: }
732: isbf = (PetscBool)((supportSize == 1 && has_owned) || (supportSize > 1 && has_owned && has_ghost));
733: } else {
734: isbf = (PetscBool)(supportSize == 1);
735: }
736: if (!isbf && perLabel) {
737: const PetscInt *cone;
738: PetscInt coneSize,i;
740: DMPlexGetCone(dm,p,&cone);
741: DMPlexGetConeSize(dm,p,&coneSize);
742: isbf = PETSC_TRUE;
743: for (i=0;i<coneSize;i++) {
744: PetscInt v,d;
746: DMLabelGetValue(perLabel,cone[i],&v);
747: DMLabelGetDefaultValue(perLabel,&d);
748: isbf = (PetscBool)(isbf && v != d);
749: }
750: }
751: if (isbf) {
752: PetscBTSet(bfaces,p-fStart);
753: }
754: }
755: /* count boundary faces */
756: for (p=fStart,bf=0;p<fEnd;p++) {
757: if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) {
758: const PetscInt *support;
759: PetscInt supportSize,c;
761: DMPlexGetSupportSize(dm,p,&supportSize);
762: DMPlexGetSupport(dm,p,&support);
763: if (pown) {
764: for (c=0;c<supportSize;c++) {
765: if (PetscLikely(PetscBTLookup(pown,support[c]-cStart))) {
766: bf++;
767: }
768: }
769: } else bf += supportSize;
770: }
771: }
772: DMGetLabel(dm,bmark,&label);
773: PetscViewerASCIIPrintf(viewer,"%D\n",bf);
774: for (p=fStart;p<fEnd;p++) {
775: if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) {
776: const PetscInt *support;
777: PetscInt supportSize,c,nc = 0;
779: DMPlexGetSupportSize(dm,p,&supportSize);
780: DMPlexGetSupport(dm,p,&support);
781: if (pown) {
782: for (c=0;c<supportSize;c++) {
783: if (PetscLikely(PetscBTLookup(pown,support[c]-cStart))) {
784: fcells[nc++] = support[c];
785: }
786: }
787: } else for (c=0;c<supportSize;c++) fcells[nc++] = support[c];
788: for (c=0;c<nc;c++) {
789: const PetscInt *cone;
790: int vids[8];
791: PetscInt i,cell,cl,nv,cid = -1,mid = -1;
793: cell = fcells[c];
794: DMPlexGetCone(dm,cell,&cone);
795: for (cl=0;cl<fpc;cl++)
796: if (cone[cl] == p)
797: break;
799: /* face material id and type */
800: DMPlexGetPointMFEMCellID_Internal(dm,label,p,&mid,&cid);
801: PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);
802: /* vertex ids */
803: DMPlexGetPointMFEMVertexIDs_Internal(dm,cell,(localized && !hovec) ? coordSection : NULL,&nv,vids);
804: for (i=0;i<vpf;i++) {
805: PetscViewerASCIIPrintf(viewer," %d",vids[faces[cl*vpf+i]]);
806: }
807: PetscViewerASCIIPrintf(viewer,"\n");
808: }
809: bf = bf-nc;
810: }
811: }
812: if (bf) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Remaining boundary faces %D",bf);
813: PetscBTDestroy(&bfaces);
814: PetscFree(fcells);
815: }
817: /* mark owned vertices */
818: vown = NULL;
819: if (pown) {
820: PetscBTCreate(vEnd-vStart,&vown);
821: for (p=cStart;p<cEnd;p++) {
822: PetscInt i,closureSize,*closure = NULL;
824: if (PetscUnlikely(!PetscBTLookup(pown,p-cStart))) continue;
825: DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
826: for (i=0;i<closureSize;i++) {
827: const PetscInt pp = closure[2*i];
829: if (pp >= vStart && pp < vEnd) {
830: PetscBTSet(vown,pp-vStart);
831: }
832: }
833: DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
834: }
835: }
837: /* vertex_parents (Non-conforming meshes) */
838: parentSection = NULL;
839: if (enable_ncmesh) {
840: DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
841: }
842: if (parentSection) {
843: PetscInt vp,gvp;
845: for (vp=0,p=vStart;p<vEnd;p++) {
846: DMLabel dlabel;
847: PetscInt parent,depth;
849: if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue;
850: DMPlexGetDepthLabel(dm,&dlabel);
851: DMLabelGetValue(dlabel,p,&depth);
852: DMPlexGetTreeParent(dm,p,&parent,NULL);
853: if (parent != p) vp++;
854: }
855: MPIU_Allreduce(&vp,&gvp,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
856: if (gvp) {
857: PetscInt maxsupp;
858: PetscBool *skip = NULL;
860: PetscViewerASCIIPrintf(viewer,"\nvertex_parents\n");
861: PetscViewerASCIIPrintf(viewer,"%D\n",vp);
862: DMPlexGetMaxSizes(dm,NULL,&maxsupp);
863: PetscMalloc1(maxsupp,&skip);
864: for (p=vStart;p<vEnd;p++) {
865: DMLabel dlabel;
866: PetscInt parent;
868: if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue;
869: DMPlexGetDepthLabel(dm,&dlabel);
870: DMPlexGetTreeParent(dm,p,&parent,NULL);
871: if (parent != p) {
872: int vids[8] = { -1, -1, -1, -1, -1, -1, -1, -1 }; /* silent overzealous clang static analyzer */
873: PetscInt i,nv,ssize,n,numChildren,depth = -1;
874: const PetscInt *children;
876: DMPlexGetConeSize(dm,parent,&ssize);
877: switch (ssize) {
878: case 2: /* edge */
879: nv = 0;
880: DMPlexGetPointMFEMVertexIDs_Internal(dm,parent,localized ? coordSection : NULL,&nv,vids);
881: DMPlexInvertCell(dim,nv,vids);
882: PetscViewerASCIIPrintf(viewer,"%D",p-vStart);
883: for (i=0;i<nv;i++) {
884: PetscViewerASCIIPrintf(viewer," %D",(PetscInt)vids[i]);
885: }
886: PetscViewerASCIIPrintf(viewer,"\n");
887: vp--;
888: break;
889: case 4: /* face */
890: DMPlexGetTreeChildren(dm,parent,&numChildren,&children);
891: for (n=0;n<numChildren;n++) {
892: DMLabelGetValue(dlabel,children[n],&depth);
893: if (!depth) {
894: const PetscInt *hvsupp,*hesupp,*cone;
895: PetscInt hvsuppSize,hesuppSize,coneSize;
896: PetscInt hv = children[n],he = -1,f;
898: PetscMemzero(skip,maxsupp*sizeof(PetscBool));
899: DMPlexGetSupportSize(dm,hv,&hvsuppSize);
900: DMPlexGetSupport(dm,hv,&hvsupp);
901: for (i=0;i<hvsuppSize;i++) {
902: PetscInt ep;
903: DMPlexGetTreeParent(dm,hvsupp[i],&ep,NULL);
904: if (ep != hvsupp[i]) {
905: he = hvsupp[i];
906: } else {
907: skip[i] = PETSC_TRUE;
908: }
909: }
910: if (he == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vertex %D support size %D: hanging edge not found",hv,hvsuppSize);
911: DMPlexGetCone(dm,he,&cone);
912: vids[0] = (int)((cone[0] == hv) ? cone[1] : cone[0]);
913: DMPlexGetSupportSize(dm,he,&hesuppSize);
914: DMPlexGetSupport(dm,he,&hesupp);
915: for (f=0;f<hesuppSize;f++) {
916: PetscInt j;
918: DMPlexGetCone(dm,hesupp[f],&cone);
919: DMPlexGetConeSize(dm,hesupp[f],&coneSize);
920: for (j=0;j<coneSize;j++) {
921: PetscInt k;
922: for (k=0;k<hvsuppSize;k++) {
923: if (hvsupp[k] == cone[j]) {
924: skip[k] = PETSC_TRUE;
925: break;
926: }
927: }
928: }
929: }
930: for (i=0;i<hvsuppSize;i++) {
931: if (!skip[i]) {
932: DMPlexGetCone(dm,hvsupp[i],&cone);
933: vids[1] = (int)((cone[0] == hv) ? cone[1] : cone[0]);
934: }
935: }
936: PetscViewerASCIIPrintf(viewer,"%D",hv-vStart);
937: for (i=0;i<2;i++) {
938: PetscViewerASCIIPrintf(viewer," %D",(PetscInt)(vids[i]-vStart));
939: }
940: PetscViewerASCIIPrintf(viewer,"\n");
941: vp--;
942: }
943: }
944: break;
945: default:
946: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Don't know how to deal with support size %D",ssize);
947: }
948: }
949: }
950: PetscFree(skip);
951: }
952: if (vp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %D hanging vertices",vp);
953: }
954: PetscBTDestroy(&pown);
955: PetscBTDestroy(&vown);
957: /* vertices */
958: if (hovec) { /* higher-order meshes */
959: const char *fec;
960: PetscInt i,n;
962: PetscViewerASCIIPrintf(viewer,"\nvertices\n");
963: PetscViewerASCIIPrintf(viewer,"%D\n",vEnd-vStart);
964: PetscViewerASCIIPrintf(viewer,"nodes\n");
965: PetscObjectGetName((PetscObject)hovec,&fec);
966: PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");
967: PetscViewerASCIIPrintf(viewer,"%s\n",fec);
968: PetscViewerASCIIPrintf(viewer,"VDim: %D\n",sdim);
969: PetscViewerASCIIPrintf(viewer,"Ordering: 1\n\n"); /*Ordering::byVDIM*/
970: VecGetArrayRead(hovec,&array);
971: VecGetLocalSize(hovec,&n);
972: if (n%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Size of local coordinate vector %D incompatible with space dimension %D",n,sdim);
973: for (i=0;i<n/sdim;i++) {
974: PetscInt s;
975: for (s=0;s<sdim;s++) {
976: PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[i*sdim+s]));
977: }
978: PetscViewerASCIIPrintf(viewer,"\n");
979: }
980: VecRestoreArrayRead(hovec,&array);
981: } else {
982: VecGetLocalSize(coordinates,&nvert);
983: PetscViewerASCIIPrintf(viewer,"\nvertices\n");
984: PetscViewerASCIIPrintf(viewer,"%D\n",nvert/sdim);
985: PetscViewerASCIIPrintf(viewer,"%D\n",sdim);
986: VecGetArrayRead(coordinates,&array);
987: for (p=0;p<nvert/sdim;p++) {
988: PetscInt s;
989: for (s=0;s<sdim;s++) {
990: PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[p*sdim+s]));
991: }
992: PetscViewerASCIIPrintf(viewer,"\n");
993: }
994: VecRestoreArrayRead(coordinates,&array);
995: }
996: return(0);
997: }
999: /* dispatching, prints through the socket by prepending the mesh keyword to the usual ASCII dump */
1000: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM dm, PetscViewer viewer)
1001: {
1003: PetscBool isglvis,isascii;
1008: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);
1009: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1010: if (!isglvis && !isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERGLVIS or VIEWERASCII");
1011: if (isglvis) {
1012: PetscViewer view;
1013: PetscViewerGLVisType type;
1015: PetscViewerGLVisGetType_Private(viewer,&type);
1016: PetscViewerGLVisGetDMWindow_Private(viewer,&view);
1017: if (view) { /* in the socket case, it may happen that the connection failed */
1018: if (type == PETSC_VIEWER_GLVIS_SOCKET) {
1019: PetscMPIInt size,rank;
1021: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
1022: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1023: PetscViewerASCIIPrintf(view,"parallel %d %d\nmesh\n",size,rank);
1024: }
1025: DMPlexView_GLVis_ASCII(dm,view);
1026: PetscViewerFlush(view);
1027: if (type == PETSC_VIEWER_GLVIS_SOCKET) {
1028: PetscInt dim;
1029: const char* name;
1031: PetscObjectGetName((PetscObject)dm,&name);
1032: DMGetDimension(dm,&dim);
1033: PetscViewerGLVisInitWindow_Private(view,PETSC_TRUE,dim,name);
1034: PetscBarrier((PetscObject)dm);
1035: }
1036: }
1037: } else {
1038: DMPlexView_GLVis_ASCII(dm,viewer);
1039: }
1040: return(0);
1041: }