Actual source code: plexglvis.c
petsc-3.9.4 2018-09-11
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: DMGetDefaultSection(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: }
346: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
347: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
348: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
349: DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);
350: DMGetPeriodicity(dm,&periodic,NULL,NULL,NULL);
351: DMGetCoordinatesLocalized(dm,&localized);
352: if (periodic && !localized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Coordinates need to be localized");
353: DMGetCoordinateSection(dm,&coordSection);
354: DMGetCoordinateDim(dm,&sdim);
355: DMGetCoordinatesLocal(dm,&coordinates);
356: if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
358: /* Users can attach a coordinate vector to the DM in case they have a higher-order mesh
359: DMPlex does not currently support HO meshes, so there's no API for this */
360: PetscObjectQuery((PetscObject)dm,"_glvis_mesh_coords",(PetscObject*)&hovec);
362: /*
363: a couple of sections of the mesh specification are disabled
364: - boundary: unless we want to visualize boundary attributes or we have a periodic mesh
365: the boundary is not needed for proper mesh visualization
366: - vertex_parents: used for non-conforming meshes only when we want to use MFEM as a discretization package
367: and be able to derefine the mesh (MFEM does not currently have to ability to read ncmeshes in parallel)
368: */
369: enable_boundary = periodic;
370: enable_ncmesh = PETSC_FALSE;
371: enable_mfem = PETSC_FALSE;
372: /* I'm tired of problems with negative values in the markers, disable them */
373: PetscOptionsBegin(PetscObjectComm((PetscObject)dm),((PetscObject)dm)->prefix,"GLVis PetscViewer DMPlex Options","PetscViewer");
374: PetscOptionsBool("-viewer_glvis_dm_plex_enable_boundary","Enable boundary section in mesh representation",NULL,enable_boundary,&enable_boundary,NULL);
375: PetscOptionsBool("-viewer_glvis_dm_plex_enable_ncmesh","Enable vertex_parents section in mesh representation (allows derefinement)",NULL,enable_ncmesh,&enable_ncmesh,NULL);
376: PetscOptionsBool("-viewer_glvis_dm_plex_enable_mfem","Dump a mesh that can be used with MFEM's FiniteElementSpaces",NULL,enable_mfem,&enable_mfem,NULL);
377: PetscOptionsString("-viewer_glvis_dm_plex_emarker","String for the material id label",NULL,emark,emark,sizeof(emark),NULL);
378: PetscOptionsString("-viewer_glvis_dm_plex_bmarker","String for the boundary id label",NULL,bmark,bmark,sizeof(bmark),NULL);
379: PetscOptionsEnd();
380: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
381: if (enable_ncmesh && size > 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not supported in parallel");
382: DMPlexGetDepth(dm,&depth);
383: if (enable_boundary && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. "
384: "Alternatively, run with -viewer_glvis_dm_plex_enable_boundary 0");
385: if (enable_ncmesh && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. "
386: "Alternatively, run with -viewer_glvis_dm_plex_enable_ncmesh 0");
387: if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
388: if (depth != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported depth %D. You should interpolate the mesh first",depth);
389: cellvertex = PETSC_TRUE;
390: }
392: /* Identify possible cells in the overlap */
393: novl = 0;
394: pown = NULL;
395: if (size > 1) {
396: IS globalNum = NULL;
397: const PetscInt *gNum;
398: PetscBool ovl = PETSC_FALSE;
400: DMPlexGetCellNumbering(dm,&globalNum);
401: ISGetIndices(globalNum,&gNum);
402: for (p=cStart; p<cEnd; p++) {
403: if (gNum[p-cStart] < 0) {
404: ovl = PETSC_TRUE;
405: novl++;
406: }
407: }
408: if (ovl) {
409: /* it may happen that pown get not destroyed, if the user closes the window while this function is running.
410: TODO: garbage collector? attach pown to dm? */
411: PetscBTCreate(cEnd-cStart,&pown);
412: for (p=cStart; p<cEnd; p++) {
413: if (gNum[p-cStart] < 0) continue;
414: else {
415: PetscBTSet(pown,p-cStart);
416: }
417: }
418: }
419: ISRestoreIndices(globalNum,&gNum);
420: }
422: /* return if this process is disabled */
423: if (!enabled) {
424: PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");
425: PetscViewerASCIIPrintf(viewer,"\ndimension\n");
426: PetscViewerASCIIPrintf(viewer,"%D\n",dim);
427: PetscViewerASCIIPrintf(viewer,"\nelements\n");
428: PetscViewerASCIIPrintf(viewer,"%D\n",0);
429: PetscViewerASCIIPrintf(viewer,"\nboundary\n");
430: PetscViewerASCIIPrintf(viewer,"%D\n",0);
431: PetscViewerASCIIPrintf(viewer,"\nvertices\n");
432: PetscViewerASCIIPrintf(viewer,"%D\n",0);
433: PetscViewerASCIIPrintf(viewer,"%D\n",sdim);
434: PetscBTDestroy(&pown);
435: return(0);
436: }
438: if (enable_mfem) {
439: if (periodic && !hovec) { /* we need to generate a vector of L2 coordinates, as this is how MFEM handles periodic meshes */
440: PetscInt vpc = 0;
441: char fec[64];
442: int vids[8] = {0,1,2,3,4,5,6,7};
443: int hexv[8] = {0,1,3,2,4,5,7,6}, tetv[4] = {0,1,2,3};
444: int quadv[8] = {0,1,3,2}, triv[3] = {0,1,2};
445: int *dof = NULL;
446: PetscScalar *array,*ptr;
448: PetscSNPrintf(fec,sizeof(fec),"FiniteElementCollection: L2_T1_%DD_P1",dim);
449: if (cEnd-cStart) {
450: PetscInt fpc;
452: DMPlexGetConeSize(dm,cStart,&fpc);
453: switch(dim) {
454: case 1:
455: vpc = 2;
456: dof = hexv;
457: break;
458: case 2:
459: switch (fpc) {
460: case 3:
461: vpc = 3;
462: dof = triv;
463: break;
464: case 4:
465: vpc = 4;
466: dof = quadv;
467: break;
468: default:
469: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
470: break;
471: }
472: break;
473: case 3:
474: switch (fpc) {
475: case 4: /* TODO: still need to understand L2 ordering for tets */
476: vpc = 4;
477: dof = tetv;
478: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled tethraedral case");
479: break;
480: case 6:
481: if (cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: vertices per cell %D",fpc);
482: vpc = 8;
483: dof = hexv;
484: break;
485: case 8:
486: if (!cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
487: vpc = 8;
488: dof = hexv;
489: break;
490: default:
491: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
492: break;
493: }
494: break;
495: default:
496: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim");
497: break;
498: }
499: DMPlexInvertCell(dim,vpc,vids);
500: }
501: if (!dof) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing dofs");
502: VecCreateSeq(PETSC_COMM_SELF,(cEnd-cStart-novl)*vpc*sdim,&hovec);
503: PetscObjectCompose((PetscObject)dm,"_glvis_mesh_coords",(PetscObject)hovec);
504: PetscObjectDereference((PetscObject)hovec);
505: PetscObjectSetName((PetscObject)hovec,fec);
506: VecGetArray(hovec,&array);
507: ptr = array;
508: for (p=cStart;p<cEnd;p++) {
509: PetscInt csize,v,d;
510: PetscScalar *vals = NULL;
512: if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue;
513: DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);
514: 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);
515: for (v=0;v<vpc;v++) {
516: for (d=0;d<sdim;d++) {
517: ptr[sdim*dof[v]+d] = vals[sdim*vids[v]+d];
518: }
519: }
520: ptr += vpc*sdim;
521: DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);
522: }
523: VecRestoreArray(hovec,&array);
524: }
525: }
528: /* header */
529: PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");
531: /* topological dimension */
532: PetscViewerASCIIPrintf(viewer,"\ndimension\n");
533: PetscViewerASCIIPrintf(viewer,"%D\n",dim);
535: /* elements */
536: DMGetLabel(dm,emark,&label);
537: PetscViewerASCIIPrintf(viewer,"\nelements\n");
538: PetscViewerASCIIPrintf(viewer,"%D\n",cEnd-cStart-novl);
539: for (p=cStart;p<cEnd;p++) {
540: int vids[8];
541: PetscInt i,nv = 0,cid = -1,mid = 1;
543: if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue;
544: DMPlexGetPointMFEMCellID_Internal(dm,label,p,&mid,&cid);
545: DMPlexGetPointMFEMVertexIDs_Internal(dm,p,(localized && !hovec) ? coordSection : NULL,&nv,vids);
546: DMPlexInvertCell(dim,nv,vids);
547: PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);
548: for (i=0;i<nv;i++) {
549: PetscViewerASCIIPrintf(viewer," %D",(PetscInt)vids[i]);
550: }
551: PetscViewerASCIIPrintf(viewer,"\n");
552: }
554: /* boundary */
555: PetscViewerASCIIPrintf(viewer,"\nboundary\n");
556: if (!enable_boundary) {
557: PetscViewerASCIIPrintf(viewer,"%D\n",0);
558: } else {
559: DMLabel perLabel;
560: PetscBT bfaces;
561: PetscInt fStart,fEnd,fEndInterior,*fcells;
562: PetscInt *faces = NULL,fpc = 0,vpf = 0, vpc = 0;
563: PetscInt fv1[] = {0,1},
564: fv2tri[] = {0,1,
565: 1,2,
566: 2,0},
567: fv2quad[] = {0,1,
568: 1,2,
569: 2,3,
570: 3,0},
571: fv3tet[] = {0,1,2,
572: 0,3,1,
573: 0,2,3,
574: 2,1,3},
575: fv3hex[] = {0,1,2,3,
576: 4,5,6,7,
577: 0,3,5,4,
578: 2,1,7,6,
579: 3,2,6,5,
580: 0,4,7,1};
582: /* determine orientation of boundary mesh */
583: if (cEnd-cStart) {
584: DMPlexGetConeSize(dm,cStart,&fpc);
585: switch(dim) {
586: case 1:
587: if (fpc != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case faces per cell %D",fpc);
588: faces = fv1;
589: vpf = 1;
590: vpc = 2;
591: break;
592: case 2:
593: switch (fpc) {
594: case 3:
595: faces = fv2tri;
596: vpf = 2;
597: vpc = 3;
598: break;
599: case 4:
600: faces = fv2quad;
601: vpf = 2;
602: vpc = 4;
603: break;
604: default:
605: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
606: break;
607: }
608: break;
609: case 3:
610: switch (fpc) {
611: case 4:
612: faces = fv3tet;
613: vpf = 3;
614: vpc = 4;
615: break;
616: case 6:
617: faces = fv3hex;
618: vpf = 4;
619: vpc = 8;
620: break;
621: default:
622: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc);
623: break;
624: }
625: break;
626: default:
627: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim");
628: break;
629: }
630: }
631: DMPlexGetHybridBounds(dm,NULL,&fEndInterior,NULL,NULL);
632: DMPlexGetHeightStratum(dm,1,&fStart,&fEnd);
633: fEnd = fEndInterior < 0 ? fEnd : fEndInterior;
634: PetscBTCreate(fEnd-fStart,&bfaces);
635: DMPlexGetMaxSizes(dm,NULL,&p);
636: PetscMalloc1(p,&fcells);
637: DMGetLabel(dm,"glvis_periodic_cut",&perLabel);
638: if (!perLabel && localized) { /* this periodic cut can be moved up to DMPlex setup */
639: DMCreateLabel(dm,"glvis_periodic_cut");
640: DMGetLabel(dm,"glvis_periodic_cut",&perLabel);
641: DMLabelSetDefaultValue(perLabel,1);
642: for (p=cStart;p<cEnd;p++) {
643: PetscInt dof;
644: PetscSectionGetDof(coordSection,p,&dof);
645: if (dof) {
646: PetscInt v,csize,cellClosureSize,*cellClosure = NULL,*vidxs = NULL;
647: PetscScalar *vals = NULL;
649: if (dof%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Incompatible number of cell dofs %D and space dimension %D",dof,sdim);
650: 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);
651: DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);
652: DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);
653: for (v=0;v<cellClosureSize;v++)
654: if (cellClosure[2*v] >= vStart && cellClosure[2*v] < vEnd) {
655: vidxs = cellClosure + 2*v;
656: break;
657: }
658: if (!vidxs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing vertices");
659: for (v=0;v<vpc;v++) {
660: PetscInt s;
662: for (s=0;s<sdim;s++) {
663: if (PetscAbsScalar(vals[v*sdim+s]-vals[v*sdim+s+vpc*sdim])>PETSC_MACHINE_EPSILON) {
664: DMLabelSetValue(perLabel,vidxs[2*v],2);
665: }
666: }
667: }
668: DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);
669: DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);
670: }
671: }
672: if (dim > 1) {
673: PetscInt eEnd,eStart,eEndInterior;
675: DMPlexGetHybridBounds(dm,NULL,NULL,&eEndInterior,NULL);
676: DMPlexGetDepthStratum(dm,1,&eStart,&eEnd);
677: eEnd = eEndInterior < 0 ? eEnd : eEndInterior;
678: for (p=eStart;p<eEnd;p++) {
679: const PetscInt *cone;
680: PetscInt coneSize,i;
681: PetscBool ispe = PETSC_TRUE;
683: DMPlexGetCone(dm,p,&cone);
684: DMPlexGetConeSize(dm,p,&coneSize);
685: for (i=0;i<coneSize;i++) {
686: PetscInt v;
688: DMLabelGetValue(perLabel,cone[i],&v);
689: ispe = (PetscBool)(ispe && (v==2));
690: }
691: if (ispe && coneSize) {
692: DMLabelSetValue(perLabel,p,2);
693: }
694: }
695: if (dim > 2) {
696: for (p=fStart;p<fEnd;p++) {
697: const PetscInt *cone;
698: PetscInt coneSize,i;
699: PetscBool ispe = PETSC_TRUE;
701: DMPlexGetCone(dm,p,&cone);
702: DMPlexGetConeSize(dm,p,&coneSize);
703: for (i=0;i<coneSize;i++) {
704: PetscInt v;
706: DMLabelGetValue(perLabel,cone[i],&v);
707: ispe = (PetscBool)(ispe && (v==2));
708: }
709: if (ispe && coneSize) {
710: DMLabelSetValue(perLabel,p,2);
711: }
712: }
713: }
714: }
715: }
716: for (p=fStart;p<fEnd;p++) {
717: const PetscInt *support;
718: PetscInt supportSize;
719: PetscBool isbf = PETSC_FALSE;
721: DMPlexGetSupportSize(dm,p,&supportSize);
722: if (pown) {
723: PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE;
724: PetscInt i;
726: DMPlexGetSupport(dm,p,&support);
727: for (i=0;i<supportSize;i++) {
728: if (PetscLikely(PetscBTLookup(pown,support[i]-cStart))) has_owned = PETSC_TRUE;
729: else has_ghost = PETSC_TRUE;
730: }
731: isbf = (PetscBool)((supportSize == 1 && has_owned) || (supportSize > 1 && has_owned && has_ghost));
732: } else {
733: isbf = (PetscBool)(supportSize == 1);
734: }
735: if (!isbf && perLabel) {
736: const PetscInt *cone;
737: PetscInt coneSize,i;
739: DMPlexGetCone(dm,p,&cone);
740: DMPlexGetConeSize(dm,p,&coneSize);
741: isbf = PETSC_TRUE;
742: for (i=0;i<coneSize;i++) {
743: PetscInt v,d;
745: DMLabelGetValue(perLabel,cone[i],&v);
746: DMLabelGetDefaultValue(perLabel,&d);
747: isbf = (PetscBool)(isbf && v != d);
748: }
749: }
750: if (isbf) {
751: PetscBTSet(bfaces,p-fStart);
752: }
753: }
754: /* count boundary faces */
755: for (p=fStart,bf=0;p<fEnd;p++) {
756: if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) {
757: const PetscInt *support;
758: PetscInt supportSize,c;
760: DMPlexGetSupportSize(dm,p,&supportSize);
761: DMPlexGetSupport(dm,p,&support);
762: if (pown) {
763: for (c=0;c<supportSize;c++) {
764: if (PetscLikely(PetscBTLookup(pown,support[c]-cStart))) {
765: bf++;
766: }
767: }
768: } else bf += supportSize;
769: }
770: }
771: DMGetLabel(dm,bmark,&label);
772: PetscViewerASCIIPrintf(viewer,"%D\n",bf);
773: for (p=fStart;p<fEnd;p++) {
774: if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) {
775: const PetscInt *support;
776: PetscInt supportSize,c,nc = 0;
778: DMPlexGetSupportSize(dm,p,&supportSize);
779: DMPlexGetSupport(dm,p,&support);
780: if (pown) {
781: for (c=0;c<supportSize;c++) {
782: if (PetscLikely(PetscBTLookup(pown,support[c]-cStart))) {
783: fcells[nc++] = support[c];
784: }
785: }
786: } else for (c=0;c<supportSize;c++) fcells[nc++] = support[c];
787: for (c=0;c<nc;c++) {
788: const PetscInt *cone;
789: int vids[8];
790: PetscInt i,cell,cl,nv,cid = -1,mid = -1;
792: cell = fcells[c];
793: DMPlexGetCone(dm,cell,&cone);
794: for (cl=0;cl<fpc;cl++)
795: if (cone[cl] == p)
796: break;
798: /* face material id and type */
799: DMPlexGetPointMFEMCellID_Internal(dm,label,p,&mid,&cid);
800: PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);
801: /* vertex ids */
802: DMPlexGetPointMFEMVertexIDs_Internal(dm,cell,(localized && !hovec) ? coordSection : NULL,&nv,vids);
803: for (i=0;i<vpf;i++) {
804: PetscViewerASCIIPrintf(viewer," %d",vids[faces[cl*vpf+i]]);
805: }
806: PetscViewerASCIIPrintf(viewer,"\n");
807: }
808: bf = bf-nc;
809: }
810: }
811: if (bf) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Remaining boundary faces %D",bf);
812: PetscBTDestroy(&bfaces);
813: PetscFree(fcells);
814: }
816: /* mark owned vertices */
817: vown = NULL;
818: if (pown) {
819: PetscBTCreate(vEnd-vStart,&vown);
820: for (p=cStart;p<cEnd;p++) {
821: PetscInt i,closureSize,*closure = NULL;
823: if (PetscUnlikely(!PetscBTLookup(pown,p-cStart))) continue;
824: DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
825: for (i=0;i<closureSize;i++) {
826: const PetscInt pp = closure[2*i];
828: if (pp >= vStart && pp < vEnd) {
829: PetscBTSet(vown,pp-vStart);
830: }
831: }
832: DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
833: }
834: }
836: /* vertex_parents (Non-conforming meshes) */
837: parentSection = NULL;
838: if (enable_ncmesh) {
839: DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
840: }
841: if (parentSection) {
842: PetscInt vp,gvp;
844: for (vp=0,p=vStart;p<vEnd;p++) {
845: DMLabel dlabel;
846: PetscInt parent,depth;
848: if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue;
849: DMPlexGetDepthLabel(dm,&dlabel);
850: DMLabelGetValue(dlabel,p,&depth);
851: DMPlexGetTreeParent(dm,p,&parent,NULL);
852: if (parent != p) vp++;
853: }
854: MPIU_Allreduce(&vp,&gvp,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
855: if (gvp) {
856: PetscInt maxsupp;
857: PetscBool *skip = NULL;
859: PetscViewerASCIIPrintf(viewer,"\nvertex_parents\n");
860: PetscViewerASCIIPrintf(viewer,"%D\n",vp);
861: DMPlexGetMaxSizes(dm,NULL,&maxsupp);
862: PetscMalloc1(maxsupp,&skip);
863: for (p=vStart;p<vEnd;p++) {
864: DMLabel dlabel;
865: PetscInt parent;
867: if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue;
868: DMPlexGetDepthLabel(dm,&dlabel);
869: DMPlexGetTreeParent(dm,p,&parent,NULL);
870: if (parent != p) {
871: int vids[8] = { -1, -1, -1, -1, -1, -1, -1, -1 }; /* silent overzealous clang static analyzer */
872: PetscInt i,nv,ssize,n,numChildren,depth = -1;
873: const PetscInt *children;
875: DMPlexGetConeSize(dm,parent,&ssize);
876: switch (ssize) {
877: case 2: /* edge */
878: nv = 0;
879: DMPlexGetPointMFEMVertexIDs_Internal(dm,parent,localized ? coordSection : NULL,&nv,vids);
880: DMPlexInvertCell(dim,nv,vids);
881: PetscViewerASCIIPrintf(viewer,"%D",p-vStart);
882: for (i=0;i<nv;i++) {
883: PetscViewerASCIIPrintf(viewer," %D",(PetscInt)vids[i]);
884: }
885: PetscViewerASCIIPrintf(viewer,"\n");
886: vp--;
887: break;
888: case 4: /* face */
889: DMPlexGetTreeChildren(dm,parent,&numChildren,&children);
890: for (n=0;n<numChildren;n++) {
891: DMLabelGetValue(dlabel,children[n],&depth);
892: if (!depth) {
893: const PetscInt *hvsupp,*hesupp,*cone;
894: PetscInt hvsuppSize,hesuppSize,coneSize;
895: PetscInt hv = children[n],he = -1,f;
897: PetscMemzero(skip,maxsupp*sizeof(PetscBool));
898: DMPlexGetSupportSize(dm,hv,&hvsuppSize);
899: DMPlexGetSupport(dm,hv,&hvsupp);
900: for (i=0;i<hvsuppSize;i++) {
901: PetscInt ep;
902: DMPlexGetTreeParent(dm,hvsupp[i],&ep,NULL);
903: if (ep != hvsupp[i]) {
904: he = hvsupp[i];
905: } else {
906: skip[i] = PETSC_TRUE;
907: }
908: }
909: if (he == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vertex %D support size %D: hanging edge not found",hv,hvsuppSize);
910: DMPlexGetCone(dm,he,&cone);
911: vids[0] = (int)((cone[0] == hv) ? cone[1] : cone[0]);
912: DMPlexGetSupportSize(dm,he,&hesuppSize);
913: DMPlexGetSupport(dm,he,&hesupp);
914: for (f=0;f<hesuppSize;f++) {
915: PetscInt j;
917: DMPlexGetCone(dm,hesupp[f],&cone);
918: DMPlexGetConeSize(dm,hesupp[f],&coneSize);
919: for (j=0;j<coneSize;j++) {
920: PetscInt k;
921: for (k=0;k<hvsuppSize;k++) {
922: if (hvsupp[k] == cone[j]) {
923: skip[k] = PETSC_TRUE;
924: break;
925: }
926: }
927: }
928: }
929: for (i=0;i<hvsuppSize;i++) {
930: if (!skip[i]) {
931: DMPlexGetCone(dm,hvsupp[i],&cone);
932: vids[1] = (int)((cone[0] == hv) ? cone[1] : cone[0]);
933: }
934: }
935: PetscViewerASCIIPrintf(viewer,"%D",hv-vStart);
936: for (i=0;i<2;i++) {
937: PetscViewerASCIIPrintf(viewer," %D",(PetscInt)(vids[i]-vStart));
938: }
939: PetscViewerASCIIPrintf(viewer,"\n");
940: vp--;
941: }
942: }
943: break;
944: default:
945: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Don't know how to deal with support size %D",ssize);
946: }
947: }
948: }
949: PetscFree(skip);
950: }
951: if (vp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %D hanging vertices",vp);
952: }
953: PetscBTDestroy(&pown);
954: PetscBTDestroy(&vown);
956: /* vertices */
957: if (hovec) { /* higher-order meshes */
958: const char *fec;
959: PetscInt i,n;
961: PetscViewerASCIIPrintf(viewer,"\nvertices\n");
962: PetscViewerASCIIPrintf(viewer,"%D\n",vEnd-vStart);
963: PetscViewerASCIIPrintf(viewer,"nodes\n");
964: PetscObjectGetName((PetscObject)hovec,&fec);
965: PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");
966: PetscViewerASCIIPrintf(viewer,"%s\n",fec);
967: PetscViewerASCIIPrintf(viewer,"VDim: %D\n",sdim);
968: PetscViewerASCIIPrintf(viewer,"Ordering: 1\n\n"); /*Ordering::byVDIM*/
969: VecGetArrayRead(hovec,&array);
970: VecGetLocalSize(hovec,&n);
971: if (n%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Size of local coordinate vector %D incompatible with space dimension %D",n,sdim);
972: for (i=0;i<n/sdim;i++) {
973: PetscInt s;
974: for (s=0;s<sdim;s++) {
975: PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[i*sdim+s]));
976: }
977: PetscViewerASCIIPrintf(viewer,"\n");
978: }
979: VecRestoreArrayRead(hovec,&array);
980: } else {
981: VecGetLocalSize(coordinates,&nvert);
982: PetscViewerASCIIPrintf(viewer,"\nvertices\n");
983: PetscViewerASCIIPrintf(viewer,"%D\n",nvert/sdim);
984: PetscViewerASCIIPrintf(viewer,"%D\n",sdim);
985: VecGetArrayRead(coordinates,&array);
986: for (p=0;p<nvert/sdim;p++) {
987: PetscInt s;
988: for (s=0;s<sdim;s++) {
989: PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[p*sdim+s]));
990: }
991: PetscViewerASCIIPrintf(viewer,"\n");
992: }
993: VecRestoreArrayRead(coordinates,&array);
994: }
995: return(0);
996: }
998: /* dispatching, prints through the socket by prepending the mesh keyword to the usual ASCII dump */
999: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM dm, PetscViewer viewer)
1000: {
1002: PetscBool isglvis,isascii;
1007: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);
1008: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1009: if (!isglvis && !isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERGLVIS or VIEWERASCII");
1010: if (isglvis) {
1011: PetscViewer view;
1012: PetscViewerGLVisType type;
1014: PetscViewerGLVisGetType_Private(viewer,&type);
1015: PetscViewerGLVisGetDMWindow_Private(viewer,&view);
1016: if (view) { /* in the socket case, it may happen that the connection failed */
1017: if (type == PETSC_VIEWER_GLVIS_SOCKET) {
1018: PetscMPIInt size,rank;
1020: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
1021: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1022: PetscViewerASCIIPrintf(view,"parallel %d %d\nmesh\n",size,rank);
1023: }
1024: DMPlexView_GLVis_ASCII(dm,view);
1025: PetscViewerFlush(view);
1026: if (type == PETSC_VIEWER_GLVIS_SOCKET) {
1027: PetscInt dim;
1028: const char* name;
1030: PetscObjectGetName((PetscObject)dm,&name);
1031: DMGetDimension(dm,&dim);
1032: PetscViewerGLVisInitWindow_Private(view,PETSC_TRUE,dim,name);
1033: PetscBarrier((PetscObject)dm);
1034: }
1035: }
1036: } else {
1037: DMPlexView_GLVis_ASCII(dm,viewer);
1038: }
1039: return(0);
1040: }