Actual source code: plexglvis.c
1: #include <petsc/private/glvisviewerimpl.h>
2: #include <petsc/private/petscimpl.h>
3: #include <petsc/private/dmpleximpl.h>
4: #include <petscbt.h>
5: #include <petscdmplex.h>
6: #include <petscsf.h>
7: #include <petscds.h>
9: typedef struct {
10: PetscInt nf;
11: VecScatter *scctx;
12: } GLVisViewerCtx;
14: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
15: {
16: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
17: PetscInt i;
19: for (i = 0; i < ctx->nf; i++) VecScatterDestroy(&ctx->scctx[i]);
20: PetscFree(ctx->scctx);
21: PetscFree(vctx);
22: return 0;
23: }
25: static PetscErrorCode DMPlexSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
26: {
27: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
28: PetscInt f;
30: for (f = 0; f < nf; f++) {
31: VecScatterBegin(ctx->scctx[f], (Vec)oX, (Vec)oXfield[f], INSERT_VALUES, SCATTER_FORWARD);
32: VecScatterEnd(ctx->scctx[f], (Vec)oX, (Vec)oXfield[f], INSERT_VALUES, SCATTER_FORWARD);
33: }
34: return 0;
35: }
37: /* for FEM, it works for H1 fields only and extracts dofs at cell vertices, discarding any other dof */
38: PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject odm, PetscViewer viewer)
39: {
40: DM dm = (DM)odm;
41: Vec xlocal, xfield, *Ufield;
42: PetscDS ds;
43: IS globalNum, isfield;
44: PetscBT vown;
45: char **fieldname = NULL, **fec_type = NULL;
46: const PetscInt *gNum;
47: PetscInt *nlocal, *bs, *idxs, *dims;
48: PetscInt f, maxfields, nfields, c, totc, totdofs, Nv, cum, i;
49: PetscInt dim, cStart, cEnd, vStart, vEnd;
50: GLVisViewerCtx *ctx;
51: PetscSection s;
53: DMGetDimension(dm, &dim);
54: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
55: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
56: PetscObjectQuery((PetscObject)dm, "_glvis_plex_gnum", (PetscObject *)&globalNum);
57: if (!globalNum) {
58: DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, &globalNum);
59: PetscObjectCompose((PetscObject)dm, "_glvis_plex_gnum", (PetscObject)globalNum);
60: PetscObjectDereference((PetscObject)globalNum);
61: }
62: ISGetIndices(globalNum, &gNum);
63: PetscBTCreate(vEnd - vStart, &vown);
64: for (c = cStart, totc = 0; c < cEnd; c++) {
65: if (gNum[c - cStart] >= 0) {
66: PetscInt i, numPoints, *points = NULL;
68: totc++;
69: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &numPoints, &points);
70: for (i = 0; i < numPoints * 2; i += 2) {
71: if ((points[i] >= vStart) && (points[i] < vEnd)) PetscBTSet(vown, points[i] - vStart);
72: }
73: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &numPoints, &points);
74: }
75: }
76: for (f = 0, Nv = 0; f < vEnd - vStart; f++)
77: if (PetscLikely(PetscBTLookup(vown, f))) Nv++;
79: DMCreateLocalVector(dm, &xlocal);
80: VecGetLocalSize(xlocal, &totdofs);
81: DMGetLocalSection(dm, &s);
82: PetscSectionGetNumFields(s, &nfields);
83: for (f = 0, maxfields = 0; f < nfields; f++) {
84: PetscInt bs;
86: PetscSectionGetFieldComponents(s, f, &bs);
87: maxfields += bs;
88: }
89: PetscCalloc7(maxfields, &fieldname, maxfields, &nlocal, maxfields, &bs, maxfields, &dims, maxfields, &fec_type, totdofs, &idxs, maxfields, &Ufield);
90: PetscNew(&ctx);
91: PetscCalloc1(maxfields, &ctx->scctx);
92: DMGetDS(dm, &ds);
93: if (ds) {
94: for (f = 0; f < nfields; f++) {
95: const char *fname;
96: char name[256];
97: PetscObject disc;
98: size_t len;
100: PetscSectionGetFieldName(s, f, &fname);
101: PetscStrlen(fname, &len);
102: if (len) {
103: PetscStrcpy(name, fname);
104: } else {
105: PetscSNPrintf(name, 256, "Field%" PetscInt_FMT, f);
106: }
107: PetscDSGetDiscretization(ds, f, &disc);
108: if (disc) {
109: PetscClassId id;
110: PetscInt Nc;
111: char fec[64];
113: PetscObjectGetClassId(disc, &id);
114: if (id == PETSCFE_CLASSID) {
115: PetscFE fem = (PetscFE)disc;
116: PetscDualSpace sp;
117: PetscDualSpaceType spname;
118: PetscInt order;
119: PetscBool islag, continuous, H1 = PETSC_TRUE;
121: PetscFEGetNumComponents(fem, &Nc);
122: PetscFEGetDualSpace(fem, &sp);
123: PetscDualSpaceGetType(sp, &spname);
124: PetscStrcmp(spname, PETSCDUALSPACELAGRANGE, &islag);
126: PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
127: PetscDualSpaceGetOrder(sp, &order);
128: if (continuous && order > 0) { /* no support for high-order viz, still have to figure out the numbering */
129: PetscSNPrintf(fec, 64, "FiniteElementCollection: H1_%" PetscInt_FMT "D_P1", dim);
130: } else {
132: H1 = PETSC_FALSE;
133: PetscSNPrintf(fec, 64, "FiniteElementCollection: L2_%" PetscInt_FMT "D_P%" PetscInt_FMT, dim, order);
134: }
135: PetscStrallocpy(name, &fieldname[ctx->nf]);
136: bs[ctx->nf] = Nc;
137: dims[ctx->nf] = dim;
138: if (H1) {
139: nlocal[ctx->nf] = Nc * Nv;
140: PetscStrallocpy(fec, &fec_type[ctx->nf]);
141: VecCreateSeq(PETSC_COMM_SELF, Nv * Nc, &xfield);
142: for (i = 0, cum = 0; i < vEnd - vStart; i++) {
143: PetscInt j, off;
145: if (PetscUnlikely(!PetscBTLookup(vown, i))) continue;
146: PetscSectionGetFieldOffset(s, i + vStart, f, &off);
147: for (j = 0; j < Nc; j++) idxs[cum++] = off + j;
148: }
149: ISCreateGeneral(PetscObjectComm((PetscObject)xlocal), Nv * Nc, idxs, PETSC_USE_POINTER, &isfield);
150: } else {
151: nlocal[ctx->nf] = Nc * totc;
152: PetscStrallocpy(fec, &fec_type[ctx->nf]);
153: VecCreateSeq(PETSC_COMM_SELF, Nc * totc, &xfield);
154: for (i = 0, cum = 0; i < cEnd - cStart; i++) {
155: PetscInt j, off;
157: if (PetscUnlikely(gNum[i] < 0)) continue;
158: PetscSectionGetFieldOffset(s, i + cStart, f, &off);
159: for (j = 0; j < Nc; j++) idxs[cum++] = off + j;
160: }
161: ISCreateGeneral(PetscObjectComm((PetscObject)xlocal), totc * Nc, idxs, PETSC_USE_POINTER, &isfield);
162: }
163: VecScatterCreate(xlocal, isfield, xfield, NULL, &ctx->scctx[ctx->nf]);
164: VecDestroy(&xfield);
165: ISDestroy(&isfield);
166: ctx->nf++;
167: } else if (id == PETSCFV_CLASSID) {
168: PetscInt c;
170: PetscFVGetNumComponents((PetscFV)disc, &Nc);
171: PetscSNPrintf(fec, 64, "FiniteElementCollection: L2_%" PetscInt_FMT "D_P0", dim);
172: for (c = 0; c < Nc; c++) {
173: char comp[256];
174: PetscSNPrintf(comp, 256, "%s-Comp%" PetscInt_FMT, name, c);
175: PetscStrallocpy(comp, &fieldname[ctx->nf]);
176: bs[ctx->nf] = 1; /* Does PetscFV support components with different block size? */
177: nlocal[ctx->nf] = totc;
178: dims[ctx->nf] = dim;
179: PetscStrallocpy(fec, &fec_type[ctx->nf]);
180: VecCreateSeq(PETSC_COMM_SELF, totc, &xfield);
181: for (i = 0, cum = 0; i < cEnd - cStart; i++) {
182: PetscInt off;
184: if (PetscUnlikely(gNum[i]) < 0) continue;
185: PetscSectionGetFieldOffset(s, i + cStart, f, &off);
186: idxs[cum++] = off + c;
187: }
188: ISCreateGeneral(PetscObjectComm((PetscObject)xlocal), totc, idxs, PETSC_USE_POINTER, &isfield);
189: VecScatterCreate(xlocal, isfield, xfield, NULL, &ctx->scctx[ctx->nf]);
190: VecDestroy(&xfield);
191: ISDestroy(&isfield);
192: ctx->nf++;
193: }
194: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
195: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing discretization for field %" PetscInt_FMT, f);
196: }
197: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Needs a DS attached to the DM");
198: PetscBTDestroy(&vown);
199: VecDestroy(&xlocal);
200: ISRestoreIndices(globalNum, &gNum);
202: /* create work vectors */
203: for (f = 0; f < ctx->nf; f++) {
204: VecCreateMPI(PetscObjectComm((PetscObject)dm), nlocal[f], PETSC_DECIDE, &Ufield[f]);
205: PetscObjectSetName((PetscObject)Ufield[f], fieldname[f]);
206: VecSetBlockSize(Ufield[f], bs[f]);
207: VecSetDM(Ufield[f], dm);
208: }
210: /* customize the viewer */
211: PetscViewerGLVisSetFields(viewer, ctx->nf, (const char **)fec_type, dims, DMPlexSampleGLVisFields_Private, (PetscObject *)Ufield, ctx, DestroyGLVisViewerCtx_Private);
212: for (f = 0; f < ctx->nf; f++) {
213: PetscFree(fieldname[f]);
214: PetscFree(fec_type[f]);
215: VecDestroy(&Ufield[f]);
216: }
217: PetscFree7(fieldname, nlocal, bs, dims, fec_type, idxs, Ufield);
218: return 0;
219: }
221: typedef enum {
222: MFEM_POINT = 0,
223: MFEM_SEGMENT,
224: MFEM_TRIANGLE,
225: MFEM_SQUARE,
226: MFEM_TETRAHEDRON,
227: MFEM_CUBE,
228: MFEM_PRISM,
229: MFEM_UNDEF
230: } MFEM_cid;
232: MFEM_cid mfem_table_cid[4][7] = {
233: {MFEM_POINT, MFEM_UNDEF, MFEM_UNDEF, MFEM_UNDEF, MFEM_UNDEF, MFEM_UNDEF, MFEM_UNDEF},
234: {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_UNDEF, MFEM_UNDEF, MFEM_UNDEF, MFEM_UNDEF},
235: {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_TRIANGLE, MFEM_SQUARE, MFEM_UNDEF, MFEM_UNDEF},
236: {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_UNDEF, MFEM_TETRAHEDRON, MFEM_PRISM, MFEM_CUBE }
237: };
239: MFEM_cid mfem_table_cid_unint[4][9] = {
240: {MFEM_POINT, MFEM_UNDEF, MFEM_UNDEF, MFEM_UNDEF, MFEM_UNDEF, MFEM_UNDEF, MFEM_PRISM, MFEM_UNDEF, MFEM_UNDEF},
241: {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_UNDEF, MFEM_UNDEF, MFEM_UNDEF, MFEM_PRISM, MFEM_UNDEF, MFEM_UNDEF},
242: {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_TRIANGLE, MFEM_SQUARE, MFEM_UNDEF, MFEM_PRISM, MFEM_UNDEF, MFEM_UNDEF},
243: {MFEM_POINT, MFEM_UNDEF, MFEM_SEGMENT, MFEM_UNDEF, MFEM_TETRAHEDRON, MFEM_UNDEF, MFEM_PRISM, MFEM_UNDEF, MFEM_CUBE }
244: };
246: static PetscErrorCode DMPlexGetPointMFEMCellID_Internal(DM dm, DMLabel label, PetscInt minl, PetscInt p, PetscInt *mid, PetscInt *cid)
247: {
248: DMLabel dlabel;
249: PetscInt depth, csize, pdepth, dim;
251: DMPlexGetDepthLabel(dm, &dlabel);
252: DMLabelGetValue(dlabel, p, &pdepth);
253: DMPlexGetConeSize(dm, p, &csize);
254: DMPlexGetDepth(dm, &depth);
255: DMGetDimension(dm, &dim);
256: if (label) {
257: DMLabelGetValue(label, p, mid);
258: *mid = *mid - minl + 1; /* MFEM does not like negative markers */
259: } else *mid = 1;
260: if (depth >= 0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
264: *cid = mfem_table_cid_unint[dim][csize];
265: } else {
268: *cid = mfem_table_cid[pdepth][csize];
269: }
270: return 0;
271: }
273: static PetscErrorCode DMPlexGetPointMFEMVertexIDs_Internal(DM dm, PetscInt p, PetscSection csec, PetscInt *nv, PetscInt vids[])
274: {
275: 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) PetscSectionGetDof(csec, p, &dof);
288: }
289: if (!dof) {
290: DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &numPoints, &points);
291: for (i = 0, q = 0; i < numPoints * 2; i += 2)
292: if ((points[i] >= vStart) && (points[i] < vEnd)) vids[q++] = points[i] - vStart + off;
293: DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &numPoints, &points);
294: } else {
295: PetscSectionGetOffset(csec, p, &off);
296: PetscSectionGetDof(csec, p, &dof);
297: for (q = 0; q < dof / sdim; q++) vids[q] = off / sdim + q;
298: }
299: *nv = q;
300: return 0;
301: }
303: static PetscErrorCode GLVisCreateFE(PetscFE femIn, char name[32], PetscFE *fem, IS *perm)
304: {
305: DM K;
306: PetscSpace P;
307: PetscDualSpace Q;
308: PetscQuadrature q, fq;
309: PetscInt dim, deg, dof;
310: DMPolytopeType ptype;
311: PetscBool isSimplex, isTensor;
312: PetscBool continuity = PETSC_FALSE;
313: PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI;
314: PetscBool endpoint = PETSC_TRUE;
315: MPI_Comm comm;
317: PetscObjectGetComm((PetscObject)femIn, &comm);
318: PetscFEGetBasisSpace(femIn, &P);
319: PetscFEGetDualSpace(femIn, &Q);
320: PetscDualSpaceGetDM(Q, &K);
321: DMGetDimension(K, &dim);
322: PetscSpaceGetDegree(P, °, NULL);
323: PetscSpaceGetNumComponents(P, &dof);
324: DMPlexGetCellType(K, 0, &ptype);
325: switch (ptype) {
326: case DM_POLYTOPE_QUADRILATERAL:
327: case DM_POLYTOPE_HEXAHEDRON:
328: isSimplex = PETSC_FALSE;
329: break;
330: default:
331: isSimplex = PETSC_TRUE;
332: break;
333: }
334: isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
335: if (isSimplex) deg = PetscMin(deg, 3); /* Permutation not coded for degree higher than 3 */
336: /* Create space */
337: PetscSpaceCreate(comm, &P);
338: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
339: PetscSpacePolynomialSetTensor(P, isTensor);
340: PetscSpaceSetNumComponents(P, dof);
341: PetscSpaceSetNumVariables(P, dim);
342: PetscSpaceSetDegree(P, deg, PETSC_DETERMINE);
343: PetscSpaceSetUp(P);
344: /* Create dual space */
345: PetscDualSpaceCreate(comm, &Q);
346: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
347: PetscDualSpaceLagrangeSetTensor(Q, isTensor);
348: PetscDualSpaceLagrangeSetContinuity(Q, continuity);
349: PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0);
350: PetscDualSpaceSetNumComponents(Q, dof);
351: PetscDualSpaceSetOrder(Q, deg);
352: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K);
353: PetscDualSpaceSetDM(Q, K);
354: DMDestroy(&K);
355: PetscDualSpaceSetUp(Q);
356: /* Create quadrature */
357: if (isSimplex) {
358: PetscDTStroudConicalQuadrature(dim, 1, deg + 1, -1, +1, &q);
359: PetscDTStroudConicalQuadrature(dim - 1, 1, deg + 1, -1, +1, &fq);
360: } else {
361: PetscDTGaussTensorQuadrature(dim, 1, deg + 1, -1, +1, &q);
362: PetscDTGaussTensorQuadrature(dim - 1, 1, deg + 1, -1, +1, &fq);
363: }
364: /* Create finite element */
365: PetscFECreate(comm, fem);
366: PetscSNPrintf(name, 32, "L2_T1_%" PetscInt_FMT "D_P%" PetscInt_FMT, dim, deg);
367: PetscObjectSetName((PetscObject)*fem, name);
368: PetscFESetType(*fem, PETSCFEBASIC);
369: PetscFESetNumComponents(*fem, dof);
370: PetscFESetBasisSpace(*fem, P);
371: PetscFESetDualSpace(*fem, Q);
372: PetscFESetQuadrature(*fem, q);
373: PetscFESetFaceQuadrature(*fem, fq);
374: PetscFESetUp(*fem);
376: /* Both MFEM and PETSc are lexicographic, but PLEX stores the swapped cone */
377: *perm = NULL;
378: if (isSimplex && dim == 3) {
379: PetscInt celldofs, *pidx;
381: PetscDualSpaceGetDimension(Q, &celldofs);
382: celldofs /= dof;
383: PetscMalloc1(celldofs, &pidx);
384: switch (celldofs) {
385: case 4:
386: pidx[0] = 2;
387: pidx[1] = 0;
388: pidx[2] = 1;
389: pidx[3] = 3;
390: break;
391: case 10:
392: pidx[0] = 5;
393: pidx[1] = 3;
394: pidx[2] = 0;
395: pidx[3] = 4;
396: pidx[4] = 1;
397: pidx[5] = 2;
398: pidx[6] = 8;
399: pidx[7] = 6;
400: pidx[8] = 7;
401: pidx[9] = 9;
402: break;
403: case 20:
404: pidx[0] = 9;
405: pidx[1] = 7;
406: pidx[2] = 4;
407: pidx[3] = 0;
408: pidx[4] = 8;
409: pidx[5] = 5;
410: pidx[6] = 1;
411: pidx[7] = 6;
412: pidx[8] = 2;
413: pidx[9] = 3;
414: pidx[10] = 15;
415: pidx[11] = 13;
416: pidx[12] = 10;
417: pidx[13] = 14;
418: pidx[14] = 11;
419: pidx[15] = 12;
420: pidx[16] = 18;
421: pidx[17] = 16;
422: pidx[18] = 17;
423: pidx[19] = 19;
424: break;
425: default:
426: SETERRQ(comm, PETSC_ERR_SUP, "Unhandled degree,dof pair %" PetscInt_FMT ",%" PetscInt_FMT, deg, celldofs);
427: break;
428: }
429: ISCreateBlock(PETSC_COMM_SELF, dof, celldofs, pidx, PETSC_OWN_POINTER, perm);
430: }
432: /* Cleanup */
433: PetscSpaceDestroy(&P);
434: PetscDualSpaceDestroy(&Q);
435: PetscQuadratureDestroy(&q);
436: PetscQuadratureDestroy(&fq);
437: return 0;
438: }
440: /*
441: ASCII visualization/dump: full support for simplices and tensor product cells. It supports AMR
442: Higher order meshes are also supported
443: */
444: static PetscErrorCode DMPlexView_GLVis_ASCII(DM dm, PetscViewer viewer)
445: {
446: DMLabel label;
447: PetscSection coordSection, coordSectionCell, parentSection, hoSection = NULL;
448: Vec coordinates, coordinatesCell, hovec;
449: const PetscScalar *array;
450: PetscInt bf, p, sdim, dim, depth, novl, minl;
451: PetscInt cStart, cEnd, vStart, vEnd, nvert;
452: PetscMPIInt size;
453: PetscBool localized, isascii;
454: PetscBool enable_mfem, enable_boundary, enable_ncmesh, view_ovl = PETSC_FALSE;
455: PetscBT pown, vown;
456: PetscContainer glvis_container;
457: PetscBool cellvertex = PETSC_FALSE, enabled = PETSC_TRUE;
458: PetscBool enable_emark, enable_bmark;
459: const char *fmt;
460: char emark[64] = "", bmark[64] = "";
464: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
466: MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size);
468: DMGetDimension(dm, &dim);
469: DMPlexGetDepth(dm, &depth);
471: /* get container: determines if a process visualizes is portion of the data or not */
472: PetscObjectQuery((PetscObject)viewer, "_glvis_info_container", (PetscObject *)&glvis_container);
474: {
475: PetscViewerGLVisInfo glvis_info;
476: PetscContainerGetPointer(glvis_container, (void **)&glvis_info);
477: enabled = glvis_info->enabled;
478: fmt = glvis_info->fmt;
479: }
481: /* Users can attach a coordinate vector to the DM in case they have a higher-order mesh */
482: PetscObjectQuery((PetscObject)dm, "_glvis_mesh_coords", (PetscObject *)&hovec);
483: PetscObjectReference((PetscObject)hovec);
484: if (!hovec) {
485: DM cdm;
486: PetscFE disc;
487: PetscClassId classid;
489: DMGetCoordinateDM(dm, &cdm);
490: DMGetField(cdm, 0, NULL, (PetscObject *)&disc);
491: PetscObjectGetClassId((PetscObject)disc, &classid);
492: if (classid == PETSCFE_CLASSID) {
493: DM hocdm;
494: PetscFE hodisc;
495: Vec vec;
496: Mat mat;
497: char name[32], fec_type[64];
498: IS perm = NULL;
500: GLVisCreateFE(disc, name, &hodisc, &perm);
501: DMClone(cdm, &hocdm);
502: DMSetField(hocdm, 0, NULL, (PetscObject)hodisc);
503: PetscFEDestroy(&hodisc);
504: DMCreateDS(hocdm);
506: DMGetCoordinates(dm, &vec);
507: DMCreateGlobalVector(hocdm, &hovec);
508: DMCreateInterpolation(cdm, hocdm, &mat, NULL);
509: MatInterpolate(mat, vec, hovec);
510: MatDestroy(&mat);
511: DMGetLocalSection(hocdm, &hoSection);
512: PetscSectionSetClosurePermutation(hoSection, (PetscObject)hocdm, depth, perm);
513: ISDestroy(&perm);
514: DMDestroy(&hocdm);
515: PetscSNPrintf(fec_type, sizeof(fec_type), "FiniteElementCollection: %s", name);
516: PetscObjectSetName((PetscObject)hovec, fec_type);
517: }
518: }
520: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
521: DMPlexGetGhostCellStratum(dm, &p, NULL);
522: if (p >= 0) cEnd = p;
523: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
524: DMGetCoordinatesLocalized(dm, &localized);
525: DMGetCoordinateSection(dm, &coordSection);
526: DMGetCoordinateDim(dm, &sdim);
527: DMGetCoordinatesLocal(dm, &coordinates);
530: /*
531: a couple of sections of the mesh specification are disabled
532: - 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)
533: - vertex_parents: used for non-conforming meshes only when we want to use MFEM as a discretization package
534: and be able to derefine the mesh (MFEM does not currently have to ability to read ncmeshes in parallel)
535: */
536: enable_boundary = PETSC_FALSE;
537: enable_ncmesh = PETSC_FALSE;
538: enable_mfem = PETSC_FALSE;
539: enable_emark = PETSC_FALSE;
540: enable_bmark = PETSC_FALSE;
541: /* I'm tired of problems with negative values in the markers, disable them */
542: PetscOptionsBegin(PetscObjectComm((PetscObject)dm), ((PetscObject)dm)->prefix, "GLVis PetscViewer DMPlex Options", "PetscViewer");
543: PetscOptionsBool("-viewer_glvis_dm_plex_enable_boundary", "Enable boundary section in mesh representation", NULL, enable_boundary, &enable_boundary, NULL);
544: PetscOptionsBool("-viewer_glvis_dm_plex_enable_ncmesh", "Enable vertex_parents section in mesh representation (allows derefinement)", NULL, enable_ncmesh, &enable_ncmesh, NULL);
545: PetscOptionsBool("-viewer_glvis_dm_plex_enable_mfem", "Dump a mesh that can be used with MFEM's FiniteElementSpaces", NULL, enable_mfem, &enable_mfem, NULL);
546: PetscOptionsBool("-viewer_glvis_dm_plex_overlap", "Include overlap region in local meshes", NULL, view_ovl, &view_ovl, NULL);
547: PetscOptionsString("-viewer_glvis_dm_plex_emarker", "String for the material id label", NULL, emark, emark, sizeof(emark), &enable_emark);
548: PetscOptionsString("-viewer_glvis_dm_plex_bmarker", "String for the boundary id label", NULL, bmark, bmark, sizeof(bmark), &enable_bmark);
549: PetscOptionsEnd();
550: if (enable_bmark) enable_boundary = PETSC_TRUE;
552: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
555: "Mesh must be interpolated. "
556: "Alternatively, run with -viewer_glvis_dm_plex_enable_boundary 0");
558: "Mesh must be interpolated. "
559: "Alternatively, run with -viewer_glvis_dm_plex_enable_ncmesh 0");
560: if (depth >= 0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */
562: cellvertex = PETSC_TRUE;
563: }
565: /* Identify possible cells in the overlap */
566: novl = 0;
567: pown = NULL;
568: if (size > 1) {
569: IS globalNum = NULL;
570: const PetscInt *gNum;
571: PetscBool ovl = PETSC_FALSE;
573: PetscObjectQuery((PetscObject)dm, "_glvis_plex_gnum", (PetscObject *)&globalNum);
574: if (!globalNum) {
575: if (view_ovl) {
576: ISCreateStride(PetscObjectComm((PetscObject)dm), cEnd - cStart, 0, 1, &globalNum);
577: } else {
578: DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, &globalNum);
579: }
580: PetscObjectCompose((PetscObject)dm, "_glvis_plex_gnum", (PetscObject)globalNum);
581: PetscObjectDereference((PetscObject)globalNum);
582: }
583: ISGetIndices(globalNum, &gNum);
584: for (p = cStart; p < cEnd; p++) {
585: if (gNum[p - cStart] < 0) {
586: ovl = PETSC_TRUE;
587: novl++;
588: }
589: }
590: if (ovl) {
591: /* it may happen that pown get not destroyed, if the user closes the window while this function is running.
592: TODO: garbage collector? attach pown to dm? */
593: PetscBTCreate(cEnd - cStart, &pown);
594: for (p = cStart; p < cEnd; p++) {
595: if (gNum[p - cStart] < 0) continue;
596: else PetscBTSet(pown, p - cStart);
597: }
598: }
599: ISRestoreIndices(globalNum, &gNum);
600: }
602: /* vertex_parents (Non-conforming meshes) */
603: parentSection = NULL;
604: if (enable_ncmesh) {
605: DMPlexGetTree(dm, &parentSection, NULL, NULL, NULL, NULL);
606: enable_ncmesh = (PetscBool)(enable_ncmesh && parentSection);
607: }
608: /* return if this process is disabled */
609: if (!enabled) {
610: PetscViewerASCIIPrintf(viewer, "MFEM mesh %s\n", enable_ncmesh ? "v1.1" : "v1.0");
611: PetscViewerASCIIPrintf(viewer, "\ndimension\n");
612: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", dim);
613: PetscViewerASCIIPrintf(viewer, "\nelements\n");
614: PetscViewerASCIIPrintf(viewer, "0\n");
615: PetscViewerASCIIPrintf(viewer, "\nboundary\n");
616: PetscViewerASCIIPrintf(viewer, "0\n");
617: PetscViewerASCIIPrintf(viewer, "\nvertices\n");
618: PetscViewerASCIIPrintf(viewer, "0\n");
619: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", sdim);
620: PetscBTDestroy(&pown);
621: VecDestroy(&hovec);
622: return 0;
623: }
625: if (enable_mfem) {
626: if (localized && !hovec) { /* we need to generate a vector of L2 coordinates, as this is how MFEM handles periodic meshes */
627: PetscInt vpc = 0;
628: char fec[64];
629: PetscInt vids[8] = {0, 1, 2, 3, 4, 5, 6, 7};
630: PetscInt hexv[8] = {0, 1, 3, 2, 4, 5, 7, 6}, tetv[4] = {0, 1, 2, 3};
631: PetscInt quadv[8] = {0, 1, 3, 2}, triv[3] = {0, 1, 2};
632: PetscInt *dof = NULL;
633: PetscScalar *array, *ptr;
635: PetscSNPrintf(fec, sizeof(fec), "FiniteElementCollection: L2_T1_%" PetscInt_FMT "D_P1", dim);
636: if (cEnd - cStart) {
637: PetscInt fpc;
639: DMPlexGetConeSize(dm, cStart, &fpc);
640: switch (dim) {
641: case 1:
642: vpc = 2;
643: dof = hexv;
644: break;
645: case 2:
646: switch (fpc) {
647: case 3:
648: vpc = 3;
649: dof = triv;
650: break;
651: case 4:
652: vpc = 4;
653: dof = quadv;
654: break;
655: default:
656: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unhandled case: faces per cell %" PetscInt_FMT, fpc);
657: }
658: break;
659: case 3:
660: switch (fpc) {
661: case 4: /* TODO: still need to understand L2 ordering for tets */
662: vpc = 4;
663: dof = tetv;
664: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unhandled tethraedral case");
665: case 6:
667: vpc = 8;
668: dof = hexv;
669: break;
670: case 8:
672: vpc = 8;
673: dof = hexv;
674: break;
675: default:
676: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unhandled case: faces per cell %" PetscInt_FMT, fpc);
677: }
678: break;
679: default:
680: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unhandled dim");
681: }
682: DMPlexReorderCell(dm, cStart, vids);
683: }
685: VecCreateSeq(PETSC_COMM_SELF, (cEnd - cStart - novl) * vpc * sdim, &hovec);
686: PetscObjectSetName((PetscObject)hovec, fec);
687: VecGetArray(hovec, &array);
688: ptr = array;
689: for (p = cStart; p < cEnd; p++) {
690: PetscInt csize, v, d;
691: PetscScalar *vals = NULL;
693: if (PetscUnlikely(pown && !PetscBTLookup(pown, p - cStart))) continue;
694: DMPlexVecGetClosure(dm, coordSection, coordinates, p, &csize, &vals);
696: for (v = 0; v < vpc; v++) {
697: for (d = 0; d < sdim; d++) ptr[sdim * dof[v] + d] = vals[sdim * vids[v] + d];
698: }
699: ptr += vpc * sdim;
700: DMPlexVecRestoreClosure(dm, coordSection, coordinates, p, &csize, &vals);
701: }
702: VecRestoreArray(hovec, &array);
703: }
704: }
705: /* if we have high-order coordinates in 3D, we need to specify the boundary */
706: if (hovec && dim == 3) enable_boundary = PETSC_TRUE;
708: /* header */
709: PetscViewerASCIIPrintf(viewer, "MFEM mesh %s\n", enable_ncmesh ? "v1.1" : "v1.0");
711: /* topological dimension */
712: PetscViewerASCIIPrintf(viewer, "\ndimension\n");
713: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", dim);
715: /* elements */
716: minl = 1;
717: label = NULL;
718: if (enable_emark) {
719: PetscInt lminl = PETSC_MAX_INT;
721: DMGetLabel(dm, emark, &label);
722: if (label) {
723: IS vals;
724: PetscInt ldef;
726: DMLabelGetDefaultValue(label, &ldef);
727: DMLabelGetValueIS(label, &vals);
728: ISGetMinMax(vals, &lminl, NULL);
729: ISDestroy(&vals);
730: lminl = PetscMin(ldef, lminl);
731: }
732: MPIU_Allreduce(&lminl, &minl, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
733: if (minl == PETSC_MAX_INT) minl = 1;
734: }
735: PetscViewerASCIIPrintf(viewer, "\nelements\n");
736: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", cEnd - cStart - novl);
737: for (p = cStart; p < cEnd; p++) {
738: PetscInt vids[8];
739: PetscInt i, nv = 0, cid = -1, mid = 1;
741: if (PetscUnlikely(pown && !PetscBTLookup(pown, p - cStart))) continue;
742: DMPlexGetPointMFEMCellID_Internal(dm, label, minl, p, &mid, &cid);
743: DMPlexGetPointMFEMVertexIDs_Internal(dm, p, (localized && !hovec) ? coordSection : NULL, &nv, vids);
744: DMPlexReorderCell(dm, p, vids);
745: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT, mid, cid);
746: for (i = 0; i < nv; i++) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, vids[i]);
747: PetscViewerASCIIPrintf(viewer, "\n");
748: }
750: /* boundary */
751: PetscViewerASCIIPrintf(viewer, "\nboundary\n");
752: if (!enable_boundary) {
753: PetscViewerASCIIPrintf(viewer, "0\n");
754: } else {
755: DMLabel perLabel;
756: PetscBT bfaces;
757: PetscInt fStart, fEnd, *fcells;
759: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
760: PetscBTCreate(fEnd - fStart, &bfaces);
761: DMPlexGetMaxSizes(dm, NULL, &p);
762: PetscMalloc1(p, &fcells);
763: DMGetLabel(dm, "glvis_periodic_cut", &perLabel);
764: if (!perLabel && localized) { /* this periodic cut can be moved up to DMPlex setup */
765: DMCreateLabel(dm, "glvis_periodic_cut");
766: DMGetLabel(dm, "glvis_periodic_cut", &perLabel);
767: DMLabelSetDefaultValue(perLabel, 1);
768: DMGetCellCoordinateSection(dm, &coordSectionCell);
769: DMGetCellCoordinatesLocal(dm, &coordinatesCell);
770: for (p = cStart; p < cEnd; p++) {
771: DMPolytopeType cellType;
772: PetscInt dof;
774: DMPlexGetCellType(dm, p, &cellType);
775: PetscSectionGetDof(coordSectionCell, p, &dof);
776: if (dof) {
777: PetscInt uvpc, v, csize, csizeCell, cellClosureSize, *cellClosure = NULL, *vidxs = NULL;
778: PetscScalar *vals = NULL, *valsCell = NULL;
780: uvpc = DMPolytopeTypeGetNumVertices(cellType);
782: DMPlexVecGetClosure(dm, coordSection, coordinates, p, &csize, &vals);
783: DMPlexVecGetClosure(dm, coordSectionCell, coordinatesCell, p, &csizeCell, &valsCell);
785: DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &cellClosureSize, &cellClosure);
786: for (v = 0; v < cellClosureSize; v++)
787: if (cellClosure[2 * v] >= vStart && cellClosure[2 * v] < vEnd) {
788: vidxs = cellClosure + 2 * v;
789: break;
790: }
792: for (v = 0; v < uvpc; v++) {
793: PetscInt s;
795: for (s = 0; s < sdim; s++) {
796: if (PetscAbsScalar(vals[v * sdim + s] - valsCell[v * sdim + s]) > PETSC_MACHINE_EPSILON) DMLabelSetValue(perLabel, vidxs[2 * v], 2);
797: }
798: }
799: DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &cellClosureSize, &cellClosure);
800: DMPlexVecRestoreClosure(dm, coordSection, coordinates, p, &csize, &vals);
801: DMPlexVecRestoreClosure(dm, coordSectionCell, coordinatesCell, p, &csizeCell, &valsCell);
802: }
803: }
804: if (dim > 1) {
805: PetscInt eEnd, eStart;
807: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
808: for (p = eStart; p < eEnd; p++) {
809: const PetscInt *cone;
810: PetscInt coneSize, i;
811: PetscBool ispe = PETSC_TRUE;
813: DMPlexGetCone(dm, p, &cone);
814: DMPlexGetConeSize(dm, p, &coneSize);
815: for (i = 0; i < coneSize; i++) {
816: PetscInt v;
818: DMLabelGetValue(perLabel, cone[i], &v);
819: ispe = (PetscBool)(ispe && (v == 2));
820: }
821: if (ispe && coneSize) {
822: PetscInt ch, numChildren;
823: const PetscInt *children;
825: DMLabelSetValue(perLabel, p, 2);
826: DMPlexGetTreeChildren(dm, p, &numChildren, &children);
827: for (ch = 0; ch < numChildren; ch++) DMLabelSetValue(perLabel, children[ch], 2);
828: }
829: }
830: if (dim > 2) {
831: for (p = fStart; p < fEnd; p++) {
832: const PetscInt *cone;
833: PetscInt coneSize, i;
834: PetscBool ispe = PETSC_TRUE;
836: DMPlexGetCone(dm, p, &cone);
837: DMPlexGetConeSize(dm, p, &coneSize);
838: for (i = 0; i < coneSize; i++) {
839: PetscInt v;
841: DMLabelGetValue(perLabel, cone[i], &v);
842: ispe = (PetscBool)(ispe && (v == 2));
843: }
844: if (ispe && coneSize) {
845: PetscInt ch, numChildren;
846: const PetscInt *children;
848: DMLabelSetValue(perLabel, p, 2);
849: DMPlexGetTreeChildren(dm, p, &numChildren, &children);
850: for (ch = 0; ch < numChildren; ch++) DMLabelSetValue(perLabel, children[ch], 2);
851: }
852: }
853: }
854: }
855: }
856: for (p = fStart; p < fEnd; p++) {
857: const PetscInt *support;
858: PetscInt supportSize;
859: PetscBool isbf = PETSC_FALSE;
861: DMPlexGetSupportSize(dm, p, &supportSize);
862: if (pown) {
863: PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE;
864: PetscInt i;
866: DMPlexGetSupport(dm, p, &support);
867: for (i = 0; i < supportSize; i++) {
868: if (PetscLikely(PetscBTLookup(pown, support[i] - cStart))) has_owned = PETSC_TRUE;
869: else has_ghost = PETSC_TRUE;
870: }
871: isbf = (PetscBool)((supportSize == 1 && has_owned) || (supportSize > 1 && has_owned && has_ghost));
872: } else {
873: isbf = (PetscBool)(supportSize == 1);
874: }
875: if (!isbf && perLabel) {
876: const PetscInt *cone;
877: PetscInt coneSize, i;
879: DMPlexGetCone(dm, p, &cone);
880: DMPlexGetConeSize(dm, p, &coneSize);
881: isbf = PETSC_TRUE;
882: for (i = 0; i < coneSize; i++) {
883: PetscInt v, d;
885: DMLabelGetValue(perLabel, cone[i], &v);
886: DMLabelGetDefaultValue(perLabel, &d);
887: isbf = (PetscBool)(isbf && v != d);
888: }
889: }
890: if (isbf) PetscBTSet(bfaces, p - fStart);
891: }
892: /* count boundary faces */
893: for (p = fStart, bf = 0; p < fEnd; p++) {
894: if (PetscUnlikely(PetscBTLookup(bfaces, p - fStart))) {
895: const PetscInt *support;
896: PetscInt supportSize, c;
898: DMPlexGetSupportSize(dm, p, &supportSize);
899: DMPlexGetSupport(dm, p, &support);
900: for (c = 0; c < supportSize; c++) {
901: const PetscInt *cone;
902: PetscInt cell, cl, coneSize;
904: cell = support[c];
905: if (pown && PetscUnlikely(!PetscBTLookup(pown, cell - cStart))) continue;
906: DMPlexGetCone(dm, cell, &cone);
907: DMPlexGetConeSize(dm, cell, &coneSize);
908: for (cl = 0; cl < coneSize; cl++) {
909: if (cone[cl] == p) {
910: bf += 1;
911: break;
912: }
913: }
914: }
915: }
916: }
917: minl = 1;
918: label = NULL;
919: if (enable_bmark) {
920: PetscInt lminl = PETSC_MAX_INT;
922: DMGetLabel(dm, bmark, &label);
923: if (label) {
924: IS vals;
925: PetscInt ldef;
927: DMLabelGetDefaultValue(label, &ldef);
928: DMLabelGetValueIS(label, &vals);
929: ISGetMinMax(vals, &lminl, NULL);
930: ISDestroy(&vals);
931: lminl = PetscMin(ldef, lminl);
932: }
933: MPIU_Allreduce(&lminl, &minl, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
934: if (minl == PETSC_MAX_INT) minl = 1;
935: }
936: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", bf);
937: for (p = fStart; p < fEnd; p++) {
938: if (PetscUnlikely(PetscBTLookup(bfaces, p - fStart))) {
939: const PetscInt *support;
940: PetscInt supportSize, c, nc = 0;
942: DMPlexGetSupportSize(dm, p, &supportSize);
943: DMPlexGetSupport(dm, p, &support);
944: if (pown) {
945: for (c = 0; c < supportSize; c++) {
946: if (PetscLikely(PetscBTLookup(pown, support[c] - cStart))) fcells[nc++] = support[c];
947: }
948: } else
949: for (c = 0; c < supportSize; c++) fcells[nc++] = support[c];
950: for (c = 0; c < nc; c++) {
951: const DMPolytopeType *faceTypes;
952: DMPolytopeType cellType;
953: const PetscInt *faceSizes, *cone;
954: PetscInt vids[8], *faces, st, i, coneSize, cell, cl, nv, cid = -1, mid = -1;
956: cell = fcells[c];
957: DMPlexGetCone(dm, cell, &cone);
958: DMPlexGetConeSize(dm, cell, &coneSize);
959: for (cl = 0; cl < coneSize; cl++)
960: if (cone[cl] == p) break;
961: if (cl == coneSize) continue;
963: /* face material id and type */
964: DMPlexGetPointMFEMCellID_Internal(dm, label, minl, p, &mid, &cid);
965: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT, mid, cid);
966: /* vertex ids */
967: DMPlexGetCellType(dm, cell, &cellType);
968: DMPlexGetPointMFEMVertexIDs_Internal(dm, cell, (localized && !hovec) ? coordSection : NULL, &nv, vids);
969: DMPlexGetRawFaces_Internal(dm, cellType, vids, NULL, &faceTypes, &faceSizes, (const PetscInt **)&faces);
970: st = 0;
971: for (i = 0; i < cl; i++) st += faceSizes[i];
972: DMPlexInvertCell(faceTypes[cl], faces + st);
973: for (i = 0; i < faceSizes[cl]; i++) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, faces[st + i]);
974: PetscViewerASCIIPrintf(viewer, "\n");
975: DMPlexRestoreRawFaces_Internal(dm, cellType, vids, NULL, &faceTypes, &faceSizes, (const PetscInt **)&faces);
976: bf -= 1;
977: }
978: }
979: }
981: PetscBTDestroy(&bfaces);
982: PetscFree(fcells);
983: }
985: /* mark owned vertices */
986: vown = NULL;
987: if (pown) {
988: PetscBTCreate(vEnd - vStart, &vown);
989: for (p = cStart; p < cEnd; p++) {
990: PetscInt i, closureSize, *closure = NULL;
992: if (PetscUnlikely(!PetscBTLookup(pown, p - cStart))) continue;
993: DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
994: for (i = 0; i < closureSize; i++) {
995: const PetscInt pp = closure[2 * i];
997: if (pp >= vStart && pp < vEnd) PetscBTSet(vown, pp - vStart);
998: }
999: DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
1000: }
1001: }
1003: if (parentSection) {
1004: PetscInt vp, gvp;
1006: for (vp = 0, p = vStart; p < vEnd; p++) {
1007: DMLabel dlabel;
1008: PetscInt parent, depth;
1010: if (PetscUnlikely(vown && !PetscBTLookup(vown, p - vStart))) continue;
1011: DMPlexGetDepthLabel(dm, &dlabel);
1012: DMLabelGetValue(dlabel, p, &depth);
1013: DMPlexGetTreeParent(dm, p, &parent, NULL);
1014: if (parent != p) vp++;
1015: }
1016: MPIU_Allreduce(&vp, &gvp, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm));
1017: if (gvp) {
1018: PetscInt maxsupp;
1019: PetscBool *skip = NULL;
1021: PetscViewerASCIIPrintf(viewer, "\nvertex_parents\n");
1022: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", vp);
1023: DMPlexGetMaxSizes(dm, NULL, &maxsupp);
1024: PetscMalloc1(maxsupp, &skip);
1025: for (p = vStart; p < vEnd; p++) {
1026: DMLabel dlabel;
1027: PetscInt parent;
1029: if (PetscUnlikely(vown && !PetscBTLookup(vown, p - vStart))) continue;
1030: DMPlexGetDepthLabel(dm, &dlabel);
1031: DMPlexGetTreeParent(dm, p, &parent, NULL);
1032: if (parent != p) {
1033: PetscInt vids[8] = {-1, -1, -1, -1, -1, -1, -1, -1}; /* silent overzealous clang static analyzer */
1034: PetscInt i, nv, ssize, n, numChildren, depth = -1;
1035: const PetscInt *children;
1037: DMPlexGetConeSize(dm, parent, &ssize);
1038: switch (ssize) {
1039: case 2: /* edge */
1040: nv = 0;
1041: DMPlexGetPointMFEMVertexIDs_Internal(dm, parent, localized ? coordSection : NULL, &nv, vids);
1042: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, p - vStart);
1043: for (i = 0; i < nv; i++) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, vids[i]);
1044: PetscViewerASCIIPrintf(viewer, "\n");
1045: vp--;
1046: break;
1047: case 4: /* face */
1048: DMPlexGetTreeChildren(dm, parent, &numChildren, &children);
1049: for (n = 0; n < numChildren; n++) {
1050: DMLabelGetValue(dlabel, children[n], &depth);
1051: if (!depth) {
1052: const PetscInt *hvsupp, *hesupp, *cone;
1053: PetscInt hvsuppSize, hesuppSize, coneSize;
1054: PetscInt hv = children[n], he = -1, f;
1056: PetscArrayzero(skip, maxsupp);
1057: DMPlexGetSupportSize(dm, hv, &hvsuppSize);
1058: DMPlexGetSupport(dm, hv, &hvsupp);
1059: for (i = 0; i < hvsuppSize; i++) {
1060: PetscInt ep;
1061: DMPlexGetTreeParent(dm, hvsupp[i], &ep, NULL);
1062: if (ep != hvsupp[i]) {
1063: he = hvsupp[i];
1064: } else {
1065: skip[i] = PETSC_TRUE;
1066: }
1067: }
1069: DMPlexGetCone(dm, he, &cone);
1070: vids[0] = (cone[0] == hv) ? cone[1] : cone[0];
1071: DMPlexGetSupportSize(dm, he, &hesuppSize);
1072: DMPlexGetSupport(dm, he, &hesupp);
1073: for (f = 0; f < hesuppSize; f++) {
1074: PetscInt j;
1076: DMPlexGetCone(dm, hesupp[f], &cone);
1077: DMPlexGetConeSize(dm, hesupp[f], &coneSize);
1078: for (j = 0; j < coneSize; j++) {
1079: PetscInt k;
1080: for (k = 0; k < hvsuppSize; k++) {
1081: if (hvsupp[k] == cone[j]) {
1082: skip[k] = PETSC_TRUE;
1083: break;
1084: }
1085: }
1086: }
1087: }
1088: for (i = 0; i < hvsuppSize; i++) {
1089: if (!skip[i]) {
1090: DMPlexGetCone(dm, hvsupp[i], &cone);
1091: vids[1] = (cone[0] == hv) ? cone[1] : cone[0];
1092: }
1093: }
1094: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, hv - vStart);
1095: for (i = 0; i < 2; i++) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, vids[i] - vStart);
1096: PetscViewerASCIIPrintf(viewer, "\n");
1097: vp--;
1098: }
1099: }
1100: break;
1101: default:
1102: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Don't know how to deal with support size %" PetscInt_FMT, ssize);
1103: }
1104: }
1105: }
1106: PetscFree(skip);
1107: }
1109: }
1110: PetscBTDestroy(&vown);
1112: /* vertices */
1113: if (hovec) { /* higher-order meshes */
1114: const char *fec;
1115: PetscInt i, n, s;
1116: PetscViewerASCIIPrintf(viewer, "\nvertices\n");
1117: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", vEnd - vStart);
1118: PetscViewerASCIIPrintf(viewer, "nodes\n");
1119: PetscObjectGetName((PetscObject)hovec, &fec);
1120: PetscViewerASCIIPrintf(viewer, "FiniteElementSpace\n");
1121: PetscViewerASCIIPrintf(viewer, "%s\n", fec);
1122: PetscViewerASCIIPrintf(viewer, "VDim: %" PetscInt_FMT "\n", sdim);
1123: PetscViewerASCIIPrintf(viewer, "Ordering: 1\n\n"); /*Ordering::byVDIM*/
1124: if (hoSection) {
1125: DM cdm;
1127: VecGetDM(hovec, &cdm);
1128: for (p = cStart; p < cEnd; p++) {
1129: PetscScalar *vals = NULL;
1130: PetscInt csize;
1132: if (PetscUnlikely(pown && !PetscBTLookup(pown, p - cStart))) continue;
1133: DMPlexVecGetClosure(cdm, hoSection, hovec, p, &csize, &vals);
1135: for (i = 0; i < csize / sdim; i++) {
1136: for (s = 0; s < sdim; s++) PetscViewerASCIIPrintf(viewer, fmt, (double)PetscRealPart(vals[i * sdim + s]));
1137: PetscViewerASCIIPrintf(viewer, "\n");
1138: }
1139: DMPlexVecRestoreClosure(cdm, hoSection, hovec, p, &csize, &vals);
1140: }
1141: } else {
1142: VecGetArrayRead(hovec, &array);
1143: VecGetLocalSize(hovec, &n);
1145: for (i = 0; i < n / sdim; i++) {
1146: for (s = 0; s < sdim; s++) PetscViewerASCIIPrintf(viewer, fmt, (double)PetscRealPart(array[i * sdim + s]));
1147: PetscViewerASCIIPrintf(viewer, "\n");
1148: }
1149: VecRestoreArrayRead(hovec, &array);
1150: }
1151: } else {
1152: VecGetLocalSize(coordinates, &nvert);
1153: PetscViewerASCIIPrintf(viewer, "\nvertices\n");
1154: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", nvert / sdim);
1155: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", sdim);
1156: VecGetArrayRead(coordinates, &array);
1157: for (p = 0; p < nvert / sdim; p++) {
1158: PetscInt s;
1159: for (s = 0; s < sdim; s++) {
1160: PetscReal v = PetscRealPart(array[p * sdim + s]);
1162: PetscViewerASCIIPrintf(viewer, fmt, PetscIsInfOrNanReal(v) ? 0.0 : (double)v);
1163: }
1164: PetscViewerASCIIPrintf(viewer, "\n");
1165: }
1166: VecRestoreArrayRead(coordinates, &array);
1167: }
1168: PetscBTDestroy(&pown);
1169: VecDestroy(&hovec);
1170: return 0;
1171: }
1173: PetscErrorCode DMPlexView_GLVis(DM dm, PetscViewer viewer)
1174: {
1175: DMView_GLVis(dm, viewer, DMPlexView_GLVis_ASCII);
1176: return 0;
1177: }