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